static void eigsrt (double d[dimension], double v[dimension][dimension]) { int k, j, i; double p; for (i = 0; i < dimension - 1; i++) { p = d[k = i]; for (j = i + 1; j < dimension; j++) if (d[j] >= p) p = d[k = j]; if (k != i) { d[k] = d[i]; d[i] = p; for (j = 0; j < dimension; j++) { p = v[j][i]; v[j][i] = v[j][k]; v[j][k] = p; } } } } #define ROTATE(a,i,j,k,l) {\ g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau);} /** * eigenvalues: * @a: a symmetric matrix. * @d: a vector. * @v: another matrix. * * Fills @d (resp. @v) with the eigenvalues (resp. eigenvectors) of * matrix @a. */ void eigenvalues (double a[dimension][dimension], double d[dimension], double v[dimension][dimension]) { int j, iq, ip, i; double tresh, theta, tau, t, sm, s, h, g, c, b[dimension], z[dimension]; for (ip = 0; ip < dimension; ip++) { for (iq = 0; iq < dimension; iq++) v[ip][iq] = 0.0; v[ip][ip] = 1.0; } for (ip = 0; ip < dimension; ip++) { b[ip] = d[ip] = a[ip][ip]; z[ip] = 0.0; } for (i = 1; i <= 50; i++) { sm = 0.0; for (ip = 0; ip < dimension - 1; ip++) { for (iq = ip + 1; iq < dimension; iq++) sm += fabs (a[ip][iq]); } if (sm == 0.0) { eigsrt (d, v); return; } if (i < 4) tresh = 0.2*sm/(dimension*dimension); else tresh = 0.0; for (ip = 0; ip < dimension - 1; ip++) { for (iq = ip + 1; iq < dimension; iq++) { g = 100.0*fabs (a[ip][iq]); if (i > 4 && fabs(d[ip]) + g == fabs(d[ip]) && fabs(d[iq]) + g == fabs(d[iq])) a[ip][iq] = 0.0; else if (fabs (a[ip][iq]) > tresh) { h = d[iq] - d[ip]; if (fabs(h) + g == fabs(h)) t = a[ip][iq]/h; else { theta = 0.5*h/a[ip][iq]; t = 1.0/(fabs (theta) + sqrt (1.0 + theta*theta)); if (theta < 0.0) t = -t; } c = 1.0/sqrt (1 + t*t); s = t*c; tau = s/(1.0 + c); h = t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq] = 0.0; for (j = 0; j <= ip - 1; j++) ROTATE (a, j, ip, j, iq); for (j = ip + 1; j <= iq - 1; j++) ROTATE (a, ip, j, j, iq); for (j = iq + 1; j < dimension; j++) ROTATE(a, ip, j, iq, j); for (j = 0; j < dimension; j++) ROTATE(v, j, ip, j, iq); } } } for (ip = 0; ip < dimension; ip++) { b[ip] += z[ip]; d[ip] = b[ip]; z[ip] = 0.0; } } /* Too many iterations */ for (i = 0; i < dimension; i++) { for (j = 0; j < dimension; j++) fprintf (stderr, "%10.3g ", a[i][j]); fprintf (stderr, "\n"); } assert (false); } void lambda2 (const vector u, scalar l2) { foreach() { double JJ[dimension][dimension]; scalar s = u.x; int i = 0; foreach_dimension() JJ[0][i++] = center_gradient (s); s = u.y; i = 0; foreach_dimension() JJ[1][i++] = center_gradient (s); s = u.z; i = 0; foreach_dimension() JJ[2][i++] = center_gradient (s); double S2O2[dimension][dimension]; for (int i = 0; i < dimension; i++) for (int j = 0; j < dimension; j++) { S2O2[i][j] = 0.; for (int k = 0; k < dimension; k++) S2O2[i][j] += JJ[i][k]*JJ[k][j] + JJ[k][i]*JJ[j][k]; } double lambda[dimension], ev[dimension][dimension]; eigenvalues (S2O2, lambda, ev); l2[] = lambda[1]/2.; } boundary ({l2}); }