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"
12 typedef enum { DSDPNoMatrix=1, DSDPUnfactoredMatrix=2, DSDPFactoredMatrix=3} DSDPCGType;
13 
14 typedef struct{
15  DSDPCGType type;
16  DSDPSchurMat M;
17  DSDPVec Diag;
18  DSDP dsdp;
19 } DSDPCGMat;
20 
21 #undef __FUNCT__
22 #define __FUNCT__ "DSDPCGMatMult"
23 int 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"
40 int 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"
57 int 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"
73 int 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"
90 int 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"
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){
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"
239 int 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"
312 int 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"
324 int 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"
342 int 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 }
DSDPTruth
Boolean variables.
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.
Definition: dsdpschurmat.h:35
int DSDPSchurMatSolve(DSDPSchurMat M, DSDPVec b, DSDPVec x)
Solve the linear system.
Definition: dsdpschurmat.c:466
Error handling, printing, and profiling.
Internal structures for the DSDP solver.
Definition: dsdp.h:65
int DSDPSchurMatMultiply(DSDPSchurMat M, DSDPVec x, DSDPVec y)
Multiply M by a vector. y = M x.
Definition: dsdpschurmat.c:231
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 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.
Definition: dsdpcops.c:188
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.
Definition: dsdpsetdata.c:645