/* 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(double *Amat, int n, double rho, double gapchange, 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 d1,sig1,d2,sig2,norma12,mu,Ntheo,L;	double alpha,beta,buf,gapk;	double dmax=0.0,fmu,lambda,tol=.01;	int n2=n*n,incx=1,precision_flag=0,iteration_flag=0;	int lwork,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],uplo[1],transa[1],transb[1];	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 checkgap_count=0,firstiter=0; // added for test variables		// Start...	if (info>=1){		mexPrintf("DSPCA starting... Sparse 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;}	// First, compute some local params	d1=rho*rho*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;norma12=1.0;mu=tol/(2.0*d2);	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);		*jobz='V';*uplo='U';lwork=work_size;		dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack); // call LAPACK (most CPU time is here) 		// compute fmu(X) = mu*log(trace((exp(A+X)/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=0.0;cblas_dscal(n2,alpha,bufmatb,incx);		for (i=0;i<n;i++) {gvec[i]=hvec[i]/buf;bufmatb[i*n+i]=gvec[i];}		alpha=1.0;beta=0.0;*transa='N';*transb='T';		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Vmat,n,bufmatb,n,beta,bufmata,n);		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,n,n,n,alpha,bufmata,n,Vmat,n,beta,Umat,n);		// 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/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);		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;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);				alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);cblas_dscal(n2,alpha,Fmat,incx);}					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	// eigenvalue decomposition of A+X 	alpha=0.0;cblas_dscal(n2,alpha,Vmat,incx);	cblas_dcopy(n2,Umat,incx,Vmat,incx);	*jobz='V';*uplo='U';lwork=work_size;	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(Dvec);	free(workvec);	free(gvec);	free(hvec);}