DSDP
dsdpcg.c
Go to the documentation of this file.
1#include "dsdpcg.h"
2#include "dsdpschurmat.h"
3#include "dsdpvec.h"
4#include "dsdpsys.h"
5#include "dsdp.h"
12typedef enum { DSDPNoMatrix=1, DSDPUnfactoredMatrix=2, DSDPFactoredMatrix=3} DSDPCGType;
13
14typedef struct{
15 DSDPCGType type;
17 DSDPVec Diag;
18 DSDP dsdp;
19} DSDPCGMat;
20
21#undef __FUNCT__
22#define __FUNCT__ "DSDPCGMatMult"
23int DSDPCGMatMult(DSDPCGMat M, DSDPVec X, DSDPVec Y){
24 int info;
25 DSDPFunctionBegin;
26 info=DSDPVecZero(Y); DSDPCHKERR(info);
27 if (M.type==DSDPUnfactoredMatrix){
28 info=DSDPSchurMatMultiply(M.M,X,Y); DSDPCHKERR(info);
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){
33 info=DSDPHessianMultiplyAdd(M.dsdp,X,Y);DSDPCHKERR(info);
34 }
35 DSDPFunctionReturn(0);
36}
37
38#undef __FUNCT__
39#define __FUNCT__ "DSDPCGMatPre"
40int DSDPCGMatPre(DSDPCGMat M, DSDPVec X, DSDPVec Y){
41 int info;
42 DSDPFunctionBegin;
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){
48 info=DSDPSchurMatSolve(M.M,X,Y);DSDPCHKERR(info);
49 } else if (M.type==DSDPNoMatrix){
50 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
51 }
52 DSDPFunctionReturn(0);
53}
54
55#undef __FUNCT__
56#define __FUNCT__ "DSDPCGMatPreLeft"
57int DSDPCGMatPreLeft(DSDPCGMat M, DSDPVec X, DSDPVec Y){
58 int info;
59 DSDPFunctionBegin;
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){
64 info=DSDPSchurMatSolve(M.M,X,Y);DSDPCHKERR(info);
65 } else if (M.type==DSDPNoMatrix){
66 info=DSDPVecCopy(X,Y);DSDPCHKERR(info);
67 }
68 DSDPFunctionReturn(0);
69}
70
71#undef __FUNCT__
72#define __FUNCT__ "DSDPCGMatPreRight"
73int DSDPCGMatPreRight(DSDPCGMat M, DSDPVec X, DSDPVec Y){
74 int info;
75 DSDPFunctionBegin;
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);
83 }
84 DSDPFunctionReturn(0);
85}
86
87
88#undef __FUNCT__
89#define __FUNCT__ "DSDPConjugateResidual"
90int DSDPConjugateResidual(DSDPCGMat B, DSDPVec X, DSDPVec D, DSDPVec R, DSDPVec BR, DSDPVec P, DSDPVec BP, DSDPVec TT3, int maxits, int *iter){
91
92 int i,n,info;
93 double zero=0.0,minus_one=-1.0;
94 double alpha,beta,bpbp,rBr,rBrOld;
95 double r0,rerr=1.0e20;
96
97 DSDPFunctionBegin;
98 info=DSDPVecNorm2(X,&rBr); DSDPCHKERR(info);
99 if (rBr>0){
100 info=DSDPVecCopy(X,P); DSDPCHKERR(info);
101 info=DSDPCGMatPreRight(B,P,X);DSDPCHKERR(info);
102 info=DSDPCGMatMult(B,X,R); DSDPCHKERR(info);
103 } else {
104 info=DSDPVecSet(zero,R); DSDPCHKERR(info);
105 }
106 info=DSDPVecAYPX(minus_one,D,R); DSDPCHKERR(info);
107
108 info=DSDPCGMatPreLeft(B,D,R);DSDPCHKERR(info);
109 info=DSDPVecCopy(R,P); DSDPCHKERR(info);
110
111 info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
112 info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
113 info=DSDPCGMatPreRight(B,TT3,BR);DSDPCHKERR(info);
114
115 info=DSDPVecCopy(BR,BP); DSDPCHKERR(info);
116 info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
117 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
118 r0=rBr;
119
120 for (i=0;i<maxits;i++){
121
122 if (rerr/n < 1.0e-30 || rBr/n <= 1.0e-30 || rBr< 1.0e-12 * r0 ) break;
123
124 info=DSDPVecDot(BP,BP,&bpbp); DSDPCHKERR(info);
125 alpha = rBr / bpbp;
126 info= DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
127 alpha = -alpha;
128 info=DSDPVecAXPY(alpha,BP,R); DSDPCHKERR(info);
129
130 info=DSDPCGMatPreRight(B,R,BR);DSDPCHKERR(info);
131 info=DSDPCGMatMult(B,BR,TT3); DSDPCHKERR(info);
132 info=DSDPCGMatPreLeft(B,TT3,BR);DSDPCHKERR(info);
133
134 rBrOld=rBr;
135 info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
136 info=DSDPVecDot(BR,R,&rBr); DSDPCHKERR(info);
137
138 DSDPLogInfo(0,11,"CG: rerr: %4.4e, rBr: %4.4e \n",rerr,rBr);
139
140 beta = rBr/rBrOld;
141 info=DSDPVecAYPX(beta,R,P); DSDPCHKERR(info);
142 info=DSDPVecAYPX(beta,BR,BP); DSDPCHKERR(info);
143
144 }
145 info=DSDPVecCopy(X,BR);DSDPCHKERR(info);
146 info=DSDPCGMatPreRight(B,BR,X);DSDPCHKERR(info);
147
148 DSDPLogInfo(0,2,"Conjugate Residual, Initial rMr, %4.2e, Final rMr: %4.2e, Iterates: %d\n",r0,rBr,i);
149
150 *iter=i;
151
152 DSDPFunctionReturn(0);
153}
154
155#undef __FUNCT__
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){
158
159 int i,n,info;
160 double alpha,beta=0,bpbp;
161 double r0,rerr=1.0e20;
162 double rz,rzold,rz0;
163
164 DSDPFunctionBegin;
165 if (maxits<=0){
166 *iter=0;
167 DSDPFunctionReturn(0);
168 }
169 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
170 info=DSDPVecNorm2(X,&rerr); DSDPCHKERR(info);
171 info=DSDPVecZero(R); DSDPCHKERR(info);
172 if (rerr>0){
173 info=DSDPCGMatMult(B,X,R);DSDPCHKERR(info);
174 }
175 alpha=-1.0;
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){
179 *iter=1;
180 DSDPFunctionReturn(0);
181 }
182
183 if (rerr>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);
188 }
189
190 info=DSDPVecCopy(Z,P); DSDPCHKERR(info);
191 info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
192 r0=rerr;rz0=rz;
193
194 for (i=0;i<maxits;i++){
195
196 info=DSDPCGMatMult(B,P,BP);DSDPCHKERR(info);
197 info=DSDPVecDot(BP,P,&bpbp); DSDPCHKERR(info);
198 if (bpbp==0) break;
199 alpha = rz/bpbp;
200 info=DSDPVecAXPY(alpha,P,X); DSDPCHKERR(info);
201 info=DSDPVecAXPY(-alpha,BP,R); DSDPCHKERR(info);
202 info=DSDPVecNorm2(R,&rerr); DSDPCHKERR(info);
203
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);
205 if (i>1){
206 if (rerr*sqrt(1.0*n) < cgtol){ break;}
207 if (rz*n < cgtol*cgtol){ break;}
208 if (rerr/r0 < cgtol){ break;}
209 }
210 if (rerr<=0){ break;}
211 info=DSDPCGMatPre(B,R,Z);DSDPCHKERR(info);
212 rzold=rz;
213 info=DSDPVecDot(R,Z,&rz); DSDPCHKERR(info);
214 beta=rz/rzold;
215 info= DSDPVecAYPX(beta,Z,P); DSDPCHKERR(info);
216 }
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);
218
219 *iter=i+1;
220
221 DSDPFunctionReturn(0);
222}
223
237#undef __FUNCT__
238#define __FUNCT__ "DSDPCGSolve"
239int DSDPCGSolve(DSDP dsdp, DSDPSchurMat MM, DSDPVec RHS, DSDPVec X,double cgtol, DSDPTruth *success){
240 int iter=0,n,info,maxit=10;
241 double dd,ymax;
242 DSDPCG *sles=dsdp->sles;
243 DSDPCGMat CGM;
244 DSDPFunctionBegin;
245
246 info=DSDPEventLogBegin(dsdp->cgtime);
247 info=DSDPVecZero(X);DSDPCHKERR(info);
248 info=DSDPVecGetSize(X,&n); DSDPCHKERR(info);
249 *success=DSDP_TRUE;
250 if (0){
251 maxit=0;
252 } else if (dsdp->slestype==1){
253
254 CGM.type=DSDPNoMatrix;
255 CGM.M=MM;
256 CGM.dsdp=dsdp;
257 cgtol*=1000;
258 maxit=5;
259
260 } else if (dsdp->slestype==2){
261 CGM.type=DSDPUnfactoredMatrix;
262 CGM.M=MM;
263 CGM.Diag=sles->Diag;
264 CGM.dsdp=dsdp;
265 cgtol*=100;
266 maxit=(int)sqrt(1.0*n)+10;
267 if (maxit>20) maxit=20;
268 info=DSDPVecSet(1.0,sles->Diag);DSDPCHKERR(info);
269 /*
270 info=DSDPSchurMatGetDiagonal(MM,sles->Diag);DSDPCHKERR(info);
271 info=DSDPVecReciprocalSqrt(sles->Diag); DSDPCHKERR(info);
272 */
273
274 } else if (dsdp->slestype==3){
275 CGM.type=DSDPFactoredMatrix;
276 CGM.M=MM;
277 CGM.dsdp=dsdp;
278 maxit=0;
279 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
280 if (ymax > 1e5 && dsdp->rgap<1e-1) maxit=3;
281 if (0 && ymax > 1e5 && dsdp->rgap<1e-2){
282 maxit=6;
283 } else if (dsdp->rgap<1e-5){
284 maxit=3;
285 }
286
287 info=DSDPSchurMatSolve(MM,RHS,X);DSDPCHKERR(info);
288
289 } else if (dsdp->slestype==4){
290 CGM.type=DSDPFactoredMatrix;
291 CGM.M=MM;
292 CGM.dsdp=dsdp;
293 maxit=3;
294 info=DSDPSchurMatSolve(MM,RHS,X);DSDPCHKERR(info);
295 }
296 if (n<maxit) maxit=n;
297
298 info=DSDPConjugateGradient(CGM,X,RHS,
299 sles->R,sles->BR,sles->P,sles->BP,
300 sles->TTT,cgtol,maxit,&iter);DSDPCHKERR(info);
301
302 if (iter>=maxit) *success=DSDP_FALSE;
303 info=DSDPVecDot(RHS,X,&dd);DSDPCHKERR(info);
304 if (dd<0) *success=DSDP_FALSE;
305 info=DSDPEventLogEnd(dsdp->cgtime);
306 DSDPFunctionReturn(0);
307}
308
309
310#undef __FUNCT__
311#define __FUNCT__ "DSDPCGInitialize"
312int DSDPCGInitialize(DSDPCG **sles){
313 DSDPCG *ssles;
314 int info;
315 DSDPFunctionBegin;
316 DSDPCALLOC1(&ssles,DSDPCG,&info);DSDPCHKERR(info);
317 ssles->setup2=0;
318 *sles=ssles;
319 DSDPFunctionReturn(0);
320}
321
322#undef __FUNCT__
323#define __FUNCT__ "DSDPCGSetup"
324int DSDPCGSetup(DSDPCG *sles,DSDPVec X){
325 int info;
326 DSDPFunctionBegin;
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);
335 }
336 sles->setup2=1;
337 DSDPFunctionReturn(0);
338}
339
340#undef __FUNCT__
341#define __FUNCT__ "DSDPCGDestroy"
342int DSDPCGDestroy(DSDPCG **ssles){
343 int info;
344 DSDPCG *sles=*ssles;
345 DSDPFunctionBegin;
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);
354 }
355 DSDPFREE(ssles,&info);DSDPCHKERR(info);
356 DSDPFunctionReturn(0);
357}
Internal data structure for the DSDP solver.
int DSDPHessianMultiplyAdd(DSDP, DSDPVec, DSDPVec)
Add the product of Schur matrix with v.
Definition dsdpcops.c:188
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
int DSDPCGSolve(DSDP dsdp, DSDPSchurMat MM, DSDPVec RHS, DSDPVec X, double cgtol, DSDPTruth *success)
Apply CG to solve for the step directions.
Definition dsdpcg.c:239
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.
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.
Definition dsdpvec.h:25
Schur complement matrix whose solution is the Newton direction.
Internal structures for the DSDP solver.
Definition dsdp.h:65