#include <stdio.h>#include <math.h>#include <time.h>#include <mex.h>#include <matrix.h>#include <string.h>#ifdef WIN32#include <malloc.h>extern void dsaupd();extern void dseupd();#endif#ifdef mac#include "/System/Library/Frameworks/vecLib.framework/Headers/cblas.h"#endif// Calling LAPACK on the mac#ifdef macextern void dsaupd_();extern void dseupd_();#define dsyev dsyev_ #define dgesvd dgesvd_#define dsaupd dsaupd_#define dseupd dseupd_#endif// Calling LAPACK on linux#ifdef linuxpextern void dsaupd_();extern void dseupd_();#define dsyev dsyev_#define dgesvd dgesvd_#define dsaupd dsaupd_#define dseupd dseupd_#endif#ifdef WIN32enum CBLAS_ORDER 	{CblasRowMajor=101, CblasColMajor=102};enum CBLAS_TRANSPOSE 	{CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};enum CBLAS_UPLO		{CblasUpper=121, CblasLower=122};enum CBLAS_DIAG		{CblasNonUnit=131, CblasUnit=132};enum CBLAS_SIDE		{CblasLeft=141, CblasRight=142};void cblas_dscal(int N,double alpha, double *X,int incX);void cblas_dcopy(int N,double *X,int incX,double *Y,int 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);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);void cblas_daxpy(int N,double alpha,double *X,int incX,double *Y,int 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);#endif// Our cblas implementation on linux#ifdef linuxpenum CBLAS_ORDER 	{CblasRowMajor=101, CblasColMajor=102};enum CBLAS_TRANSPOSE 	{CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};enum CBLAS_UPLO		{CblasUpper=121, CblasLower=122};enum CBLAS_DIAG		{CblasNonUnit=131, CblasUnit=132};enum CBLAS_SIDE		{CblasLeft=141, CblasRight=142};void cblas_dscal(int N,double alpha, double *X,int incX);void cblas_dcopy(int N,double *X,int incX,double *Y,int 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);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);void cblas_daxpy(int N,double alpha,double *X,int incX,double *Y,int incY);#endif// Local functionsdouble doubsum(double *xmat, int n);double doubdot(double *xvec, double *yvec, int n);double doubasum(double *xmat, int n);double doubnorm2(double *xmat, int n);double infnorm(double *xmat, int n);double frobnorm(double *xmat, int n);int idxmax(double *xmat, int n);void sparse_rank_one(double *Amat, int n, double rho, double tol, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, int info, int checkgap, double *dualitygap_alliter, double *cputime_alliter);void sparse_geneig(double *Amat, double *Bmat, double *Rmat, int n, int m, double rho, double tol, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, int info, int checkgap, double *dualitygap_alliter, double *cputime_alliter);void sparse_rank_one_partialeig(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, int info, int numeigs, int addeigs, int checkgap, double perceigs, int check_for_more_eigs, double *dualitygap_alliter, double *cputime_alliter, double *perceigs_alliter);void sparsegeneig_partial(double *Amat, double *Bmat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, int info, int numeigs, int addeigs, int checkgap, double perceigs, int check_for_more_eigs, double *dualitygap_alliter, double *cputime_alliter, double *perceigs_alliter);void sparse_rank_one_partialeig_matlab(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iters, int info, int numeigs, int addeigs, int checkgap, double perceigs, int check_for_more_eigs, double *dualitygap_alliter, double *cputime_alliter, double *perceigs_alliter);double dsignf(double x);double dminif(double x, double y);double dmaxf(double x, double y);int imaxf(int x, int y);double dabsf(double x);void dispmat(double *xmat, int n, int m);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);double partial_eig_matlab(int n,int k,double mu,double eigcut,double *bufmata,				double *bufmatb,double *numeigs_matlab,double *evector_temp,				double *evector_store,double *evalue,mxArray *input[4],mxArray *output[3],				double *hvec,double *Vmat,double *Umat,double *workvec,int *count,				double last_dmax,int addeigs, double perceigs,int check_for_more_eigs);void symmetrize(double *xmat,double *ymat,int n);double maxeig(double *xmat,double *bufveca,double *bufvecb,int n);int simarpack(double *Amat, int n, int nev, int ncv, double tol, char* which, int maxitr, int printlevel, double* evals, double* evecs, int* nconv, int* nummatvec, int* count);void dispmat(double *xmat, int n, int m);#ifdef WIN32// BLAS functions usedvoid dscal(int *n,double *a,double *x,int *incx);void dcopy(int *n,double *x,int *incx,double *y,int *incy);void dgemm(char *transa,char *transb,int *m,int *n,int *k,double *alpha,double *a,int *lda,double *b,int *ldb,double *beta,double *c,int *ldc);void dgemv(char *trans,int *m,int *n,double *alpha,double *a,int *lda,double *x,int *incx,double *beta,double *y,int *incy);void daxpy(int *n,double *alpha,double *x,int *incx,double *y,int *incy);void dger(int *m,int *n,double *alpha,double *x,int *incx,double *y,int *incy,double *A,int *lda);// LAPACK functions usedvoid dsyev(char *jobz,char *uplo,int *n,double *a,int *lda,double *w,double *work,int *lwork,int *info);void dsyevx(char *jobz,char *range,char *uplo,int *n,double *a,int *lda,double *vl,double *vu,int *il,int *iu,double *abstol,int *m,double *w,double *z,int *ldz,double *work,int *lwork,int *iwork,int *ifail,int *info);#endif