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 
11 typedef struct{
12  chfac* spsym;
13  double *sinv;
14  char UPLQ;
15  int n;
16  int dsinv;
17 } spmat;
18 
19 
20 static 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 
31 static int SMatGetSize(void *S, int *n){
32  spmat* SS=(spmat*)S;
33  *n=SS->spsym->nrow;
34  return 0;
35 }
36 
37 static int SMatView(void* S){
38  spmat* SS=(spmat*)S;
39  int info;
40  info=Mat4View(SS->spsym); DSDPCHKERR(info);
41  return 0;
42 }
43 
44 static 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 
53 static 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 
92 static 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 
129 static 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 
140 static 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 
159 static 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 
165 static 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 
171 static int SMatFull(void *S, int *full){
172  *full=0;
173  return 0;
174 }
175 
176 static 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 
192 static 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 
219 static 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 
242 static 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 
253 static 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 
259 static 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 
277 static struct DSDPDualMat_Ops sdmatops;
278 static const char* tmatname="SPARSE PSD";
279 static 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 
300 static 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 
312 static 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"
322 int 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"
341 int 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.
Definition: dsdpdualmat.c:423
Error handling, printing, and profiling.
Table of function pointers that operate on the S matrix.
Structure of function pointers that each symmetric positive definite matrix type (dense, sparse) must implement.
DSDP uses BLAS and LAPACK for many of its operations.