DSDP
cholmat2.c
Go to the documentation of this file.
1#include "numchol.h"
2#include "dsdpdualmat_impl.h"
3#include "dsdpsys.h"
4#include "dsdplapack.h"
5
11typedef struct{
12 chfac* spsym;
13 double *sinv;
14 char UPLQ;
15 int n;
16 int dsinv;
17} spmat;
18
19
20static int SMatDestroy(void*S){
21 spmat* SS=(spmat*)S;
22 int info;
23 CfcFree(&SS->spsym);
24 if (SS->dsinv){
25 DSDPFREE(&SS->sinv,&info);
26 }
27 DSDPFREE(&SS,&info);
28 return 0;
29}
30
31static int SMatGetSize(void *S, int *n){
32 spmat* SS=(spmat*)S;
33 *n=SS->spsym->nrow;
34 return 0;
35}
36
37static int SMatView(void* S){
38 spmat* SS=(spmat*)S;
39 int info;
40 info=Mat4View(SS->spsym); DSDPCHKERR(info);
41 return 0;
42}
43
44static int SMatLogDet(void* S, double *dd){
45 spmat* SS=(spmat*)S;
46 int info;
47 info=Mat4LogDet(SS->spsym,dd); DSDPCHKERR(info);
48 return 0;
49}
50
51
52
53static int SMatSetURMatP(void*S, double v[], int nn, int n){
54 spmat* SS=(spmat*)S;
55 int k,j,row,info;
56 double *rw1,*rw2,*xr=v;
57 rw1=SS->spsym->rw;rw2=rw1+n;
58 info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
59 for (k=0;k<n/2;k++){
60 row = 2*k;
61
62 xr=v+row*(row+1)/2;
63 memcpy((void*)rw1,(void*)(xr),(row+1)*sizeof(double));
64 xr+=row+1;
65 rw1[row+1]=xr[row];
66 memcpy((void*)(rw2),(void*)(xr),(row+2)*sizeof(double));
67 xr+=row+2;
68
69 /* memset((void*)rw,0,(n-row)*sizeof(double)); */
70 for (j=row+2;j<n;j++){
71 rw1[j]=xr[row];
72 rw2[j]=xr[row+1];
73 xr+=j+1;
74 }
75
76 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
77 info=MatSetColumn4(SS->spsym,rw2,row+1); DSDPCHKERR(info);
78 }
79
80 for (row=2*(n/2);row<n;row++){
81 /* memset((void*)rw,0,(n-row)*sizeof(double)); */
82 xr=v+row*(row+1)/2;
83 memcpy((void*)(rw1),(void*)(xr),(row+1)*sizeof(double));
84 xr+=row+1;
85 for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=(j+2); }
86 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
87 xr+=(n-row);
88 }
89 return 0;
90}
91
92static int SMatSetURMatU(void*S, double v[], int nn, int n){
93 spmat* SS=(spmat*)S;
94 int k,j,row,info;
95 double *rw1,*rw2,*xr=v;
96 rw1=SS->spsym->rw;rw2=rw1+n;
97 info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
98 for (k=0;k<n/2;k++){
99 row = 2*k;
100
101 xr=v+row*n;
102 memcpy((void*)rw1,(void*)(xr),(row+1)*sizeof(double));
103 xr+=n;
104 rw1[row+1]=xr[row];
105 memcpy((void*)(rw2),(void*)(xr),(row+2)*sizeof(double));
106 xr+=n;
107 /* memset((void*)rw,0,(n-row)*sizeof(double)); */
108 for (j=row+2;j<n;j++){
109 rw1[j]=xr[row];
110 rw2[j]=xr[row+1];
111 xr+=n;
112 }
113
114 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
115 info=MatSetColumn4(SS->spsym,rw2,row+1); DSDPCHKERR(info);
116 }
117
118 for (row=2*(n/2);row<n;row++){
119 /* memset((void*)rw,0,(n-row)*sizeof(double)); */
120 xr=v+row*n;
121 memcpy((void*)(rw1),(void*)(xr),(row+1)*sizeof(double));
122 xr+=n;
123 for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=n; }
124 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
125 }
126 return 0;
127}
128
129static int SMatSetURMat(void*S, double v[], int nn, int n){
130 spmat* SS=(spmat*)S;
131 int info;
132 if (SS->UPLQ=='P'){
133 info=SMatSetURMatP(S,v,nn,n);DSDPCHKERR(info);
134 } else if (SS->UPLQ=='U'){
135 info=SMatSetURMatU(S,v,nn,n);DSDPCHKERR(info);
136 }
137 return 0;
138}
139
140static int SMatSolve(void *S, int indx[], int nind, double b[], double x[],int n){
141 spmat* SS=(spmat*)S;
142 int i,ii;
143 double alpha,*s1=SS->sinv,*s2;
144 ffinteger nn,ione;
145 if (SS->sinv && nind < n/4){
146 memset((void*)x,0,n*sizeof(double));
147 for (i=0;i<nind;i++){
148 ii=indx[i];
149 ione=1;nn=n;alpha=b[ii];s2=s1+n*ii;
150 daxpy(&nn,&alpha,s2,&ione,x,&ione);
151 }
152 } else {
153 memcpy(x,b,n*sizeof(double));
154 ChlSolve(SS->spsym, b, x);
155 }
156 return 0;
157}
158
159static int SMatCholeskySolveBackward(void *S, double b[], double x[],int n){
160 spmat* SS=(spmat*)S;
161 ChlSolveBackward(SS->spsym, b, x);
162 return 0;
163}
164
165static int SMatCholeskySolveForward(void *S, double b[], double x[], int n){
166 spmat* SS=(spmat*)S;
167 ChlSolveForward(SS->spsym, b, x);
168 return 0;
169}
170
171static int SMatFull(void *S, int *full){
172 *full=0;
173 return 0;
174}
175
176static int SMatCholeskyFactor(void *S, int *flag){
177 spmat* SS=(spmat*)S;
178 int *iw;
179 double *rw;
180 cfc_sta Cfact;
181 iw=SS->spsym->iw;
182 rw=SS->spsym->rw;
183 Cfact=(cfc_sta)ChlFact(SS->spsym,iw,rw,TRUE);
184 if (CfcOk!= Cfact){
185 *flag=1;
186 } else {
187 *flag=0;
188 }
189 return 0;
190}
191
192static int SMatInverseAddP(void *S, double alpha, double v[], int nn, int n){
193 spmat* SS=(spmat*)S;
194 ffinteger ii,ione=1;
195 double *x,*b,*ss=SS->sinv;
196 int i,j,k=0;
197
198 if (ss){
199 for (i=0;i<n;i++){
200 v+=i; ii=i+1;
201 daxpy(&ii,&alpha,ss,&ione,v,&ione);
202 ss+=n;
203 }
204 } else {
205 b=SS->spsym->rw;x=b+n;
206 for (i=0;i<n;i++){
207 memset((void*)b,0,n*sizeof(double));
208 b[i]=alpha;
209 ChlSolve(SS->spsym, b, x);
210 k=k+i;
211 for (j=0;j<=i;j++){
212 v[k+j]+=x[j];
213 }
214 }
215 }
216 return 0;
217}
218
219static int SMatInverseAddU(void *S, double alpha, double v[], int nn, int n){
220 spmat* SS=(spmat*)S;
221 ffinteger n2=n*n,ione=1;
222 double *x,*b,*ss=SS->sinv;
223 int i,j,k=0;
224
225 if (ss){
226 daxpy(&n2,&alpha,ss,&ione,v,&ione);
227 } else {
228 b=SS->spsym->rw;x=b+n;
229 for (i=0;i<n;i++){
230 memset((void*)b,0,n*sizeof(double));
231 b[i]=alpha;
232 ChlSolve(SS->spsym, b, x);
233 k=i*n;
234 for (j=0;j<n;j++){
235 v[k+j]+=x[j];
236 }
237 }
238 }
239 return 0;
240}
241
242static int SMatInverseAdd(void *S, double alpha, double v[], int nn, int n){
243 spmat* SS=(spmat*)S;
244 int info;
245 if (SS->UPLQ=='P'){
246 info=SMatInverseAddP(S,alpha,v,nn,n);DSDPCHKERR(info);
247 } else if (SS->UPLQ=='U'){
248 info=SMatInverseAddU(S,alpha,v,nn,n);DSDPCHKERR(info);
249 }
250 return 0;
251}
252
253static int SMatCholeskyForwardMultiply(void *S, double b[], double x[], int n){
254 spmat* SS=(spmat*)S;
255 GetUhat(SS->spsym,b,x);
256 return 0;
257}
258
259static int SMatInvert(void *S){
260 spmat* SS=(spmat*)S;
261 double *x,*b,*v=SS->sinv;
262 int i,n=SS->n;
263
264 b=SS->spsym->rw;x=b+n;
265
266 if (v){
267 for (i=0;i<n;i++){
268 memset((void*)b,0,n*sizeof(double));
269 b[i]=1.0;
270 ChlSolve(SS->spsym, b, x);
271 memcpy((void*)(v+i*n),(void*)x,n*sizeof(double));
272 }
273 }
274 return 0;
275}
276
277static struct DSDPDualMat_Ops sdmatops;
278static const char* tmatname="SPARSE PSD";
279static int SDualOpsInitialize(struct DSDPDualMat_Ops* sops){
280 int info;
281 if (sops==NULL) return 0;
282 info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
283 sops->matcholesky=SMatCholeskyFactor;
284 sops->matsolveforward=SMatCholeskySolveForward;
285 sops->matsolvebackward=SMatCholeskySolveBackward;
286 sops->matinversemultiply=SMatSolve;
287 sops->matinvert=SMatInvert;
288 sops->matinverseadd=SMatInverseAdd;
289 sops->matforwardmultiply=SMatCholeskyForwardMultiply;
290 sops->matseturmat=SMatSetURMat;
291 sops->matfull=SMatFull;
292 sops->matdestroy=SMatDestroy;
293 sops->matgetsize=SMatGetSize;
294 sops->matview=SMatView;
295 sops->matlogdet=SMatLogDet;
296 sops->matname=tmatname;
297 return 0;
298 }
299
300static int dcholmatcreate(int n, char UPLQ, chfac *sp,
301 struct DSDPDualMat_Ops **sops, void**smat){
302 spmat *S;
303 int info;
304 DSDPCALLOC1(&S,spmat,&info);DSDPCHKERR(info);
305 S->UPLQ=UPLQ; S->n=n; S->sinv=0; S->dsinv=0; S->spsym=sp;
306 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
307 *sops=&sdmatops;
308 *smat=(void*)S;
309 return 0;
310 }
311
312static int dcholmatsinverse(int n, spmat *S1, spmat *S2){
313 int info;
314 double *ssinv;
315 DSDPCALLOC2(&ssinv,double,n*n,&info);
316 S1->sinv=ssinv; S2->sinv=ssinv; S2->dsinv=1;
317 return 0;
318}
319
320#undef __FUNCT__
321#define __FUNCT__ "DSDPDenseDualMatCreate"
322int DSDPDenseDualMatCreate(int n, char UPLQ,
323 struct DSDPDualMat_Ops **sops1, void**smat1,
324 struct DSDPDualMat_Ops **sops2, void**smat2){
325 int info=0;
326 chfac *sp;
327
328 DSDPFunctionBegin;
329 info=MchlSetup2(n,&sp); DSDPCHKERR(info);
330 info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
331 info=MchlSetup2(n,&sp); DSDPCHKERR(info);
332 info=dcholmatcreate(n,UPLQ,sp,sops1,smat2);DSDPCHKERR(info);
333 info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
334
335 DSDPFunctionReturn(0);
336}
337
338
339#undef __FUNCT__
340#define __FUNCT__ "DSDPSparseDualMatCreate"
341int DSDPSparseDualMatCreate(int n, int *rnnz, int *snnz,
342 int trank,char UPLQ,int*nnzz,
343 struct DSDPDualMat_Ops **sops1, void**smat1,
344 struct DSDPDualMat_Ops **sops2, void**smat2){
345 int nnz,info=0;
346 chfac *sp;
347
348 DSDPFunctionBegin;
349 SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
350 info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
351 SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
352 info=dcholmatcreate(n,UPLQ,sp,sops2,smat2);DSDPCHKERR(info);
353 nnz=sp->unnz;*nnzz=nnz;
354 if (trank>2*n+2){
355 info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
356 }
357 DSDPFunctionReturn(0);
358}
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Structure of function pointers that each symmetric positive definite matrix type (dense,...
DSDP uses BLAS and LAPACK for many of its operations.
Error handling, printing, and profiling.
Table of function pointers that operate on the S matrix.