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