12 typedef enum { DSDPNoMatrix=1, DSDPUnfactoredMatrix=2, DSDPFactoredMatrix=3} DSDPCGType;
22 #define __FUNCT__ "DSDPCGMatMult" 26 info=DSDPVecZero(Y); DSDPCHKERR(info);
27 if (M.type==DSDPUnfactoredMatrix){
29 }
else if (M.type==DSDPFactoredMatrix){
30 info=DSDPSchurMatMultR(M.M,X,Y); DSDPCHKERR(info);
31 info=DSDPVecAXPY(-0*M.dsdp->Mshift,X,Y); DSDPCHKERR(info);
32 }
else if (M.type==DSDPNoMatrix){
35 DSDPFunctionReturn(0);
39 #define __FUNCT__ "DSDPCGMatPre" 43 info=DSDPVecZero(Y); DSDPCHKERR(info);
44 if (M.type==DSDPUnfactoredMatrix){
45 info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
46 info=DSDPVecPointwiseMult(Y,M.Diag,Y);DSDPCHKERR(info);
47 }
else if (M.type==DSDPFactoredMatrix){
49 }
else if (M.type==DSDPNoMatrix){
50 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
52 DSDPFunctionReturn(0);
56 #define __FUNCT__ "DSDPCGMatPreLeft" 60 info=DSDPVecZero(Y); DSDPCHKERR(info);
61 if (M.type==DSDPUnfactoredMatrix){
62 info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
63 }
else if (M.type==DSDPFactoredMatrix){
65 }
else if (M.type==DSDPNoMatrix){
66 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
68 DSDPFunctionReturn(0);
72 #define __FUNCT__ "DSDPCGMatPreRight" 76 info=DSDPVecZero(Y); DSDPCHKERR(info);
77 if (M.type==DSDPNoMatrix){
78 info=DSDPVecPointwiseMult(X,M.Diag,Y);DSDPCHKERR(info);
79 }
else if (M.type==DSDPFactoredMatrix){
80 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
81 }
else if (M.type==DSDPUnfactoredMatrix){
82 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
84 DSDPFunctionReturn(0);
89 #define __FUNCT__ "DSDPConjugateResidual" 93 double zero=0.0,minus_one=-1.0;
94 double alpha,beta,bpbp,rBr,rBrOld;
95 double r0,rerr=1.0e20;
98 info=DSDPVecNorm2(X,&rBr); DSDPCHKERR(info);
100 info=DSDPVecCopy(X,P); DSDPCHKERR(info);
101 info=DSDPCGMatPreRight(B,P,X);DSDPCHKERR(info);
102 info=DSDPCGMatMult(B,X,R); DSDPCHKERR(info);
104 info=DSDPVecSet(zero,R); DSDPCHKERR(info);
106 info=DSDPVecAYPX(minus_one,D,R); DSDPCHKERR(info);
108 info=DSDPCGMatPreLeft(B,D,R);DSDPCHKERR(info);
109 info=DSDPVecCopy(R,P); DSDPCHKERR(info);
111 info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
112 info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
113 info=DSDPCGMatPreRight(B,TT3,BR);DSDPCHKERR(info);
115 info=DSDPVecCopy(BR,BP); DSDPCHKERR(info);
116 info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
117 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
120 for (i=0;i<maxits;i++){
122 if (rerr/n < 1.0e-30 || rBr/n <= 1.0e-30 || rBr< 1.0e-12 * r0 )
break;
124 info=DSDPVecDot(BP,BP,&bpbp); DSDPCHKERR(info);
126 info= DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
128 info=DSDPVecAXPY(alpha,BP,R); DSDPCHKERR(info);
130 info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
131 info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
132 info=DSDPCGMatPreLeft(B,TT3,BR);DSDPCHKERR(info);
135 info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
136 info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
138 DSDPLogInfo(0,11,
"CG: rerr: %4.4e, rBr: %4.4e \n",rerr,rBr);
141 info=DSDPVecAYPX(beta,R,P); DSDPCHKERR(info);
142 info=DSDPVecAYPX(beta,BR,BP); DSDPCHKERR(info);
145 info=DSDPVecCopy(X,BR);DSDPCHKERR(info);
146 info=DSDPCGMatPreRight(B,BR,X);DSDPCHKERR(info);
148 DSDPLogInfo(0,2,
"Conjugate Residual, Initial rMr, %4.2e, Final rMr: %4.2e, Iterates: %d\n",r0,rBr,i);
152 DSDPFunctionReturn(0);
156 #define __FUNCT__ "DSDPConjugateGradient" 157 int DSDPConjugateGradient(DSDPCGMat B,
DSDPVec X,
DSDPVec D,
DSDPVec R,
DSDPVec BR,
DSDPVec P,
DSDPVec BP,
DSDPVec Z,
double cgtol,
int maxits,
int *iter){
160 double alpha,beta=0,bpbp;
161 double r0,rerr=1.0e20;
167 DSDPFunctionReturn(0);
169 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
170 info=DSDPVecNorm2(X,&rerr); DSDPCHKERR(info);
171 info=DSDPVecZero(R); DSDPCHKERR(info);
173 info=DSDPCGMatMult(B,X,R);DSDPCHKERR(info);
176 info=DSDPVecAYPX(alpha,D,R); DSDPCHKERR(info);
177 info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
178 if (rerr*sqrt(1.0*n)<1e-11 +0*cgtol*cgtol){
180 DSDPFunctionReturn(0);
184 info=DSDPVecCopy(R,P); DSDPCHKERR(info);
185 info=DSDPVecSetR(P,0);DSDPCHKERR(info);
186 info=DSDPVecNorm2(P,&rerr); DSDPCHKERR(info);
187 info=DSDPCGMatPre(B,R,Z);DSDPCHKERR(info);
190 info=DSDPVecCopy(Z,P); DSDPCHKERR(info);
191 info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
194 for (i=0;i<maxits;i++){
196 info=DSDPCGMatMult(B,P,BP);DSDPCHKERR(info);
197 info=DSDPVecDot(BP,P,&bpbp); DSDPCHKERR(info);
200 info=DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
201 info=DSDPVecAXPY(-alpha,BP,R); DSDPCHKERR(info);
202 info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
204 DSDPLogInfo(0,15,
"CG: iter: %d, rerr: %4.4e, alpha: %4.4e, beta=%4.4e, rz: %4.4e \n",i,rerr,alpha,beta,rz);
206 if (rerr*sqrt(1.0*n) < cgtol){
break;}
207 if (rz*n < cgtol*cgtol){
break;}
208 if (rerr/r0 < cgtol){
break;}
210 if (rerr<=0){
break;}
211 info=DSDPCGMatPre(B,R,Z);DSDPCHKERR(info);
213 info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
215 info= DSDPVecAYPX(beta,Z,P); DSDPCHKERR(info);
217 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);
221 DSDPFunctionReturn(0);
238 #define __FUNCT__ "DSDPCGSolve" 240 int iter=0,n,info,maxit=10;
242 DSDPCG *sles=dsdp->sles;
246 info=DSDPEventLogBegin(dsdp->cgtime);
247 info=DSDPVecZero(X);DSDPCHKERR(info);
248 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
252 }
else if (dsdp->slestype==1){
254 CGM.type=DSDPNoMatrix;
260 }
else if (dsdp->slestype==2){
261 CGM.type=DSDPUnfactoredMatrix;
266 maxit=(int)sqrt(1.0*n)+10;
267 if (maxit>20) maxit=20;
268 info=DSDPVecSet(1.0,sles->Diag);DSDPCHKERR(info);
274 }
else if (dsdp->slestype==3){
275 CGM.type=DSDPFactoredMatrix;
280 if (ymax > 1e5 && dsdp->rgap<1e-1) maxit=3;
281 if (0 && ymax > 1e5 && dsdp->rgap<1e-2){
283 }
else if (dsdp->rgap<1e-5){
289 }
else if (dsdp->slestype==4){
290 CGM.type=DSDPFactoredMatrix;
296 if (n<maxit) maxit=n;
298 info=DSDPConjugateGradient(CGM,X,RHS,
299 sles->R,sles->BR,sles->P,sles->BP,
300 sles->TTT,cgtol,maxit,&iter);DSDPCHKERR(info);
303 info=DSDPVecDot(RHS,X,&dd);DSDPCHKERR(info);
305 info=DSDPEventLogEnd(dsdp->cgtime);
306 DSDPFunctionReturn(0);
311 #define __FUNCT__ "DSDPCGInitialize" 312 int DSDPCGInitialize(DSDPCG **sles){
316 DSDPCALLOC1(&ssles,DSDPCG,&info);DSDPCHKERR(info);
319 DSDPFunctionReturn(0);
323 #define __FUNCT__ "DSDPCGSetup" 324 int DSDPCGSetup(DSDPCG *sles,
DSDPVec X){
327 info=DSDPVecGetSize(X,&sles->m);DSDPCHKERR(info);
328 if (sles->setup2==0){
329 info = DSDPVecDuplicate(X,&sles->R); DSDPCHKERR(info);
330 info = DSDPVecDuplicate(X,&sles->P); DSDPCHKERR(info);
331 info = DSDPVecDuplicate(X,&sles->BP); DSDPCHKERR(info);
332 info = DSDPVecDuplicate(X,&sles->BR); DSDPCHKERR(info);
333 info = DSDPVecDuplicate(X,&sles->Diag); DSDPCHKERR(info);
334 info = DSDPVecDuplicate(X,&sles->TTT); DSDPCHKERR(info);
337 DSDPFunctionReturn(0);
341 #define __FUNCT__ "DSDPCGDestroy" 342 int DSDPCGDestroy(DSDPCG **ssles){
346 if (ssles==0 || sles==0){DSDPFunctionReturn(0);}
347 if (sles->setup2==1){
348 info = DSDPVecDestroy(&sles->R); DSDPCHKERR(info);
349 info = DSDPVecDestroy(&sles->P); DSDPCHKERR(info);
350 info = DSDPVecDestroy(&sles->BP); DSDPCHKERR(info);
351 info = DSDPVecDestroy(&sles->BR); DSDPCHKERR(info);
352 info = DSDPVecDestroy(&sles->Diag); DSDPCHKERR(info);
353 info = DSDPVecDestroy(&sles->TTT); DSDPCHKERR(info);
355 DSDPFREE(ssles,&info);DSDPCHKERR(info);
356 DSDPFunctionReturn(0);
DSDPTruth
Boolean variables.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Schur complement matrix whose solution is the Newton direction.
int DSDPSchurMatSolve(DSDPSchurMat M, DSDPVec b, DSDPVec x)
Solve the linear system.
Error handling, printing, and profiling.
Internal structures for the DSDP solver.
int DSDPSchurMatMultiply(DSDPSchurMat M, DSDPVec x, DSDPVec y)
Multiply M by a vector. y = M x.
int DSDPCGSolve(DSDP dsdp, DSDPSchurMat MM, DSDPVec RHS, DSDPVec X, double cgtol, DSDPTruth *success)
Apply CG to solve for the step directions.
Internal data structure for the DSDP solver.
Internal data structure for CG method.
Methods of a Schur Matrix.
Vector operations used by the solver.
int DSDPHessianMultiplyAdd(DSDP, DSDPVec, DSDPVec)
Add the product of Schur matrix with v.
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.