utils.c

关于有直接稀疏PCA的方法

#include "sparsesvd.h"

// BLAS wrapper for win32 version
#ifdef WIN32
void 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 linuxp
void 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 references
double 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