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
00271
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 }