/* Finds a sparse rank-one approximation to a given symmetric matrix A, by solving the SDP	min_X lambda_max(A+X) : X = X', abs(X(i,j)) <= rho, 1<=i,j<= nand its dual:	max_U Tr(UA) - rho sum_ij |U_ij| : U=U', U \succeq 0, Tr(U)=1***	inputs: ***A			nxn symmetric matrix (left unchanged)n			problem sizerho			non-negative scalar gapchange	required change in gap from first gap (default: 1e-4) MaxIter		maximum number of iterationsinfo		controls verbosity: 0 silent, n>0 frequency of progress reportWarmStart	0 if cold start, k0 if WarmStart (total number of iterations in previous run)F			Average gradient (for warm start, Fmat is updated)***	outputs: ***X			solves the primal SDP U			dual variable, solves the dual SDP u			largest eigenvector of U F			Average gradientThis code implements Nesterov's smooth minimization algorithm. See: Y. Nesterov "Smooth Minimization of NonSmooth Functions."Here, the gradient is only only computed aproximately. See A. d'Aspremont "Smooth optimization with approximate gradient."Last Modified: Alexandre d'Aspremont, Laurent El Ghaoui, Vijay Krishnamurthy, Ronny Luss November 2007.http://www.carva.org/alexandre.daspremont*/ #include "sparsesvd.h"// ****** UPDATE CODE COMMENTS ******************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){	// Hard parameters	int Nperiod=imaxf(1,info);	int work_size=3*n+n*n;	// Working variables	double norma12,mu,L;	double alpha,beta,buf,gapk;	double dmax=0.0,fmu,lambda;	int n2=n*n,m2=m*m,incx=1,precision_flag=0,iteration_flag=0;	int lwork=work_size,inflapack,indmax,k=0,i,j;	char jobz[1]="V",uplo[1]="U";	double cputime,last_time=(double)clock();double start_time=(double)clock();int left_h=0,left_m=0,left_s=0;	double *Vmat=(double *) calloc(m*m,sizeof(double));	double *bufmata=(double *) calloc(m*m,sizeof(double));	double *bufmatb=(double *) calloc(m*m,sizeof(double));	double *bufmatr=(double *) calloc(n*m,sizeof(double));	double *bufmatrb=(double *) calloc(n*m,sizeof(double));	double *avec=(double *) calloc(m,sizeof(double));	double *Dvec=(double *) calloc(n,sizeof(double));	double *qvec=(double *) calloc(n,sizeof(double));	double *workvec=(double *) calloc(work_size,sizeof(double));	double *gvec=(double *) calloc(n,sizeof(double));	double *hvec=(double *) calloc(n,sizeof(double));	int checkgap_count=0,firstiter=0; // added for test variables		// Start...	if (info>=1){		mexPrintf("DSPCA starting... Sparse generalized eig. maximization.\n");mexEvalString("drawnow;");}	// Test malloc results	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL) || (Dvec==NULL) || (workvec==NULL) || (gvec==NULL) || (hvec==NULL)){		mexPrintf("DSPCA: memory allocation failed ... \n");mexEvalString("drawnow;");return;}	// Replace B by B^{-1/2}. DEBUG, this saves space but changes the value of B in the code, dangerous !!!!!!!	dsyev(jobz,uplo,&n,Bmat,&n,qvec,workvec,&lwork,&inflapack);	alpha=0.0;cblas_dscal(n2,alpha,bufmatb,incx);	for (i=0;i<n;i++) {bufmatb[i*n+i]=1.0/sqrt(qvec[i]);}	alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Bmat,n,bufmatb,n,beta,bufmata,n);	cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,n,n,n,alpha,bufmata,n,Bmat,n,beta,bufmatb,n);		cblas_dcopy(n2,bufmatb,incx,Bmat,incx);	// Compute alpha vector	alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,m,n,n,alpha,Rmat,m,Bmat,n,beta,bufmatr,m);	for(i=0;i<m;i++){		avec[i]=0;		for(j=0;j<n;j++){avec[i]+=bufmatr[m*j+i]*bufmatr[m*j+i];}		avec[i]=sqrt(avec[i]);}	// Compute intial complexity params	norma12=0.0;for(i=0;i<n;i++){norma12+=1.0/qvec[i];};norma12=sqrt(norma12);	L=0.001*(2.0*log(n)*norma12*norma12)/tol;cputime=start_time;mu=tol/(2.0*log(n));	alpha=0.0;cblas_dscal(m2,alpha,Xmat,incx);	while ((precision_flag+iteration_flag)==0){ 		// eigenvalue decomposition of M*(A+R'*X*R)*M		cblas_dcopy(m2,Xmat,incx,Vmat,incx);		alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,n,m,m,alpha,Rmat,m,Vmat,m,beta,bufmatr,n);		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,m,alpha,bufmatr,n,Rmat,m,beta,Vmat,n);		alpha=1.0;cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);		alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Bmat,n,Vmat,n,beta,bufmata,n);		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,bufmata,n,Bmat,n,beta,Vmat,n);		dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack); // call LAPACK (most CPU time is here)		// compute fmu(X) = mu*log(trace((exp(M*(A+X)*M)/mu)))-mu*log(n) reliably 		indmax=idxmax(Dvec,n);dmax=Dvec[indmax];		for (i=0;i<n;i++) {hvec[i]=exp((Dvec[i]-dmax)/mu);}		buf=doubsum(hvec,n);fmu=dmax+mu*log(buf/n);		// compute gradient of fmu w.r.t. X, which is the dual variable U 		alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Bmat,n,Vmat,n,beta,bufmatb,n);		alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,m,n,n,alpha,Rmat,m,bufmatb,n,beta,bufmatr,m);		alpha=0.0;cblas_dscal(n2,alpha,Vmat,incx);for(i=0;i<n;i++){gvec[i]=hvec[i]/buf;Vmat[i*n+i]=gvec[i];}		alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,m,n,n,alpha,bufmatr,m,Vmat,n,beta,bufmatrb,m);		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,m,m,n,alpha,bufmatrb,m,bufmatr,m,beta,Umat,m);		// update gradient's weighted average 		alpha=((double)(k)+1.0)/2.0;		cblas_daxpy(m2,alpha,Umat,incx,Fmat,incx);		// find a projection of X-Gmu/L on feasible set 		cblas_dcopy(m2,Xmat,incx,bufmata,incx);		alpha=-1.0/L;		cblas_daxpy(m2,alpha,Umat,incx,bufmata,incx);		// project again		alpha=-(1.0/L);		cblas_dcopy(m2,Fmat,incx,bufmatb,incx);		cblas_dscal(m2,alpha,bufmatb,incx);		// update X		lambda=2.0/((double)(k)+3.0);		for (j=0;j<m;j++){			for (i=0;i<m;i++){				Xmat[j*m+i]=lambda*dsignf(bufmatb[j*m+i])*dminif(rho*avec[i]*avec[j],dabsf(bufmatb[j*m+i]))+(1-lambda)*dsignf(bufmata[j*m+i])*dminif(rho*avec[i]*avec[j],dabsf(bufmata[j*m+i]));}}		// check convergence and gap periodically		cputime=((double)clock()-start_time)/CLOCKS_PER_SEC;		if ((k%checkgap==0)||(k%Nperiod==0)||(((double)(clock())/CLOCKS_PER_SEC-last_time)>=900)){			gapk=dmax;			for (j=0;j<m;j++){for (i=0;i<m;i++){gapk+=rho*dabsf(Umat[j*m+i])*avec[i]*avec[j];}} // DEBUG: gap not computed properly			alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Bmat,n,Vmat,n,beta,bufmatb,n);			alpha=0.0;cblas_dscal(n2,alpha,Vmat,incx);for(i=0;i<n;i++){gvec[i]=hvec[i]/buf;Vmat[i*n+i]=gvec[i];}			alpha=1.0;beta=0.0;cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,n,n,n,alpha,bufmatb,m,Vmat,n,beta,bufmata,m);			cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,n,n,n,alpha,bufmata,n,bufmatb,n,beta,Vmat,n);				gapk-=doubdot(Amat,Vmat,n2); 			if (firstiter==1) {dualitygap_alliter[checkgap_count]=gapk;cputime_alliter[checkgap_count]=cputime;checkgap_count++;}						last_time=(double)(clock())/CLOCKS_PER_SEC;			if (gapk<=tol) precision_flag=1;			if (k>=MaxIter) iteration_flag=1;			// report iteration, gap and time left		if (((info>=1)&&(k%Nperiod==0)&&(firstiter==1))||(precision_flag+iteration_flag>0)){				left_h=(int)floor(cputime/3600);left_m=(int)floor(cputime/60-left_h*60);left_s=(int)floor(cputime-left_h*3600-left_m*60);				mexPrintf("Iter: %.3e   Obj: %.4e    Gap: %.4e   CPU Time: %2dh %2dm %2ds\n",(double)(k),dmax,gapk,left_h,left_m,left_s);				mexEvalString("drawnow;");}			if (firstiter==0) {firstiter=1;k--;}}		k++;} // End of main loop...	// set dual variable and output vector 	//cblas_dcopy(m2,Umat,incx,Vmat,incx);	//dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack);	//indmax=idxmax(Dvec,n);dmax=Dvec[indmax]; 	//for (i=0;i<n;i++) {uvec[i]=Vmat[(indmax)*n+i];}	*iter=k; // return total number of iterations	// Free everything	free(Vmat);	free(bufmata);	free(bufmatb);	free(bufmatr);	free(bufmatrb);	free(avec);	free(Dvec);	free(qvec);		free(workvec);	free(gvec);	free(hvec);}