NICE_VisLearning/math/algebra/EigValues.cpp

/** 
 * @file EigValues.cpp
 * @brief EigValues Class
 * @author Michael Koch
 * @date 08/19/2008
*/
     
#include <iostream>
     
#include "EigValues.h"

using namespace OBJREC;
#define DEBUG_ARNOLDI
// #undef DEBUG_ARNOLDI
     
using namespace std;
// refactor-nice.pl: check this substitution
// old: using namespace ice;
using namespace NICE;
     
void EVArnoldi::getEigenvalues(const GenericMatrix &data,Vector &eigenvalues,Matrix &eigenvectors,uint k)
{
          
    if(data.rows()!=data.cols())
    {
	throw("EVArnoldi: matrix has to be quadratic");
    }
          
    if(k<=0)
    {
       throw("EVArnoldi: please use k>0.");
    }
#ifdef DEBUG_ARNOLDI
    fprintf(stderr,"Initialize Matrices\n");
#endif
    uint n=data.cols();
    NICE::Matrix rmatrix=Matrix(n,k,0);//=eigenvectors
    NICE::Matrix qmatrix=Matrix(n,k,0);//saves data =eigenvectors
    eigenvalues.resize(k);
    NICE::Vector q=Vector(n);
    NICE::Vector r=Vector(n);
   
#ifdef DEBUG_ARNOLDI
    fprintf(stderr,"Random Initialization\n");
#endif

    //random initialisation
    for(uint i=0;i<k;i++)
	for(uint j=0;j<n;j++)
	     rmatrix(j, i)=drand48();
   
    //reduceddim
    double delta=1.0;
    uint iteration=0;
    while (delta>mindelta &&iteration<maxiterations)
    {
	NICE::Vector rold (rmatrix.getColumn(k-1));
	for(uint reduceddim=0;reduceddim<k;reduceddim++)
	{
	     //get Vector r_k of Matrix
	     q=rmatrix.getColumn(reduceddim);
	     q.normalizeL2();
	     // q=r_k/||r_k||
	     qmatrix.getColumnRef(reduceddim) = q;
	     
	     Vector currentCol = rmatrix.getColumnRef(reduceddim); // FIXME
	     data.multiply(currentCol,q); //r_{k+1}=A*q
	     
	     for(uint j=0;j<reduceddim;j++)
		  rmatrix.getColumnRef(reduceddim) -= qmatrix.getColumn(j)*(qmatrix.getColumn(j).scalarProduct(rmatrix.getColumn(reduceddim)));
	}
	//convergence stuff
	NICE::Vector diff=rold-rmatrix.getColumn(k-1);
	delta=diff.normL2();
	iteration++;
#ifdef DEBUG_ARNOLDI
	fprintf (stderr, "EVArnoldi: [%d] delta=%f\n", iteration, delta );
#endif
    }
    eigenvectors=rmatrix;
  
    for(uint i=0;i<k;i++)
    {
	NICE::Vector tmp;
	eigenvectors.getColumnRef(i).normalizeL2();
	data.multiply(tmp,eigenvectors.getColumn(i));
	eigenvalues[i]=tmp.scalarProduct( eigenvectors.getColumn(i) );
    }
}