lambda2.h

Go to the documentation of this file.

/** @file 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)
{
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double JJ[dimension][dimension];
    scalar s = u.x;
    int i = 0;
    for (int _d = 0; _d < dimension; _d++)
      JJ[0][i++] = center_gradient (s);
    s = u.y; i = 0;
    for (int _d = 0; _d < dimension; _d++)
      JJ[1][i++] = center_gradient (s);
    s = u.z; i = 0;
    for (int _d = 0; _d < dimension; _d++)
      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.;
  }
}