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
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.