#include "sparsesvd.h"// BLAS wrapper for win32 version#ifdef WIN32void cblas_dscal( int N, double alpha, double *X, int incX){	dscal(&N,&alpha,X,&incX);}void cblas_dcopy(int N,double *X,int incX,double *Y,int incY){	dcopy(&N,X,&incX,Y,&incY);}void cblas_dgemm(enum CBLAS_ORDER Order,enum CBLAS_TRANSPOSE transA, enum CBLAS_TRANSPOSE transB,                 int M, int N, int K, double alpha, double *A, int lda,                 double *B, int ldb, double beta, double *C, int ldc){	char ta[1],tb[1];	if (transA==111)	{		*ta='N';	}	else	{		*ta='T';	};	if (transB==111) 	{		*tb='N';	}	else	{		*tb='T';	};	dgemm(ta,tb,&M,&N,&K,&alpha,A,&lda,B,&ldb,&beta,C,&ldc);}void cblas_dgemv(enum CBLAS_ORDER Order,enum CBLAS_TRANSPOSE transA,                 int M, int N, double alpha, double *A, int lda,                 double *B, int incB, double beta, double *C, int incC){	char ta[1];	if (transA==111)	{		*ta='N';	}	else	{		*ta='T';	};	dgemv(ta,&M,&N,&alpha,A,&lda,B,&incB,&beta,C,&incC);}void cblas_daxpy(int N,double alpha,double *X,int incX,double *Y,int incY){	daxpy(&N,&alpha,X,&incX,Y,&incY);}void cblas_dger(enum CBLAS_ORDER Order,int m,int n,double alpha,double *x,int incx,double *y,int incy,double *A,int lda){	dger(&m,&n,&alpha,x,&incx,y,&incy,A,&lda);}#endif WIN32// BLAS wrapper for the linux version#ifdef linuxpvoid cblas_dscal( int N, double alpha, double *X, int incX){	dscal(&N,&alpha,X,&incX);}void cblas_dcopy(int N,double *X,int incX,double *Y,int incY){	dcopy(&N,X,&incX,Y,&incY);}void cblas_dgemm(enum CBLAS_ORDER Order,enum CBLAS_TRANSPOSE transA, enum CBLAS_TRANSPOSE transB,                 int M, int N, int K, double alpha, double *A, int lda,                 double *B, int ldb, double beta, double *C, int ldc){	char ta[1],tb[1];	if (transA==111)	{		*ta='N';	}	else	{		*ta='T';	};	if (transB==111)	{		*tb='N';	}	else	{		*tb='T';	};	dgemm(ta,tb,&M,&N,&K,&alpha,A,&lda,B,&ldb,&beta,C,&ldc);}void cblas_dgemv(enum CBLAS_ORDER Order,enum CBLAS_TRANSPOSE transA,                 int M, int N, double alpha, double *A, int lda,                 double *B, int incB, double beta, double *C, int incC){	char ta[1];	if (transA==111)	{		*ta='N';	}	else	{		*ta='T';	};	dgemv(ta,&M,&N,&alpha,A,&lda,B,&incB,&beta,C,&incC);}void cblas_daxpy(int N,double alpha,double *X,int incX,double *Y,int incY){	daxpy(&N,&alpha,X,&incX,Y,&incY);}void cblas_dger(enum CBLAS_ORDER Order,int m,int n,double alpha,double *x,int incx,double *y,int incy,double *A,int lda){	dger(&m,&n,&alpha,x,&incx,y,&incy,A,&lda);}#endif// Some useful functions ...double doubsum(double *xmat, int n){	int i;	double res=0.0;	for (i=0;i<n;i++){res+=xmat[i];};	return res;}double doubdot(double *xvec, double *yvec, int n){	int i;	double res=0.0;	for (i=0;i<n;i++){res+=xvec[i]*yvec[i];};	return res;}int idxmax(double *xmat, int n){	int i;	int res=0;	for (i=0;i<n;i++)	{		if (xmat[i]>xmat[res]) {res=i;}	}	return res;}double doubasum(double *xmat, int n){	int i;	double res=0.0;	for (i=0;i<n;i++){res+=dabsf(xmat[i]);};	return res;}double doubnorm2(double *xmat, int n){	int i;	double res=0.0;	for (i=0;i<n;i++){res+=xmat[i]*xmat[i];};	return sqrt(res);}double infnorm(double *xmat, int n){	int i,j;	double res=0.0,sum;		for (j=0;j<n;j++){		sum=0.0;		for(i=0;i<n;i++)			sum+=dabsf(xmat[j+i*n]);		if(sum>=res) res=sum;	}		return res;}double frobnorm(double *xmat, int n){	int i,j;	double res=0.0;		for (i=0;i<n;i++)		for (j=0;j<n;j++)			res+=(xmat[i*n+j]*xmat[i*n+j]);		return pow(res,.5);}double dsignf(double x){	if (x>=0)		return 1.0;	else		return -1.0;}double dminif(double x, double y){	if (x>=y)		return y;	else		return x;}double dmaxf(double x, double y){	if (x>=y)		return x;	else		return y;}int imaxf(int x, int y){	if (x>=y)		return x;	else		return y;}double dabsf(double x){	if (x>=0)		return x;	else		return -x;}void dispmat(double *xmat, int n, int m){	int i,j;		for (i=0; i<n; i++)	{		for (j=0;j<m;j++)		{			printf("%+.4f ",xmat[j*n+i]);		}		printf("\n");	}	printf("\n");}// do partial eig approximation of exp(bufmata)// return fmu, and get dmax and numeigs from parameter referencesdouble partial_eig(int n,int k,double mu,double eigcut,double *bufmata,				   double *bufmatb,double *numeigs_matlab,double *evector_temp,				   double *evector_store,double *eig,double *Dvec,double *gvec,				   double *hvec,double *Vmat,double *Umat,double *workvec,int *count,				   int addeigs, double perceigs,int check_for_more_eigs, int *arcount){	int numeigs=(int)(*numeigs_matlab),nvls=0,h,i,incx=1,n2=n*n;	int lwork,inflapack,indmax,check_other_eigs=0,neceigs=0;	double alpha,beta,hs=0.0,dmax=0.0,fmu,buf,bufmata_shift=0.0;	double *evector_index;	char jobz[1],uplo[1];		double sum_sq_eigs=0.0; // variable for eigcut check, hs is the sum of the eigs	double l2normbound,tol,minDvec;	// Arpack parameters	char which[2]="LA"; // Arpack: we want largest algebraic eigs...	int ncv,info,nconv,nummatvec;	int maxitr=500; 	
	// TODO: find the optimal value for ncv	if(numeigs<n-2 && (numeigs*1.0/n)<perceigs) {  		// skip all this if we already know we want many eigs		if (k==0 || numeigs>1 || (numeigs==1 && k%check_for_more_eigs==0)) check_other_eigs=1;		bufmata_shift=frobnorm(bufmata,n);// Simple bound on largest magnitude eigenvalue