12typedef 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"
157int 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"
312int DSDPCGInitialize(DSDPCG **sles){
316 DSDPCALLOC1(&ssles,DSDPCG,&info);DSDPCHKERR(info);
319 DSDPFunctionReturn(0);
323#define __FUNCT__ "DSDPCGSetup"
324int 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"
342int 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);
Internal data structure for the DSDP solver.
int DSDPHessianMultiplyAdd(DSDP, DSDPVec, DSDPVec)
Add the product of Schur matrix with v.
DSDPTruth
Boolean variables.
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
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 CG method.
int DSDPSchurMatMultiply(DSDPSchurMat M, DSDPVec x, DSDPVec y)
Multiply M by a vector. y = M x.
int DSDPSchurMatSolve(DSDPSchurMat M, DSDPVec b, DSDPVec x)
Solve the linear system.
Methods of a Schur Matrix.
struct DSDPSchurMat_C DSDPSchurMat
This object represents the Schur Matrix. Its structure is opaque to the DSDP solver,...
Error handling, printing, and profiling.
Vector operations used by the solver.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.