tensor-eig3.cc

Go to the documentation of this file.

// http://barnesc.blogspot.com/2007/02/eigenvectors-of-3x3-symmetric-matrix.html
// http://www.connellybarnes.com/code/c/eig3-1.0.0.zip
// A C++ source and header file to compute eigenvectors/values of a 3x3 symmetric
#include "rheolef/tensor.h"
namespace rheolef { 

template<class T>
static
inline
T hypot2(const T& x, const T& y) { return sqrt(x*x+y*y); }

template<class T>
static
void
tred2 (tensor_basic<T>& V, point_basic<T>& d, point_basic<T>& e)
{
  const int n = 3;


  for (int j = 0; j < n; j++) {
    d[j] = V(n-1,j);
  }


  for (int i = n-1; i > 0; i--) {

    // Scale to avoid under/overflow.

    T scale = 0.0;
    T h = 0.0;
    for (int k = 0; k < i; k++) {
      scale = scale + fabs(d[k]);
    }
    if (scale == 0.0) {
      e[i] = d[i-1];
      for (int j = 0; j < i; j++) {
        d[j] = V(i-1,j);
        V(i,j) = 0.0;
        V(j,i) = 0.0;
      }
    } else {


      for (int k = 0; k < i; k++) {
        d[k] /= scale;
        h += d[k] * d[k];
      }
      T f = d[i-1];
      T g = sqrt(h);
      if (f > 0) {
        g = -g;
      }
      e[i] = scale * g;
      h = h - f * g;
      d[i-1] = f - g;
      for (int j = 0; j < i; j++) {
        e[j] = 0.0;
      }


      for (int j = 0; j < i; j++) {
        f = d[j];
        V(j,i) = f;
        g = e[j] + V(j,j) * f;
        for (int k = j+1; k <= i-1; k++) {
          g += V(k,j) * d[k];
          e[k] += V(k,j) * f;
        }
        e[j] = g;
      }
      f = 0.0;
      for (int j = 0; j < i; j++) {
        e[j] /= h;
        f += e[j] * d[j];
      }
      T hh = f / (h + h);
      for (int j = 0; j < i; j++) {
        e[j] -= hh * d[j];
      }
      for (int j = 0; j < i; j++) {
        f = d[j];
        g = e[j];
        for (int k = j; k <= i-1; k++) {
          V(k,j) -= (f * e[k] + g * d[k]);
        }
        d[j] = V(i-1,j);
        V(i,j) = 0.0;
      }
    }
    d[i] = h;
  }


  for (int i = 0; i < n-1; i++) {
    V(n-1,i) = V(i,i);
    V(i,i) = 1.0;
    T h = d[i+1];
    if (h != 0.0) {
      for (int k = 0; k <= i; k++) {
        d[k] = V(k,i+1) / h;
      }
      for (int j = 0; j <= i; j++) {
        T g = 0.0;
        for (int k = 0; k <= i; k++) {
          g += V(k,i+1) * V(k,j);
        }
        for (int k = 0; k <= i; k++) {
          V(k,j) -= g * d[k];
        }
      }
    }
    for (int k = 0; k <= i; k++) {
      V(k,i+1) = 0.0;
    }
  }
  for (int j = 0; j < n; j++) {
    d[j] = V(n-1,j);
    V(n-1,j) = 0.0;
  }
  V(n-1,n-1) = 1.0;
  e[0] = 0.0;
} 
template<class T>
static
void
tql2 (tensor_basic<T>& V, point_basic<T>& d, point_basic<T>& e)
{
  const int n = 3;


  for (int i = 1; i < n; i++) {
    e[i-1] = e[i];
  }
  e[n-1] = 0.0;

  T f = 0.0;
  T tst1 = 0.0;
  T eps = pow(2.0,-52.0);
  for (int l = 0; l < n; l++) {


    tst1 = std::max (tst1, fabs(d[l]) + fabs(e[l]));
    int m = l;
    while (m < n) {
      if (fabs(e[m]) <= eps*tst1) {
        break;
      }
      m++;
    }


    if (m > l) {
      int iter = 0;
      do {
        iter = iter + 1;  // (Could check iteration count here.)


        T g = d[l];
        T p = (d[l+1] - g) / (2.0 * e[l]);
        T r = hypot2(p,1.0);
        if (p < 0) {
          r = -r;
        }
        d[l] = e[l] / (p + r);
        d[l+1] = e[l] * (p + r);
        T dl1 = d[l+1];
        T h = g - d[l];
        for (int i = l+2; i < n; i++) {
          d[i] -= h;
        }
        f = f + h;


        p = d[m];
        T c = 1.0;
        T c2 = c;
        T c3 = c;
        T el1 = e[l+1];
        T s = 0.0;
        T s2 = 0.0;
        for (int i = m-1; i >= l; i--) {
          c3 = c2;
          c2 = c;
          s2 = s;
          g = c * e[i];
          h = c * p;
          r = hypot2(p,e[i]);
          e[i+1] = s * r;
          s = e[i] / r;
          c = p / r;
          p = c * d[i] - s * g;
          d[i+1] = h + s * (c * g + s * d[i]);


          for (int k = 0; k < n; k++) {
            h = V(k,i+1);
            V(k,i+1) = s * V(k,i) + c * h;
            V(k,i) = c * V(k,i) - s * h;
          }
        }
        p = -s * s2 * c3 * el1 * e[l] / dl1;
        e[l] = s * p;
        d[l] = c * p;


      } while (fabs(e[l]) > eps*tst1);
    }
    d[l] = d[l] + f;
    e[l] = 0.0;
  }
  

  for (int i = 0; i < n-1; i++) {
    int k = i;
    T p = d[i];
    for (int j = i+1; j < n; j++) {
      if (d[j] > p) {
        k = j;
        p = d[j];
      }
    }
    if (k != i) {
      d[k] = d[i];
      d[i] = p;
      for (int j = 0; j < n; j++) {
        p = V(j,i);
        V(j,i) = V(j,k);
        V(j,k) = p;
      }
    }
  }
}
template<class T>
void
eigen_decomposition (const tensor_basic<T>& A, tensor_basic<T>& V, point_basic<T>& d)
{
  V = A;
  point_basic<T> e;
  tred2 (V, d, e);
  tql2  (V, d, e);
}
#define _RHEOLEF_instanciation(T)                                               \
template void eigen_decomposition (const tensor_basic<T>& A, tensor_basic<T>& V, point_basic<T>& d);

_RHEOLEF_instanciation(Float)

}// namespace rheolef