/* 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: A. d'Aspremont, Laurent El Ghaoui, Ronny Luss November 2007.http://www.carva.org/alexandre.daspremont*/#include "sparsesvd.h"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){	// Hard parameters	int Nperiod=imaxf(1,info);	int work_size=3*n+n*n;	// Working variables	double d1,sig1,d2,sig2,norma12,mu,Ntheo,L;	double alpha,gapk,dmax=0.0,fmu,lambda;	int n2=n*n,incx=1,precision_flag=0,iteration_flag=0;	int lwork=work_size,inflapack,indmax,k=0,i,j;	double cputime,last_time=(double)clock();double start_time=(double)clock();int left_h=0,left_m=0,left_s=0;	char jobz[1]="V",uplo[1]="U";	double *Vmat=(double *) calloc(n*n,sizeof(double));	double *bufmata=(double *) calloc(n*n,sizeof(double));	double *bufmatb=(double *) calloc(n*n,sizeof(double));	double *Dvec=(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 work_size3=8*n;	double *workvec2=(double *) calloc(work_size3,sizeof(double));	int *iwork=(int *) calloc(5*n,sizeof(int));	double *numeigs_matlab=(double *) calloc(1,sizeof(double));	double *evector_store=(double *) calloc(n*n,sizeof(double));	double *evector_temp=(double *) calloc(n*n,sizeof(double)); // TODO: one of these matrices can go.	double *eig=(double *) calloc(n,sizeof(double));	double eigcut,tolweight=.75,tol=.01,numeigstemp=numeigs;	int *count=(int *) calloc(n,sizeof(int)); // Records the distribution of eig computations	int checkgap_count=0,firstiter=0; // added for test variables	int arcount;	// Start...	if (info>=1){mexPrintf("DSPCA starting... Sparse eig. maximization using approximate gradient.\n");mexEvalString("drawnow;");}	// Test malloc results	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL)|| (Dvec==NULL) || (workvec==NULL) || (gvec==NULL) || (hvec==NULL) ||(evector_temp==NULL)||(evector_store==NULL)||(iwork==NULL)||(workvec2==NULL)||(numeigs_matlab==NULL)||(eig==NULL)){		mexPrintf("DSPCA: memory allocation failed ... \n");mexEvalString("drawnow;");return;}	eigcut=(1-tolweight)*(tol/10)/(rho*n); // scale delta (tol/10) to get eig threshold	tol=tolweight*tol; // scale of .5 for partial eig precision	// First, compute some local params	norma12=1.0;d1=rho*rho*n*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;mu=tol/(2.0*d2); // TODO: can we get a less conservative d1?	Ntheo=(4.0*norma12*sqrt(d1*d2/(sig1*sig2)))/tol;Ntheo=ceil(Ntheo);	L=(d2*norma12*norma12)/(2.0*sig2*tol);	alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);	cputime=start_time;	while ((precision_flag+iteration_flag)==0){				// eigenvalue decomposition of A+X 		cblas_dcopy(n2,Xmat,incx,Vmat,incx);		alpha=1.0; 		cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);		cblas_dcopy(n2,Vmat,incx,bufmata,incx);				// do partial eigenvalue approximation to exp(A+X)		*numeigs_matlab=1.0*numeigs;*count=0;arcount=0;		fmu=partial_eig(n,k,mu,eigcut,bufmata,bufmatb,numeigs_matlab,evector_temp,evector_store,eig,Dvec,gvec,hvec,Vmat,Umat,workvec,count,addeigs,perceigs,check_for_more_eigs,&arcount); 		numeigs=(int)(*numeigs_matlab);dmax=bufmata[0];		// update gradient's weighted average 		alpha=((double)(k)+1)/2.0;		cblas_daxpy(n2,alpha,Umat,incx,Fmat,incx);		// find a projection of X-Gmu/L on feasible set 		cblas_dcopy(n2,Xmat,incx,bufmata,incx);		alpha=-1.0/L;		cblas_daxpy(n2,alpha,Umat,incx,bufmata,incx);		// project again		alpha=-(sig1/L);		cblas_dcopy(n2,Fmat,incx,bufmatb,incx);		cblas_dscal(n2,alpha,bufmatb,incx);		// update X		lambda=2.0/((double)(k)+3.0);		for (j=0;j<n;j++){			for (i=0;i<n;i++){				Xmat[j*n+i]=lambda*dsignf(bufmatb[j*n+i])*dminif(rho,dabsf(bufmatb[j*n+i]))+(1-lambda)*dsignf(bufmata[j*n+i])*dminif(rho,dabsf(bufmata[j*n+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-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);			if (firstiter==1) {dualitygap_alliter[checkgap_count]=gapk;cputime_alliter[checkgap_count]=cputime;perceigs_alliter[checkgap_count]=100.0*numeigs/n;checkgap_count++;}			if (firstiter==0){// If first iteration, reset precision targets				tol=gapk*gapchange;norma12=1.0;d1=rho*rho*n*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;mu=tol/(2.0*d2);				L=(d2*norma12*norma12)/(2.0*sig2*tol);eigcut=(1-tolweight)*(tol/10)/(rho*n);				alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);cblas_dscal(n2,alpha,Fmat,incx);numeigs=numeigstemp;}					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	 %% Eigs Used: %.2f\n",(double)(k),dmax,gapk,left_h,left_m,left_s,100.0*numeigs/n);				mexEvalString("drawnow;");}			if (firstiter==0) {firstiter=1;k--;}}		k++;}	// set dual variable and output vector	// eigenvalue decomposition of A+X 	alpha=0.0;cblas_dscal(n2,alpha,Vmat,incx);	cblas_dcopy(n2,Umat,incx,Vmat,incx);	dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack); // TODO: switch to ARPACK here...	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(Dvec);	free(workvec);	free(workvec2);	free(gvec);	free(hvec);	free(iwork);	free(numeigs_matlab);	free(evector_store);	free(evector_temp);	free(eig);	free(count);}