13 typedef struct _P_Mat3* Mat3;
17 static int ComputeStepFAST(Mat3 A,
SDPConeVec *Q,
int m,
SDPConeVec W,
double *dwork,
int*iwork,
double *maxstep,
double *mineig);
19 extern int DSDPGetEigsSTEP(
double[],
int,
double[],
int,
long int[],
int,
20 double[],
int,
double[],
int,
int[],
int);
22 int DSDPGetTriDiagonalEigs(
int n,
double *D,
double *E,
double*WORK2N,
int*);
33 int DSDPGetTriDiagonalEigs(
int N,
double D[],
double E[],
double WORK[],
int IIWORK[]){
34 ffinteger LDZ=DSDPMax(1,N),INFO,NN=N;
35 ffinteger M,IL=N-1,IU=N,*ISUPPZ=0;
36 ffinteger *IWORK=(ffinteger*)IIWORK,LIWORK,LWORK;
37 double WW[2],VL=-1e10,VU=1e10,*Z=0,ABSTOL=0;
38 char JOBZ=
'N', RANGE=
'I';
40 dstev(&JOBZ,&NN,D,E,Z,&LDZ,WORK,&INFO);
49 dstevr(&JOBZ,&RANGE,&NN,D,E,&VL,&VU,&IL,&IU,&ABSTOL,&M,WW,Z,&LDZ,ISUPPZ,WORK,&LWORK,IWORK,&LIWORK,&INFO);
64 #define __FUNCT__ "MatMult3"
68 double minus_one=-1.0;
81 DSDPFunctionReturn(0);
86 #define __FUNCT__ "DSDPLanczosInitialize"
108 DSDPFunctionReturn(0);
112 #define __FUNCT__ "DSDPSetMaximumLanczosIterations"
121 if (maxlanczos>0) LZ->lanczosm=maxlanczos;
122 DSDPFunctionReturn(0);
126 #define __FUNCT__ "DSDPFastLanczosSetup"
137 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info);
138 LZ->lanczosm=DSDPMin(LZ->maxlanczosm,n);
141 if (LZ->lanczosm<50){
142 DSDPCALLOC2(&LZ->dwork4n,
double,4*(LZ->lanczosm)+2,&info); DSDPCHKERR(info);
143 DSDPCALLOC2(&LZ->iwork10n,
int,1,&info); DSDPCHKERR(info);
145 DSDPCALLOC2(&LZ->dwork4n,
double,23*(LZ->lanczosm)+2,&info); DSDPCHKERR(info);
146 DSDPCALLOC2(&LZ->iwork10n,
int,10*(LZ->lanczosm),&info); DSDPCHKERR(info);
148 DSDPCALLOC2(&LZ->Q,
SDPConeVec,2,&info);DSDPCHKERR(info);
152 DSDPFunctionReturn(0);
156 #define __FUNCT__ "DSDPRobustLanczosSetup"
167 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info);
168 leig=DSDPMin(LZ->maxlanczosm,n);
173 DSDPCALLOC2(&LZ->dwork4n,
double,3*(LZ->lanczosm)+1,&info); DSDPCHKERR(info);
174 DSDPCALLOC2(&LZ->darray,
double,(leig*leig),&info); DSDPCHKERR(info);
175 DSDPCALLOC2(&LZ->Q,
SDPConeVec,(leig+1),&info);DSDPCHKERR(info);
177 for (i=0;i<=leig;i++){
180 info = SDPConeVecCreate(leig,&LZ->Tv);DSDPCHKERR(info);
181 DSDPFunctionReturn(0);
185 #define __FUNCT__ "DSDPLanczosDestroy"
195 for (i=0;i<=LZ->lanczosm;i++){
196 info = SDPConeVecDestroy(&LZ->Q[i]);DSDPCHKERR(info);
198 info=SDPConeVecDestroy(&LZ->Tv);DSDPCHKERR(info);
199 DSDPFREE(&LZ->darray,&info);DSDPCHKERR(info);
200 }
else if (LZ->type==1){
201 info = SDPConeVecDestroy(&LZ->Q[1]);DSDPCHKERR(info);
202 info = SDPConeVecDestroy(&LZ->Q[0]);DSDPCHKERR(info);
203 DSDPFREE(&LZ->iwork10n,&info);DSDPCHKERR(info);
205 DSDPFREE(&LZ->Q,&info);DSDPCHKERR(info);
206 DSDPFREE(&LZ->dwork4n,&info);DSDPCHKERR(info);
208 DSDPFunctionReturn(0);
212 #define __FUNCT__ "DSDPLanczosMinXEig"
226 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,mineig);DSDPCHKERR(info);
227 }
else if (LZ->type==2){
228 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray,LZ->Tv,LZ->dwork4n,&smaxstep,mineig);DSDPCHKERR(info);
230 DSDPSETERR1(1,
"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type);
232 DSDPFunctionReturn(0);
236 #define __FUNCT__ "DSDPLanczosStepSize"
249 double smaxstep,mineig;
260 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,&mineig);DSDPCHKERR(info);
262 }
else if (LZ->type==2){
263 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray,LZ->Tv,LZ->dwork4n,&smaxstep,&mineig);DSDPCHKERR(info);
266 DSDPSETERR1(1,
"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type);
268 DSDPFunctionReturn(0);
274 #define __FUNCT__ "ComputeStepROBUST"
278 double tt,wnorm, phi;
279 double lambda1=0,lambda2=0,delta=0;
280 double res1,res2,beta;
287 memset((
void*)darray,0,m*m*
sizeof(
double));
289 for (i=0; i<m; i++){ darray[i*m+i]=-1.0;}
291 for (i=0; i<m; i++){ darray[i*m+i]=1.0;}
298 info = MatMult3(A,Q[i],W);DSDPCHKERR(info);
300 if (phi!=phi){ *maxstep = 0.0;
return 0;}
310 if (wnorm <= 0.8 * phi){
313 if (tt==tt){tt*=-1.0;}
else {tt=0;}
316 if (i!=j){ darray[i*m+j]-=tt; }
322 darray[i*m+i+1]=wnorm;
323 darray[i*m+m+i]=wnorm;
325 if (fabs(wnorm)<=1.0e-14)
break;
339 LWORK=DSDPMax(3*m,1); LDA=DSDPMax(1,m); N=m;
340 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
341 info=DSDPGetEigsSTEP(darray,m,0,0,0,0,eigvec,m,dwork,LWORK,0,0);DSDPCHKERR(info);
342 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
349 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
351 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
355 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
357 lambda1=-eigvec[N-1];
358 lambda2=-eigvec[N-2];
359 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
363 tt=darray[(N-1)*m+i];
366 info = MatMult3(A,W,R);DSDPCHKERR(info);
372 tt=darray[(N-2)*m+i];
375 info = MatMult3(A,W,R);DSDPCHKERR(info);
379 tt = -lambda1 + lambda2 - res2;
382 delta = DSDPMin(res1,sqrt(res1)/beta);
388 *maxstep = 1.0/(delta-lambda1);
392 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info);
393 DSDPLogInfo(0,19,
"Robust Lanczos StepLength: Iterates %d, Max: %d, BlockSize: %d, Lambda1: %4.2e, Res1: %4.2e, Lambda2: %4.2e, Res2: %4.2e, Delta: %4.2e, MaxStep: %4.2e\n",i,m,n,lambda1,res1*res1,lambda2,res2*res2,delta,*maxstep);
395 DSDPFunctionReturn(0);
399 #define __FUNCT__ "ComputeStepFAST"
400 static int ComputeStepFAST(Mat3 A,
SDPConeVec *Q,
int m,
SDPConeVec W,
double *dwork,
int *iwork,
double *maxstep ,
double *mineig){
403 double tt,wnorm, phi;
404 double lambda1=0,lambda2=0,delta=0;
405 double res1,res2,beta;
408 double *diag,*subdiag,*ddwork;
416 for (i=0; i<m; i++){ diag[i]=-1; subdiag[i]=0;}
418 for (i=0; i<m; i++){ diag[i]=1.0; subdiag[i]=0;}
424 info = MatMult3(A,Q[0],W);DSDPCHKERR(info);
426 if (phi!=phi){ *maxstep = 0.0;
return 0;}
436 if (wnorm <= 1.0 * phi){
440 if (tt==tt){tt*=-1.0;}
else {tt=0;}
445 if (tt==tt){tt*=-1.0;}
else {tt=0;}
458 if (fabs(wnorm)<=1.0e-10){i++;
break;}
466 info=DSDPGetTriDiagonalEigs(m,diag,subdiag,ddwork,iwork); DSDPCHKERR(info);
483 tt = -lambda1 + lambda2 - res2;
486 delta = DSDPMin(res1,sqrt(res1)/beta);
493 *maxstep = 1.0/(delta-lambda1);
497 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info);
498 DSDPLogInfo(0,19,
"Step Length: Fast Lanczos Iterates: %2d, Max: %d, Block Size: %d, VNorm: %3.1e, Lambda1: %4.4e, Lambda2: %4.4e, Delta: %4.2e, Maxstep: %4.2e\n",
499 i,m,n,wnorm,lambda1,lambda2,delta,*maxstep);
501 DSDPFunctionReturn(0);