#include "Eigenproblem.h"

#include <iostream>

using namespace OBJREC;

using namespace std;

	static Eigenproblem* eigen=NULL;

	static void op(int* nrow, int* ncol, double* xin, int* ldx, double* xout, int* ldy){
		eigen->sparseMult(xin,xout,*ldy);
	}

	void Eigenproblem::sparseMult(double* xin, double* xout,int ldy){
		memset(xout,0,ldy*sizeof(double));
		for(uint i=0;i<mat.size();i++){
			xout[mat[i].x]+=mat[i].val * xin[mat[i].y];
		}
	}


	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);	

	void Eigenproblem::solve(vector<element> mat, int mat_size, vector<double> &evals, 
                             vector< vector<double> > &evecs, int k, int magnitude,
                             int restarting_scheme, double tolerance, int max_mult,
                             int mult_flops){
	
		this->mat=mat;

		if(max_mult<=0) max_mult=10*mat_size;
		if(mult_flops<=0) mult_flops=10*mat_size;

		eigen=this;

		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;
		trlan77_(op,ipar,&mat_size,&k,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]);

		    }
		//copy eigenvectors in evecbuff to stl vector
    		for(int i=0;i<k;i++){
			vector<double> vec;
			for(int j=0;j<mat_size;j++){
				vec.push_back(evecbuff[i*mat_size+j]);
			}
			evecs.push_back(vec);
		}
		for(int i=0;i<k;i++){
			evals.push_back(evalbuff[i]);
		}
		delete[] evalbuff;
		delete[] evecbuff;           
	}