Main Page | Modules | Alphabetical List | Data Structures | Directories | File List | Globals | Related Pages

dsdpstep.c

Go to the documentation of this file.
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 }

Generated on Fri Oct 21 14:28:34 2005 for DSDP by  doxygen 1.4.2