dsdpcg.c

Go to the documentation of this file.
00001 #include "dsdpcg.h"
00002 #include "dsdpschurmat.h"
00003 #include "dsdpvec.h"
00004 #include "dsdpsys.h"
00005 #include "dsdp.h"
00012 typedef enum { DSDPNoMatrix=1, DSDPUnfactoredMatrix=2, DSDPFactoredMatrix=3} DSDPCGType;
00013 
00014 typedef struct{
00015   DSDPCGType type;
00016   DSDPSchurMat M;
00017   DSDPVec Diag;
00018   DSDP dsdp;
00019 } DSDPCGMat;
00020 
00021 #undef __FUNCT__  
00022 #define __FUNCT__ "DSDPCGMatMult"
00023 int DSDPCGMatMult(DSDPCGMat M, DSDPVec X, DSDPVec Y){
00024   int info;
00025   DSDPFunctionBegin;
00026   info=DSDPVecZero(Y); DSDPCHKERR(info);
00027   if (M.type==DSDPUnfactoredMatrix){
00028     info=DSDPSchurMatMultiply(M.M,X,Y); DSDPCHKERR(info);
00029   } else if (M.type==DSDPFactoredMatrix){
00030     info=DSDPSchurMatMultR(M.M,X,Y); DSDPCHKERR(info);
00031     info=DSDPVecAXPY(-0*M.dsdp->Mshift,X,Y); DSDPCHKERR(info);
00032   } else if (M.type==DSDPNoMatrix){
00033     info=DSDPHessianMultiplyAdd(M.dsdp,X,Y);DSDPCHKERR(info);
00034   }
00035   DSDPFunctionReturn(0);
00036 }
00037 
00038 #undef __FUNCT__  
00039 #define __FUNCT__ "DSDPCGMatPre"
00040 int DSDPCGMatPre(DSDPCGMat M, DSDPVec X, DSDPVec Y){
00041   int info;
00042   DSDPFunctionBegin;
00043   info=DSDPVecZero(Y); DSDPCHKERR(info);
00044   if (M.type==DSDPUnfactoredMatrix){
00045     info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
00046     info=DSDPVecPointwiseMult(Y,M.Diag,Y);DSDPCHKERR(info);
00047   } else if (M.type==DSDPFactoredMatrix){
00048     info=DSDPSchurMatSolve(M.M,X,Y);DSDPCHKERR(info);
00049   } else if (M.type==DSDPNoMatrix){
00050     info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
00051   }
00052   DSDPFunctionReturn(0);
00053 }
00054 
00055 #undef __FUNCT__  
00056 #define __FUNCT__ "DSDPCGMatPreLeft"
00057 int DSDPCGMatPreLeft(DSDPCGMat M, DSDPVec X, DSDPVec Y){
00058   int info;
00059   DSDPFunctionBegin;
00060   info=DSDPVecZero(Y); DSDPCHKERR(info);
00061   if (M.type==DSDPUnfactoredMatrix){
00062     info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
00063   } else if (M.type==DSDPFactoredMatrix){
00064     info=DSDPSchurMatSolve(M.M,X,Y);DSDPCHKERR(info);
00065   } else if (M.type==DSDPNoMatrix){
00066     info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
00067   }
00068   DSDPFunctionReturn(0);
00069 }
00070 
00071 #undef __FUNCT__  
00072 #define __FUNCT__ "DSDPCGMatPreRight"
00073 int DSDPCGMatPreRight(DSDPCGMat M, DSDPVec X, DSDPVec Y){
00074   int info;
00075   DSDPFunctionBegin;
00076   info=DSDPVecZero(Y); DSDPCHKERR(info);
00077   if (M.type==DSDPNoMatrix){
00078     info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
00079   } else if (M.type==DSDPFactoredMatrix){
00080     info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
00081   } else if (M.type==DSDPUnfactoredMatrix){
00082     info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
00083   }
00084   DSDPFunctionReturn(0);
00085 }
00086 
00087 
00088 #undef __FUNCT__  
00089 #define __FUNCT__ "DSDPConjugateResidual"
00090 int DSDPConjugateResidual(DSDPCGMat B, DSDPVec X, DSDPVec D, DSDPVec R, DSDPVec BR, DSDPVec P, DSDPVec BP, DSDPVec TT3, int maxits, int *iter){
00091 
00092   int i,n,info;
00093   double zero=0.0,minus_one=-1.0;
00094   double alpha,beta,bpbp,rBr,rBrOld;
00095   double r0,rerr=1.0e20;
00096 
00097   DSDPFunctionBegin;
00098   info=DSDPVecNorm2(X,&rBr); DSDPCHKERR(info);
00099   if (rBr>0){
00100     info=DSDPVecCopy(X,P); DSDPCHKERR(info);
00101     info=DSDPCGMatPreRight(B,P,X);DSDPCHKERR(info);
00102     info=DSDPCGMatMult(B,X,R); DSDPCHKERR(info);
00103   } else {
00104     info=DSDPVecSet(zero,R); DSDPCHKERR(info);
00105   }
00106   info=DSDPVecAYPX(minus_one,D,R); DSDPCHKERR(info);
00107 
00108   info=DSDPCGMatPreLeft(B,D,R);DSDPCHKERR(info);
00109   info=DSDPVecCopy(R,P); DSDPCHKERR(info);
00110 
00111   info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
00112   info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
00113   info=DSDPCGMatPreRight(B,TT3,BR);DSDPCHKERR(info);
00114 
00115   info=DSDPVecCopy(BR,BP); DSDPCHKERR(info);
00116   info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
00117   info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
00118   r0=rBr;
00119 
00120   for (i=0;i<maxits;i++){
00121 
00122     if (rerr/n < 1.0e-30 || rBr/n <= 1.0e-30  || rBr< 1.0e-12 * r0 ) break;
00123 
00124     info=DSDPVecDot(BP,BP,&bpbp); DSDPCHKERR(info);
00125     alpha = rBr / bpbp;
00126     info= DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
00127     alpha = -alpha;
00128     info=DSDPVecAXPY(alpha,BP,R); DSDPCHKERR(info);
00129 
00130     info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
00131     info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
00132     info=DSDPCGMatPreLeft(B,TT3,BR);DSDPCHKERR(info);
00133     
00134     rBrOld=rBr;
00135     info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
00136     info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
00137 
00138     DSDPLogInfo(0,11,"CG: rerr: %4.4e, rBr: %4.4e \n",rerr,rBr);  
00139 
00140     beta = rBr/rBrOld;
00141     info=DSDPVecAYPX(beta,R,P); DSDPCHKERR(info);
00142     info=DSDPVecAYPX(beta,BR,BP); DSDPCHKERR(info);
00143     
00144   }
00145   info=DSDPVecCopy(X,BR);DSDPCHKERR(info);
00146   info=DSDPCGMatPreRight(B,BR,X);DSDPCHKERR(info);
00147 
00148   DSDPLogInfo(0,2,"Conjugate Residual, Initial rMr, %4.2e, Final rMr: %4.2e, Iterates: %d\n",r0,rBr,i);
00149 
00150   *iter=i;
00151 
00152   DSDPFunctionReturn(0);
00153 }
00154 
00155 #undef __FUNCT__  
00156 #define __FUNCT__ "DSDPConjugateGradient"
00157 int DSDPConjugateGradient(DSDPCGMat B, DSDPVec X, DSDPVec D, DSDPVec R, DSDPVec BR, DSDPVec P, DSDPVec BP, DSDPVec Z, double cgtol, int maxits, int *iter){
00158 
00159   int i,n,info;
00160   double alpha,beta=0,bpbp;
00161   double r0,rerr=1.0e20;
00162   double rz,rzold,rz0;
00163 
00164   DSDPFunctionBegin;
00165   if (maxits<=0){
00166     *iter=0;
00167     DSDPFunctionReturn(0);
00168   } 
00169   info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
00170   info=DSDPVecNorm2(X,&rerr); DSDPCHKERR(info);
00171   info=DSDPVecZero(R); DSDPCHKERR(info);
00172   if (rerr>0){
00173     info=DSDPCGMatMult(B,X,R);DSDPCHKERR(info);
00174   }
00175   alpha=-1.0;
00176   info=DSDPVecAYPX(alpha,D,R); DSDPCHKERR(info);
00177   info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);    
00178   if (rerr*sqrt(1.0*n)<1e-11 +0*cgtol*cgtol){
00179     *iter=1;
00180     DSDPFunctionReturn(0);
00181   } 
00182 
00183   if (rerr>0){
00184     info=DSDPVecCopy(R,P); DSDPCHKERR(info);
00185     info=DSDPVecSetR(P,0);DSDPCHKERR(info);
00186     info=DSDPVecNorm2(P,&rerr); DSDPCHKERR(info);    
00187     info=DSDPCGMatPre(B,R,Z);DSDPCHKERR(info);
00188   }
00189     
00190   info=DSDPVecCopy(Z,P); DSDPCHKERR(info);
00191   info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
00192   r0=rerr;rz0=rz;
00193   
00194   for (i=0;i<maxits;i++){
00195     
00196     info=DSDPCGMatMult(B,P,BP);DSDPCHKERR(info);
00197     info=DSDPVecDot(BP,P,&bpbp); DSDPCHKERR(info);
00198     if (bpbp==0) break;
00199     alpha = rz/bpbp;
00200     info=DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
00201     info=DSDPVecAXPY(-alpha,BP,R); DSDPCHKERR(info);
00202     info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
00203 
00204     DSDPLogInfo(0,15,"CG: iter: %d, rerr: %4.4e, alpha: %4.4e, beta=%4.4e, rz: %4.4e \n",i,rerr,alpha,beta,rz);  
00205     if (i>1){
00206       if (rerr*sqrt(1.0*n) < cgtol){ break;}
00207       if (rz*n < cgtol*cgtol){ break;}
00208       if (rerr/r0 < cgtol){ break;}
00209     }
00210     if (rerr<=0){ break;}
00211     info=DSDPCGMatPre(B,R,Z);DSDPCHKERR(info);
00212     rzold=rz;
00213     info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
00214     beta=rz/rzold;
00215     info= DSDPVecAYPX(beta,Z,P); DSDPCHKERR(info);
00216   }
00217   DSDPLogInfo(0,2,"Conjugate Gradient, Initial ||r||: %4.2e, Final ||r||: %4.2e, Initial ||rz||: %4.4e, ||rz||: %4.4e, Iterates: %d\n",r0,rerr,rz0,rz,i+1);
00218 
00219   *iter=i+1;
00220 
00221   DSDPFunctionReturn(0);
00222 }
00223 
00237 #undef __FUNCT__  
00238 #define __FUNCT__ "DSDPCGSolve"
00239 int DSDPCGSolve(DSDP dsdp, DSDPSchurMat MM, DSDPVec RHS, DSDPVec X,double cgtol, DSDPTruth *success){
00240   int iter=0,n,info,maxit=10;
00241   double dd,ymax;
00242   DSDPCG *sles=dsdp->sles; 
00243   DSDPCGMat    CGM;
00244   DSDPFunctionBegin;
00245 
00246   info=DSDPEventLogBegin(dsdp->cgtime);
00247   info=DSDPVecZero(X);DSDPCHKERR(info);
00248   info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
00249   *success=DSDP_TRUE;
00250   if (0){
00251     maxit=0;
00252   } else if (dsdp->slestype==1){
00253 
00254     CGM.type=DSDPNoMatrix;
00255     CGM.M=MM;
00256     CGM.dsdp=dsdp;
00257     cgtol*=1000;
00258     maxit=5;
00259 
00260   } else if (dsdp->slestype==2){
00261     CGM.type=DSDPUnfactoredMatrix;
00262     CGM.M=MM;
00263     CGM.Diag=sles->Diag;
00264     CGM.dsdp=dsdp;
00265     cgtol*=100;
00266     maxit=(int)sqrt(1.0*n)+10;
00267     if (maxit>20) maxit=20;
00268     info=DSDPVecSet(1.0,sles->Diag);DSDPCHKERR(info);
00269     /*
00270       info=DSDPSchurMatGetDiagonal(MM,sles->Diag);DSDPCHKERR(info);
00271       info=DSDPVecReciprocalSqrt(sles->Diag); DSDPCHKERR(info);
00272     */
00273     
00274   } else if (dsdp->slestype==3){
00275     CGM.type=DSDPFactoredMatrix;
00276     CGM.M=MM;
00277     CGM.dsdp=dsdp;
00278     maxit=0;
00279     info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
00280     if (ymax > 1e5 && dsdp->rgap<1e-1) maxit=3;
00281     if (0 && ymax > 1e5 && dsdp->rgap<1e-2){ 
00282       maxit=6;
00283     } else if (dsdp->rgap<1e-5){
00284       maxit=3;
00285     }
00286 
00287     info=DSDPSchurMatSolve(MM,RHS,X);DSDPCHKERR(info);
00288 
00289   } else if (dsdp->slestype==4){
00290     CGM.type=DSDPFactoredMatrix;
00291     CGM.M=MM;
00292     CGM.dsdp=dsdp;
00293     maxit=3;
00294     info=DSDPSchurMatSolve(MM,RHS,X);DSDPCHKERR(info);
00295   }
00296   if (n<maxit) maxit=n;
00297   
00298   info=DSDPConjugateGradient(CGM,X,RHS,
00299                              sles->R,sles->BR,sles->P,sles->BP,
00300                              sles->TTT,cgtol,maxit,&iter);DSDPCHKERR(info);
00301 
00302   if (iter>=maxit) *success=DSDP_FALSE;
00303   info=DSDPVecDot(RHS,X,&dd);DSDPCHKERR(info);
00304   if (dd<0) *success=DSDP_FALSE;
00305   info=DSDPEventLogEnd(dsdp->cgtime);
00306   DSDPFunctionReturn(0);
00307 }
00308 
00309 
00310 #undef __FUNCT__  
00311 #define __FUNCT__ "DSDPCGInitialize"
00312 int DSDPCGInitialize(DSDPCG **sles){
00313   DSDPCG *ssles;
00314   int info;
00315   DSDPFunctionBegin;
00316   DSDPCALLOC1(&ssles,DSDPCG,&info);DSDPCHKERR(info);
00317   ssles->setup2=0;
00318   *sles=ssles;
00319   DSDPFunctionReturn(0);
00320 }
00321 
00322 #undef __FUNCT__  
00323 #define __FUNCT__ "DSDPCGSetup"
00324 int DSDPCGSetup(DSDPCG *sles,DSDPVec X){
00325   int info;
00326   DSDPFunctionBegin;
00327   info=DSDPVecGetSize(X,&sles->m);DSDPCHKERR(info);
00328   if (sles->setup2==0){
00329     info = DSDPVecDuplicate(X,&sles->R); DSDPCHKERR(info);
00330     info = DSDPVecDuplicate(X,&sles->P); DSDPCHKERR(info);
00331     info = DSDPVecDuplicate(X,&sles->BP); DSDPCHKERR(info);
00332     info = DSDPVecDuplicate(X,&sles->BR); DSDPCHKERR(info);
00333     info = DSDPVecDuplicate(X,&sles->Diag); DSDPCHKERR(info);
00334     info = DSDPVecDuplicate(X,&sles->TTT); DSDPCHKERR(info);
00335   }
00336   sles->setup2=1;
00337   DSDPFunctionReturn(0);
00338 }
00339 
00340 #undef __FUNCT__  
00341 #define __FUNCT__ "DSDPCGDestroy"
00342 int DSDPCGDestroy(DSDPCG **ssles){
00343   int info;
00344   DSDPCG *sles=*ssles;
00345   DSDPFunctionBegin;
00346   if (ssles==0 || sles==0){DSDPFunctionReturn(0);}
00347   if (sles->setup2==1){
00348     info = DSDPVecDestroy(&sles->R); DSDPCHKERR(info);
00349     info = DSDPVecDestroy(&sles->P); DSDPCHKERR(info);
00350     info = DSDPVecDestroy(&sles->BP); DSDPCHKERR(info);
00351     info = DSDPVecDestroy(&sles->BR); DSDPCHKERR(info);
00352     info = DSDPVecDestroy(&sles->Diag); DSDPCHKERR(info);
00353     info = DSDPVecDestroy(&sles->TTT); DSDPCHKERR(info);
00354   }
00355   DSDPFREE(ssles,&info);DSDPCHKERR(info);
00356   DSDPFunctionReturn(0);
00357 }

Generated on Sun Mar 23 07:30:49 2008 for DSDP by  doxygen 1.5.5