NICE_VisLearning/math/algebra_trlan/EigValuesTRLAN.cpp

#include "EigValuesTRLAN.h"

#include <iostream>
#include <assert.h>
#define DEBUG_TRLAN
// #undef DEBUG_TRLAN
using namespace OBJREC;

using namespace std;

static const GenericMatrix* mat = NULL;

static void op(int* nrow, int* ncol, double* xin, int* ldx, double* xout, int* ldy)
{
    NICE::Vector xin_v(xin, *ldy, NICE::VectorBase::external);
    NICE::Vector xout_v(xout, *ldy, NICE::VectorBase::external);
    mat->multiply(xout_v, xin_v);
}

typedef void (*OpFunc)(int* nrow, int* ncol, double* xin, int* ldx, double* xout, int* ldy);
extern "C" void trlan77_(OpFunc op, int ipar[32], int* nrowl, int* mev, double* eval,
                             double* evec, int* nrow2, double* wrk, int* lwrk);

EigValuesTRLAN::EigValuesTRLAN(int magnitude,
                               int restarting_scheme,
                               double tolerance,
                               int max_mult,
                               int mult_flops)
{
    this->magnitude = magnitude;
    this->restarting_scheme = restarting_scheme;
    this->tolerance = tolerance;
    this->max_mult = max_mult;
    this->mult_flops = mult_flops;
}

// refactor-nice.pl: check this substitution
// old: void EigValuesTRLAN::getEigenvalues ( const GenericMatrix &data, Vector & eigenvalues, Matrix & eigenvectors, uint k )
void EigValuesTRLAN::getEigenvalues(const GenericMatrix &data, NICE::Vector & eigenvalues, NICE::Matrix & eigenvectors, uint k)
{
#ifdef DEBUG_TRLAN
    fprintf(stderr, "Starting eigenvalue computation with Lanczos Iteration ...\n");
#endif
    assert(data.rows() == data.cols());
    mat = &data;
    int mat_size = (int)data.rows();

    if (max_mult <= 0) max_mult = 10 * mat_size;
    if (mult_flops <= 0) mult_flops = 10 * mat_size;
    
    int computed_eigenpairs = k;
    int ipar[32] = {0, magnitude, computed_eigenpairs, 0, 0, restarting_scheme, max_mult, 0, 1, 99, -1, 98, mult_flops};
    double* evalbuff = new double[k];
    memset(evalbuff, 0, k*sizeof(double));
    double* evecbuff = new double[k*mat_size];
    memset(evecbuff, 1, k*mat_size*sizeof(double));
    int ritz_size = (computed_eigenpairs < 3 ? 2 * computed_eigenpairs : computed_eigenpairs + 6);
    int workspace_size = ritz_size * (ritz_size + 10);
    double* workspace = new double[workspace_size];
    workspace[0] = tolerance;
#ifdef DEBUG_TRLAN
    fprintf(stderr, "Initialization ready ... Starting TRLAN\n");
#endif
    // call lanczos iteration of trlan
    int kint = k;
    trlan77_(op, ipar, &mat_size, &kint, evalbuff, evecbuff, &mat_size, workspace, &workspace_size);

    if(ipar[0]!=0)
    {
    	fprintf(stderr,"EigValuesTRLAN: Error occured by calling TRLAN Eigenvalue Computation.\n");
    	fprintf(stderr,"EigValuesTRLAN: %d converged eigenvectors \n",ipar[3]);

    }
#ifdef DEBUG_TRLAN
    fprintf(stderr, "TRLAN ready\n");
#endif

    //copy eigenvectors in evecbuff to stl vector
    eigenvectors.resize(mat_size, k);
    eigenvalues.resize(k);
#ifdef DEBUG_TRLAN
    fprintf(stderr, "Store results\n");
#endif
    if (magnitude < 0)
    {
        for (uint i = 0;i < k;i++)
            for (uint j = 0;j < (uint)mat_size;j++)
                eigenvectors(j, i) = evecbuff[i*mat_size+j];
        for (uint i = 0;i < k;i++)
            eigenvalues[i] = evalbuff[i];
    }
    else
    {
        for (uint i = 0;i < k;i++)
            for (uint j = 0;j < (uint)mat_size;j++)
                eigenvectors(j, k - i - 1) = evecbuff[i*mat_size+j];
        for (uint i = 0;i < k;i++)
            eigenvalues[i] = evalbuff[k-i-1];

    }
#ifdef DEBUG_TRLAN
    fprintf(stderr, "Clean up\n");
#endif

    delete[] workspace;
    delete[] evalbuff;
    delete[] evecbuff;
#ifdef DEBUG_TRLAN
    fprintf(stderr, "Clean up Done\n");
#endif
}