src/lambda2.h

    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});
    }

    Usage

    Examples