Subversion Repositories gelsvn

#include "eigensolution.h"
#include "Mat2x2f.h"
#include "Mat3x3f.h"
#include "Mat4x4f.h"

#include <iostream>

using namespace std;


namespace
{
  // During experiments 925 iterations were observed for a hard problem
  // Ten times that should be safe.
  const int KMAX = 10000;

  // The threshold below is the smallest that seems to give reliable
  // solutions.
  const float EV_THRESH = 0.0000001;
}

namespace CGLA
{

  template <class MT>
  int power_eigensolution(const MT& Ap, MT& Q, MT& L, int max_sol)
  {
    L = MT(0);

    typedef typename MT::VectorType VT;
    MT A = Ap;
    int n = s_min(MT::get_v_dim(), max_sol);

    for(int i=0;i<n;++i)
      {
        // Seed the eigenvector estimate
        VT q(1);
        float l=0,l_old;

        // As long as we haven't reached the max iterations and the
        // eigenvalue has not converged, do
        int k=0;
        for(; k<2 || k<KMAX && (l-l_old > EV_THRESH * l) ; ++k)
          {
            // Multiply the eigenvector estimate onto A
            const VT z = A * q;

            // Check that we did not get the null vector.
            float z_len = z.length();
            if(z_len < EV_THRESH)
              return i;

            // Normalize to get the new eigenvector
            q = z/z_len;

            // Record the old eigenvalue estimate and get a new estimate.
            l_old = l;
            l = dot(q, A * q);
          }
        // If we hit the max iterations, we also don't trust the eigensolution
        if(k==KMAX)
          return i;
                                
        // Update the solution by adding the eigenvector to Q and
        // the eigenvalue to the diagonal of L.
        Q[i] = q;
        L[i][i] = l;

        // Update A by subtracting the subspace represented by the 
        // eigensolution just found. This is called the method of 
        // deflation.
        MT B;
        outer_product(q,q,B);
        A = A - l * B;
      }
    return n;
  }

  /* There is no reason to put this template in a header file, since 
     we will only use it on matrices defined in CGLA. Instead, we 
     explicitly instantiate the function for the square matrices
     of CGLA */
  template int power_eigensolution<Mat2x2f>(const Mat2x2f&,
                                            Mat2x2f&,Mat2x2f&,int);
        
  template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
                                            Mat3x3f&,Mat3x3f&,int);
  template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
                                            Mat4x4f&,Mat4x4f&,int);
}