DSDP
|
00001 #include "dsdpdualmat.h" 00002 #include "dsdpdsmat.h" 00003 #include "dsdpxmat.h" 00004 #include "dsdpsys.h" 00005 #include "dsdplanczos.h" 00006 #include "dsdplapack.h" 00007 00013 typedef struct _P_Mat3* Mat3; 00014 00015 static int MatMult3(Mat3 A, SDPConeVec x, SDPConeVec y); 00016 static int ComputeStepROBUST(Mat3 A, SDPConeVec *Q, int m, SDPConeVec W, SDPConeVec R,double*, SDPConeVec QAQTv, double *dwork, double *maxstep, double *mineig); 00017 static int ComputeStepFAST(Mat3 A, SDPConeVec *Q, int m, SDPConeVec W, double *dwork, int*iwork,double *maxstep, double *mineig); 00018 00019 extern int DSDPGetEigsSTEP(double[],int,double[],int,long int[],int, 00020 double[],int,double[],int,int[],int); 00021 00022 int DSDPGetTriDiagonalEigs(int n,double *D,double *E, double*WORK2N,int*); 00023 00024 struct _P_Mat3{ 00025 int type; 00026 DSDPDualMat ss; 00027 DSDPDSMat ds; 00028 SDPConeVec V; 00029 DSDPVMat x; 00030 }; 00031 00032 00033 int DSDPGetTriDiagonalEigs(int N,double D[],double E[], double WORK[],int IIWORK[]){ 00034 ffinteger LDZ=DSDPMax(1,N),INFO,NN=N; 00035 ffinteger M,IL=N-1,IU=N,*ISUPPZ=0; 00036 ffinteger *IWORK=(ffinteger*)IIWORK,LIWORK,LWORK; 00037 double WW[2],VL=-1e10,VU=1e10,*Z=0,ABSTOL=0; 00038 char JOBZ='N', RANGE='I'; 00039 if (N<50){ 00040 dstev(&JOBZ,&NN,D,E,Z,&LDZ,WORK,&INFO); 00041 } else { 00042 00043 if (N<=1) IL=1; 00044 if (N<=1) IU=1; 00045 00046 LWORK=20*N+1; 00047 LIWORK=10*N+1; 00048 00049 dstevr(&JOBZ,&RANGE,&NN,D,E,&VL,&VU,&IL,&IU,&ABSTOL,&M,WW,Z,&LDZ,ISUPPZ,WORK,&LWORK,IWORK,&LIWORK,&INFO); 00050 00051 if (N==1){ 00052 D[0]=WW[0]; 00053 } else if (N>=2){ 00054 D[N-2]=WW[0]; 00055 D[N-1]=WW[1]; 00056 } 00057 00058 } 00059 return INFO; 00060 } 00061 00062 /* static int id1=0, id2=0; */ 00063 #undef __FUNCT__ 00064 #define __FUNCT__ "MatMult3" 00065 static int MatMult3(Mat3 A, SDPConeVec X, SDPConeVec Y){ 00066 00067 int info=0; 00068 double minus_one=-1.0; 00069 00070 DSDPFunctionBegin; 00071 /* DSDPEventLogBegin(id2); */ 00072 if (A->type==2){ 00073 info=DSDPVMatMult(A->x,X,Y);DSDPCHKERR(info); 00074 } else { 00075 info=DSDPDualMatCholeskySolveBackward(A->ss,X,Y); DSDPCHKERR(info); 00076 info=DSDPDSMatMult(A->ds,Y,A->V); DSDPCHKERR(info); 00077 info=DSDPDualMatCholeskySolveForward(A->ss,A->V,Y); DSDPCHKERR(info); 00078 info=SDPConeVecScale(minus_one,Y); DSDPCHKERR(info); 00079 } 00080 /* DSDPEventLogEnd(id2);*/ 00081 DSDPFunctionReturn(0); 00082 } 00083 00084 00085 #undef __FUNCT__ 00086 #define __FUNCT__ "DSDPLanczosInitialize" 00087 00092 int DSDPLanczosInitialize( DSDPLanczosStepLength *LZ ){ 00093 DSDPFunctionBegin; 00094 /* Build Lanczos structures */ 00095 LZ->n=0; 00096 LZ->lanczosm=0; 00097 LZ->maxlanczosm=20; 00098 LZ->type=0; 00099 LZ->dwork4n=0; 00100 LZ->Q=0; 00101 LZ->darray=0; 00102 /* 00103 if (id1==0 && id2==0){ 00104 DSDPEventLogRegister("STEP EIGS",&id1); printf("ID1: %d\n",id1); 00105 DSDPEventLogRegister("STEP MULT",&id2); printf("ID2: %d\n",id2); 00106 } 00107 */ 00108 DSDPFunctionReturn(0); 00109 } 00110 00111 #undef __FUNCT__ 00112 #define __FUNCT__ "DSDPSetMaximumLanczosIterations" 00113 00119 int DSDPSetMaximumLanczosIterations( DSDPLanczosStepLength *LZ, int maxlanczos ){ 00120 DSDPFunctionBegin; 00121 if (maxlanczos>0) LZ->lanczosm=maxlanczos; 00122 DSDPFunctionReturn(0); 00123 } 00124 00125 #undef __FUNCT__ 00126 #define __FUNCT__ "DSDPFastLanczosSetup" 00127 00133 int DSDPFastLanczosSetup( DSDPLanczosStepLength *LZ, SDPConeVec V ){ 00134 int i,n,info; 00135 DSDPFunctionBegin; 00136 /* Build Lanczos structures */ 00137 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info); 00138 LZ->lanczosm=DSDPMin(LZ->maxlanczosm,n); 00139 LZ->type=1; 00140 LZ->n=n; 00141 if (LZ->lanczosm<50){ 00142 DSDPCALLOC2(&LZ->dwork4n,double,4*(LZ->lanczosm)+2,&info); DSDPCHKERR(info); 00143 DSDPCALLOC2(&LZ->iwork10n,int,1,&info); DSDPCHKERR(info); 00144 } else { 00145 DSDPCALLOC2(&LZ->dwork4n,double,23*(LZ->lanczosm)+2,&info); DSDPCHKERR(info); 00146 DSDPCALLOC2(&LZ->iwork10n,int,10*(LZ->lanczosm),&info); DSDPCHKERR(info); 00147 } 00148 DSDPCALLOC2(&LZ->Q,SDPConeVec,2,&info);DSDPCHKERR(info); 00149 for (i=0;i<2;i++){ 00150 info = SDPConeVecDuplicate(V,&LZ->Q[i]);DSDPCHKERR(info); 00151 } 00152 DSDPFunctionReturn(0); 00153 } 00154 00155 #undef __FUNCT__ 00156 #define __FUNCT__ "DSDPRobustLanczosSetup" 00157 00163 int DSDPRobustLanczosSetup( DSDPLanczosStepLength *LZ, SDPConeVec V ){ 00164 int i,n,leig,info; 00165 DSDPFunctionBegin; 00166 /* Build Lanczos structures */ 00167 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info); 00168 leig=DSDPMin(LZ->maxlanczosm,n); 00169 LZ->n=n; 00170 LZ->lanczosm=leig; 00171 LZ->type=2; 00172 00173 DSDPCALLOC2(&LZ->dwork4n,double,3*(LZ->lanczosm)+1,&info); DSDPCHKERR(info); 00174 DSDPCALLOC2(&LZ->darray,double,(leig*leig),&info); DSDPCHKERR(info); 00175 DSDPCALLOC2(&LZ->Q,SDPConeVec,(leig+1),&info);DSDPCHKERR(info); 00176 00177 for (i=0;i<=leig;i++){ 00178 info = SDPConeVecDuplicate(V,&LZ->Q[i]);DSDPCHKERR(info); 00179 } 00180 info = SDPConeVecCreate(leig,&LZ->Tv);DSDPCHKERR(info); 00181 DSDPFunctionReturn(0); 00182 } 00183 00184 #undef __FUNCT__ 00185 #define __FUNCT__ "DSDPLanczosDestroy" 00186 00191 int DSDPLanczosDestroy( DSDPLanczosStepLength *LZ){ 00192 int i,info; 00193 DSDPFunctionBegin; 00194 if (LZ->type==2){ 00195 for (i=0;i<=LZ->lanczosm;i++){ 00196 info = SDPConeVecDestroy(&LZ->Q[i]);DSDPCHKERR(info); 00197 } 00198 info=SDPConeVecDestroy(&LZ->Tv);DSDPCHKERR(info); 00199 DSDPFREE(&LZ->darray,&info);DSDPCHKERR(info); 00200 } else if (LZ->type==1){ 00201 info = SDPConeVecDestroy(&LZ->Q[1]);DSDPCHKERR(info); 00202 info = SDPConeVecDestroy(&LZ->Q[0]);DSDPCHKERR(info); 00203 DSDPFREE(&LZ->iwork10n,&info);DSDPCHKERR(info); 00204 } 00205 DSDPFREE(&LZ->Q,&info);DSDPCHKERR(info); 00206 DSDPFREE(&LZ->dwork4n,&info);DSDPCHKERR(info); 00207 info=DSDPLanczosInitialize(LZ);DSDPCHKERR(info); 00208 DSDPFunctionReturn(0); 00209 } 00210 00211 #undef __FUNCT__ 00212 #define __FUNCT__ "DSDPLanczosMinXEig" 00213 int DSDPLanczosMinXEig( DSDPLanczosStepLength *LZ, DSDPVMat X, SDPConeVec W1, SDPConeVec W2, double *mineig ){ 00214 int info,m; 00215 double smaxstep; 00216 struct _P_Mat3 PP; 00217 Mat3 A=&PP; 00218 00219 DSDPFunctionBegin; 00220 A->type=2; 00221 A->x=X; 00222 A->V=W2; 00223 m=LZ->lanczosm; 00224 00225 if (LZ->type==1){ 00226 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,mineig);DSDPCHKERR(info); 00227 } else if (LZ->type==2){ 00228 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray/*LZ->TT*/,LZ->Tv,LZ->dwork4n,&smaxstep,mineig);DSDPCHKERR(info); 00229 } else { 00230 DSDPSETERR1(1,"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type); 00231 } 00232 DSDPFunctionReturn(0); 00233 } 00234 00235 #undef __FUNCT__ 00236 #define __FUNCT__ "DSDPLanczosStepSize" 00237 00247 int DSDPLanczosStepSize( DSDPLanczosStepLength *LZ, SDPConeVec W1, SDPConeVec W2, DSDPDualMat S, DSDPDSMat DS, double *maxstep ){ 00248 int info,m; 00249 double smaxstep,mineig; 00250 struct _P_Mat3 PP; 00251 Mat3 A=&PP; 00252 00253 DSDPFunctionBegin; 00254 A->ss=S; 00255 A->ds=DS; A->V=W2; 00256 A->type=1; 00257 m=LZ->lanczosm; 00258 00259 if (LZ->type==1){ 00260 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,&mineig);DSDPCHKERR(info); 00261 *maxstep=smaxstep; 00262 } else if (LZ->type==2){ 00263 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray/*LZ->TT*/,LZ->Tv,LZ->dwork4n,&smaxstep,&mineig);DSDPCHKERR(info); 00264 *maxstep=smaxstep; 00265 } else { 00266 DSDPSETERR1(1,"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type); 00267 } 00268 DSDPFunctionReturn(0); 00269 } 00270 00271 00272 00273 #undef __FUNCT__ 00274 #define __FUNCT__ "ComputeStepROBUST" 00275 static int ComputeStepROBUST(Mat3 A, SDPConeVec *Q, int m, SDPConeVec W, SDPConeVec R, double*darray, SDPConeVec QAQTv, double *dwork, double *maxstep , double *mineig){ 00276 00277 int i,j,n,info; 00278 double tt,wnorm, phi; 00279 double lambda1=0,lambda2=0,delta=0; 00280 double res1,res2,beta; 00281 double one=1.0; 00282 double *eigvec; 00283 int N, LDA, LWORK; 00284 00285 DSDPFunctionBegin; 00286 00287 memset((void*)darray,0,m*m*sizeof(double)); 00288 if (A->type==1){ 00289 for (i=0; i<m; i++){ darray[i*m+i]=-1.0;} 00290 } else { 00291 for (i=0; i<m; i++){ darray[i*m+i]=1.0;} 00292 } 00293 00294 info = SDPConeVecSet(one,Q[0]);DSDPCHKERR(info); 00295 info = SDPConeVecNormalize(Q[0]);DSDPCHKERR(info); 00296 00297 for (i=0; i<m; i++){ 00298 info = MatMult3(A,Q[i],W);DSDPCHKERR(info); 00299 info = SDPConeVecNorm2(W,&phi);DSDPCHKERR(info); 00300 if (phi!=phi){ *maxstep = 0.0; return 0;} 00301 if (i>0){ 00302 tt=-darray[i*m+i-1]; 00303 info = SDPConeVecAXPY(tt,Q[i-1],W);DSDPCHKERR(info); 00304 } 00305 info = SDPConeVecDot(W,Q[i],&tt);DSDPCHKERR(info); 00306 darray[i*m+i]=tt; 00307 tt*=-1.0; 00308 info = SDPConeVecAXPY(tt,Q[i],W);DSDPCHKERR(info); 00309 info = SDPConeVecNorm2(W,&wnorm);DSDPCHKERR(info); 00310 if (wnorm <= 0.8 * phi){ 00311 for (j=0;j<=i;j++){ 00312 info = SDPConeVecDot(W,Q[j],&tt);DSDPCHKERR(info); 00313 if (tt==tt){tt*=-1.0;} else {tt=0;} 00314 info = SDPConeVecAXPY(tt,Q[j],W);DSDPCHKERR(info); 00315 darray[j*m+i]-=tt; 00316 if (i!=j){ darray[i*m+j]-=tt; } 00317 } 00318 } 00319 00320 info = SDPConeVecNorm2(W,&wnorm);DSDPCHKERR(info); 00321 if (i<m-1){ 00322 darray[i*m+i+1]=wnorm; 00323 darray[i*m+m+i]=wnorm; 00324 } 00325 if (fabs(wnorm)<=1.0e-14) break; 00326 info=SDPConeVecCopy(W,Q[i+1]);DSDPCHKERR(info); 00327 info=SDPConeVecNormalize(Q[i+1]); DSDPCHKERR(info); 00328 00329 } 00330 /* 00331 DSDPLogInfo("Step Length: Lanczos Iterates: %2.0f, ",i); 00332 DSDPLogInfo("VNorm: %3.1e, ",wnorm); 00333 */ 00334 /* 00335 printf(" --- TRI DIAGONAL MATRIX ---- \n"); 00336 */ 00337 00338 00339 LWORK=DSDPMax(3*m,1); LDA=DSDPMax(1,m); N=m; 00340 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info); 00341 info=DSDPGetEigsSTEP(darray,m,0,0,0,0,eigvec,m,dwork,LWORK,0,0);DSDPCHKERR(info); 00342 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info); 00343 00344 if (N==0){ 00345 lambda1=-0.0; 00346 delta=1.0e-20; 00347 *mineig=0; 00348 } else if (N==1){ 00349 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info); 00350 lambda1=-eigvec[0]; 00351 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info); 00352 delta=1.0e-20; 00353 *mineig=lambda1; 00354 } else if (N>1){ 00355 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info); 00356 *mineig=eigvec[0]; 00357 lambda1=-eigvec[N-1]; 00358 lambda2=-eigvec[N-2]; 00359 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info); 00360 00361 info = SDPConeVecZero(W);DSDPCHKERR(info); 00362 for (i=0;i<m;i++){ 00363 tt=darray[(N-1)*m+i]; 00364 info=SDPConeVecAXPY(tt,Q[i],W);DSDPCHKERR(info); 00365 } 00366 info = MatMult3(A,W,R);DSDPCHKERR(info); 00367 info = SDPConeVecAXPY(lambda1,W,R);DSDPCHKERR(info); 00368 info = SDPConeVecNorm2(R,&res1);DSDPCHKERR(info); 00369 00370 info = SDPConeVecZero(W);DSDPCHKERR(info); 00371 for (i=0;i<m;i++){ 00372 tt=darray[(N-2)*m+i]; 00373 info=SDPConeVecAXPY(tt,Q[i],W);DSDPCHKERR(info); 00374 } 00375 info = MatMult3(A,W,R);DSDPCHKERR(info); 00376 info = SDPConeVecAXPY(lambda2,W,R);DSDPCHKERR(info); 00377 info = SDPConeVecNorm2(R,&res2);DSDPCHKERR(info); 00378 00379 tt = -lambda1 + lambda2 - res2; 00380 if (tt>0) beta=tt; 00381 else beta=1.0e-20; 00382 delta = DSDPMin(res1,sqrt(res1)/beta); 00383 00384 } 00385 00386 00387 if (delta-lambda1>0) 00388 *maxstep = 1.0/(delta-lambda1); 00389 else 00390 *maxstep = 1.0e+30; 00391 00392 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info); 00393 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); 00394 00395 DSDPFunctionReturn(0); 00396 } 00397 00398 #undef __FUNCT__ 00399 #define __FUNCT__ "ComputeStepFAST" 00400 static int ComputeStepFAST(Mat3 A, SDPConeVec *Q, int m, SDPConeVec W, double *dwork, int *iwork,double *maxstep ,double *mineig){ 00401 00402 int i,j,n,info; 00403 double tt,wnorm, phi; 00404 double lambda1=0,lambda2=0,delta=0; 00405 double res1,res2,beta; 00406 double one=1.0; 00407 int N=m; 00408 double *diag,*subdiag,*ddwork; 00409 00410 DSDPFunctionBegin; 00411 diag=dwork; 00412 subdiag=dwork+m; 00413 ddwork=dwork+2*m; 00414 00415 if (A->type==1){ 00416 for (i=0; i<m; i++){ diag[i]=-1; subdiag[i]=0;} 00417 } else { 00418 for (i=0; i<m; i++){ diag[i]=1.0; subdiag[i]=0;} 00419 } 00420 info = SDPConeVecSet(one,Q[0]);DSDPCHKERR(info); 00421 info = SDPConeVecNormalize(Q[0]);DSDPCHKERR(info); 00422 00423 for (i=0; i<m; i++){ 00424 info = MatMult3(A,Q[0],W);DSDPCHKERR(info); 00425 info = SDPConeVecNorm2(W,&phi);DSDPCHKERR(info); 00426 if (phi!=phi){ *maxstep = 0.0; return 0;} 00427 if (i>0){ 00428 tt=-subdiag[i-1]; 00429 info = SDPConeVecAXPY(tt,Q[1],W);DSDPCHKERR(info); 00430 } 00431 info = SDPConeVecDot(W,Q[0],&tt);DSDPCHKERR(info); 00432 diag[i]=tt; 00433 tt*=-1.0; 00434 info = SDPConeVecAXPY(tt,Q[0],W);DSDPCHKERR(info); 00435 info = SDPConeVecNorm2(W,&wnorm);DSDPCHKERR(info); 00436 if (wnorm <= 1.0 * phi){ 00437 for (j=0;j<=i;j++){ 00438 if (j==i-1){ 00439 info = SDPConeVecDot(W,Q[1],&tt);DSDPCHKERR(info); 00440 if (tt==tt){tt*=-1.0;} else {tt=0;} 00441 info = SDPConeVecAXPY(tt,Q[1],W);DSDPCHKERR(info); 00442 subdiag[i-1]-=tt; 00443 } else if (j==i){ 00444 info = SDPConeVecDot(W,Q[0],&tt);DSDPCHKERR(info); 00445 if (tt==tt){tt*=-1.0;} else {tt=0;} 00446 info = SDPConeVecAXPY(tt,Q[0],W);DSDPCHKERR(info); 00447 diag[i]-=tt; 00448 } 00449 00450 } 00451 } 00452 00453 info = SDPConeVecNorm2(W,&wnorm);DSDPCHKERR(info); 00454 /* printf("PHI: %4.4e, VNORM: %4.2e Diag: %4.2e\n",phi,wnorm,diag[i]); */ 00455 if (i<m-1){ 00456 subdiag[i]=wnorm; 00457 } 00458 if (fabs(wnorm)<=1.0e-10){i++;break;} 00459 info=SDPConeVecCopy(Q[0],Q[1]);DSDPCHKERR(info); 00460 info=SDPConeVecCopy(W,Q[0]);DSDPCHKERR(info); 00461 info=SDPConeVecNormalize(Q[0]); DSDPCHKERR(info); 00462 00463 } 00464 00465 /* DSDPEventLogBegin(id1); */ 00466 info=DSDPGetTriDiagonalEigs(m,diag,subdiag,ddwork,iwork); DSDPCHKERR(info); 00467 /* DSDPEventLogEnd(id1); */ 00468 if (N==0){ 00469 lambda1=-0.0; 00470 delta=1.0e-20; 00471 *mineig=0; 00472 } else if (N==1){ 00473 lambda1=-diag[0]; 00474 delta=1.0e-20; 00475 *mineig=diag[0]; 00476 } else if (N>1){ 00477 lambda1=-diag[N-1]; 00478 lambda2=-diag[N-2]; 00479 00480 res1=1.0e-8; 00481 res2=1.0e-8; 00482 00483 tt = -lambda1 + lambda2 - res2; 00484 if (tt>0) beta=tt; 00485 else beta=1.0e-20; 00486 delta = DSDPMin(res1,sqrt(res1)/beta); 00487 00488 *mineig=diag[0]; 00489 } 00490 00491 00492 if (delta-lambda1>0) 00493 *maxstep = 1.0/(delta-lambda1); 00494 else 00495 *maxstep = 1.0e+30; 00496 00497 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info); 00498 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", 00499 i,m,n,wnorm,lambda1,lambda2,delta,*maxstep); 00500 00501 DSDPFunctionReturn(0); 00502 }