DSDP
diag.c
Go to the documentation of this file.
1 
2 #include "dsdpschurmat_impl.h"
3 #include "dsdpdualmat_impl.h"
4 #include "dsdpdsmat_impl.h"
5 #include "dsdpsys.h"
12 typedef struct {
13  int n;
14  double *val;
15  int owndata;
16 } diagmat;
17 
18 static int DiagMatCreate(int,diagmat**);
19 static int DiagMatMult(void*,double[],double[],int);
20 static int DiagMatGetSize(void*, int *);
21 static int DiagMatAddRow2(void*, int, double, double[], int);
22 static int DiagMatDestroy(void*);
23 static int DiagMatView(void*);
24 static int DiagMatLogDeterminant(void*, double *);
25 
26 /* static int DiagMatScale(double *, diagmat *); */
27 
28 static int DiagMatCreate(int n,diagmat**M){
29  int info;
30  diagmat *M7;
31 
32  DSDPCALLOC1(&M7,diagmat,&info);DSDPCHKERR(info);
33  DSDPCALLOC2(&M7->val,double,n,&info);DSDPCHKERR(info);
34  if (n>0 && M7->val==NULL) return 1;
35  M7->owndata=1; M7->n=n;
36  *M=M7;
37  return 0;
38 }
39 
40 static int DiagMatDestroy(void* AA){
41  int info;
42  diagmat* A=(diagmat*) AA;
43  if (A->owndata && A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
44  DSDPFREE(&A,&info);DSDPCHKERR(info);
45  return 0;
46 }
47 
48 
49 static int DiagMatMult(void* AA, double x[], double y[], int n){
50 
51  diagmat* A=(diagmat*) AA;
52  double *vv=A->val;
53  int i;
54 
55  if (A->n != n) return 1;
56  if (x==0 && n>0) return 3;
57  if (y==0 && n>0) return 3;
58 
59  for (i=0; i<n; i++){
60  y[i]=x[i]*vv[i];
61  }
62  return 0;
63 }
64 
65 
66 static int DiagMatGetSize(void *AA, int *n){
67  diagmat* A=(diagmat*) AA;
68  *n=A->n;
69  return 0;
70 }
71 
72 static int DiagMatView(void* AA){
73  diagmat* A=(diagmat*) AA;
74  int i;
75  for (i=0;i<A->n; i++){
76  printf(" Row: %d, Column: %d, Value: %8.4e \n",i,i,A->val[i]);
77  }
78  return 0;
79 }
80 
81 static int DiagMatAddRow2(void* AA, int nrow, double dd, double row[], int n){
82  diagmat* A=(diagmat*) AA;
83  A->val[nrow] += dd*row[nrow];
84  return 0;
85 }
86 
87 
88 static int DiagMatAddElement(void*A, int nrow, double dd){
89  diagmat* AA = (diagmat*)A;
90  AA->val[nrow] += dd;
91  return 0;
92 }
93 
94 static int DiagMatZeros(void*A){
95  diagmat* AA = (diagmat*)A;
96  int n=AA->n;
97  memset(AA->val,0,n*sizeof(double));
98  return 0;
99 }
100 
101 static int DiagMatSolve(void* A, double b[], double x[],int n){
102  diagmat* AA = (diagmat*)A;
103  double *v=AA->val;
104  int i;
105  for (i=0;i<n;i++){
106  x[i]=b[i]/v[i];
107  }
108  return 0;
109 }
110 
111 static int DiagMatSolve2(void* A, int indx[], int nindx, double b[], double x[],int n){
112  diagmat* AA = (diagmat*)A;
113  double *v=AA->val;
114  int i,j;
115  memset((void*)x,0,n*sizeof(double));
116  for (j=0;j<nindx;j++){
117  i=indx[j];
118  x[i]=b[i]/v[i];
119  }
120  return 0;
121 }
122 
123 static int DiagMatCholeskySolveBackward(void* A, double b[], double x[],int n){
124  int i;
125  for (i=0;i<n;i++){
126  x[i]=b[i];
127  }
128  return 0;
129 }
130 
131 static int DiagMatInvert(void *A){
132  return 0;
133 }
134 
135 static int DiagMatCholeskyFactor(void*A,int *flag){
136  diagmat* AA = (diagmat*)A;
137  double *v=AA->val;
138  int i,n=AA->n;
139  *flag=0;
140  for (i=0;i<n;i++){
141  if (v[i]<=0){ *flag=i+1; break;}
142  }
143  return 0;
144 }
145 
146 static int DiagMatLogDeterminant(void*A, double *dd){
147  diagmat* AA = (diagmat*)A;
148  double d=0,*val=AA->val;
149  int i,n=AA->n;
150  for (i=0;i<n;i++){
151  if (val[i]<=0) return 1;
152  d+=log(val[i]);
153  }
154  *dd=d;
155  return 0;
156 }
157 
158 static int DiagMatTakeUREntriesP(void*A, double dd[], int nn, int n){
159  diagmat* AA = (diagmat*)A;
160  int i,ii;
161  double *val=AA->val;
162  for (i=0;i<n;i++){
163  ii=(i+1)*(i+2)/2-1;
164  val[i]=dd[ii];
165  }
166  return 0;
167 }
168 static int DiagMatTakeUREntriesU(void*A, double dd[], int nn, int n){
169  diagmat* AA = (diagmat*)A;
170  int i;
171  double *val=AA->val;
172  for (i=0;i<n;i++){
173  val[i]=dd[i*n+i];
174  }
175  return 0;
176 }
177 
178 static int DiagMatInverseAddP(void*A, double alpha, double dd[], int nn, int n){
179  diagmat* AA = (diagmat*)A;
180  int i,ii;
181  double *val=AA->val;
182  for (i=0;i<n;i++){
183  ii=(i+1)*(i+2)/2-1;
184  dd[ii]+=alpha/val[i];
185  }
186  return 0;
187 }
188 static int DiagMatInverseAddU(void*A, double alpha, double dd[], int nn, int n){
189  diagmat* AA = (diagmat*)A;
190  int i;
191  double *val=AA->val;
192  for (i=0;i<n;i++){
193  dd[i*n+i]+=alpha/val[i];
194  }
195  return 0;
196 }
197 
198 static int DiagMatFull(void*A, int* dfull){
199  *dfull=1;
200  return 0;
201 }
202 
203 static struct DSDPDualMat_Ops sdmatopsp;
204 static struct DSDPDualMat_Ops sdmatopsu;
205 static const char* diagmatname="DIAGONAL";
206 
207 static int DiagDualOpsInitializeP(struct DSDPDualMat_Ops* sops){
208  int info;
209  if (sops==NULL) return 0;
210  info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
211  sops->matcholesky=DiagMatCholeskyFactor;
212  sops->matsolveforward=DiagMatSolve;
213  sops->matsolvebackward=DiagMatCholeskySolveBackward;
214  sops->matinvert=DiagMatInvert;
215  sops->matinverseadd=DiagMatInverseAddP;
216  sops->matinversemultiply=DiagMatSolve2;
217  sops->matseturmat=DiagMatTakeUREntriesP;
218  sops->matfull=DiagMatFull;
219  sops->matdestroy=DiagMatDestroy;
220  sops->matgetsize=DiagMatGetSize;
221  sops->matview=DiagMatView;
222  sops->matlogdet=DiagMatLogDeterminant;
223  sops->id=9;
224  sops->matname=diagmatname;
225  return 0;
226 }
227 static int DiagDualOpsInitializeU(struct DSDPDualMat_Ops* sops){
228  int info;
229  if (sops==NULL) return 0;
230  info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
231  sops->matcholesky=DiagMatCholeskyFactor;
232  sops->matsolveforward=DiagMatSolve;
233  sops->matsolvebackward=DiagMatCholeskySolveBackward;
234  sops->matinvert=DiagMatInvert;
235  sops->matinversemultiply=DiagMatSolve2;
236  sops->matseturmat=DiagMatTakeUREntriesU;
237  sops->matfull=DiagMatFull;
238  sops->matinverseadd=DiagMatInverseAddU;
239  sops->matdestroy=DiagMatDestroy;
240  sops->matgetsize=DiagMatGetSize;
241  sops->matview=DiagMatView;
242  sops->matlogdet=DiagMatLogDeterminant;
243  sops->id=9;
244  sops->matname=diagmatname;
245  return 0;
246 }
247 
248 #undef __FUNCT__
249 #define __FUNCT__ "DSDPDiagDualMatCreateP"
250 int DSDPDiagDualMatCreateP(int n,
251  struct DSDPDualMat_Ops **sops1, void**smat1,
252  struct DSDPDualMat_Ops **sops2, void**smat2){
253  diagmat *AA;
254  int info;
255  DSDPFunctionBegin;
256 
257  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
258  info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
259  *sops1=&sdmatopsp;
260  *smat1=(void*)AA;
261 
262  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
263  info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
264  *sops2=&sdmatopsp;
265  *smat2=(void*)AA;
266  DSDPFunctionReturn(0);
267 }
268 
269 #undef __FUNCT__
270 #define __FUNCT__ "DSDPDiagDualMatCreateU"
271 int DSDPDiagDualMatCreateU(int n,
272  struct DSDPDualMat_Ops **sops1, void**smat1,
273  struct DSDPDualMat_Ops **sops2, void**smat2){
274  diagmat *AA;
275  int info;
276  DSDPFunctionBegin;
277  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
278  info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
279  *sops1=&sdmatopsu;
280  *smat1=(void*)AA;
281  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
282  info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
283  *sops2=&sdmatopsu;
284  *smat2=(void*)AA;
285  DSDPFunctionReturn(0);
286 }
287 
288 
289 
290 static int DiagMatVecVec(void*A,double x[], int n, double *vv){
291  diagmat* AA = (diagmat*)A;
292  double *v=AA->val,vAv=0;
293  int i;
294  for (i=0;i<n;i++){
295  vAv+=x[i]*x[i]*v[i];
296  }
297  *vv=vAv;
298  return 0;
299 }
300 
301 static int DDiagDSMatOpsInitP(struct DSDPDSMat_Ops *ddiagops){
302  int info;
303  if (ddiagops==NULL) return 0;
304  info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
305  ddiagops->matseturmat=DiagMatTakeUREntriesP;
306  ddiagops->matview=DiagMatView;
307  ddiagops->matgetsize=DiagMatGetSize;
308  ddiagops->matmult=DiagMatMult;
309  ddiagops->matvecvec=DiagMatVecVec;
310  ddiagops->matzeroentries=DiagMatZeros;
311  ddiagops->matdestroy=DiagMatDestroy;
312  ddiagops->id=9;
313  ddiagops->matname=diagmatname;
314  DSDPFunctionReturn(0);
315 }
316 static int DDiagDSMatOpsInitU(struct DSDPDSMat_Ops *ddiagops){
317  int info;
318  if (ddiagops==NULL) return 0;
319  info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
320  ddiagops->matseturmat=DiagMatTakeUREntriesU;
321  ddiagops->matview=DiagMatView;
322  ddiagops->matgetsize=DiagMatGetSize;
323  ddiagops->matmult=DiagMatMult;
324  ddiagops->matvecvec=DiagMatVecVec;
325  ddiagops->matzeroentries=DiagMatZeros;
326  ddiagops->matdestroy=DiagMatDestroy;
327  ddiagops->id=9;
328  ddiagops->matname=diagmatname;
329  DSDPFunctionReturn(0);
330 }
331 
332 static struct DSDPDSMat_Ops dsdiagmatopsp;
333 static struct DSDPDSMat_Ops dsdiagmatopsu;
334 
335 #undef __FUNCT__
336 #define __FUNCT__ "DSDPDiagDSMatP"
337 int DSDPCreateDiagDSMatP(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
338 
339  int info=0;
340  diagmat *AA;
341 
342  DSDPFunctionBegin;
343  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
344  info=DDiagDSMatOpsInitP(&dsdiagmatopsp); DSDPCHKERR(info);
345  *dsmatops=&dsdiagmatopsp;
346  *dsmat=(void*)AA;
347  DSDPFunctionReturn(0);
348 }
349 #undef __FUNCT__
350 #define __FUNCT__ "DSDPDiagDSMatU"
351 int DSDPCreateDiagDSMatU(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
352 
353  int info=0;
354  diagmat *AA;
355 
356  DSDPFunctionBegin;
357  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
358  info=DDiagDSMatOpsInitU(&dsdiagmatopsu); DSDPCHKERR(info);
359  *dsmatops=&dsdiagmatopsu;
360  *dsmat=(void*)AA;
361  DSDPFunctionReturn(0);
362 }
363 
364 static int DiagRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
365  DSDPFunctionBegin;
366  *ncols = 1;
367  cols[row]=1.0;
368  DSDPFunctionReturn(0);
369 }
370 
371 static int DiagAddDiag(void*M, double diag[], int m){
372  diagmat* AA = (diagmat*)M;
373  double *v=AA->val;
374  int i;
375  DSDPFunctionBegin;
376  for (i=0;i<m;i++){
377  v[i]+=diag[i];
378  }
379  DSDPFunctionReturn(0);
380 }
381 
382 static int DiagMultiply(void*M, double xin[], double xout[], int m){
383  diagmat* AA = (diagmat*)M;
384  double *v=AA->val;
385  int i;
386  DSDPFunctionBegin;
387  for (i=0;i<m;i++){
388  xout[i]+=v[i]*xin[i];
389  }
390  DSDPFunctionReturn(0);
391 }
392 
393 static int DiagShiftDiag(void*M, double dd){
394  diagmat* AA = (diagmat*)M;
395  double *v=AA->val;
396  int i,m=AA->n;
397  DSDPFunctionBegin;
398  for (i=0;i<m;i++){
399  v[i]+=dd;
400  }
401  DSDPFunctionReturn(0);
402 }
403 
404 static int DiagAddElement(void*M, int ii, double dd){
405  diagmat* AA = (diagmat*)M;
406  DSDPFunctionBegin;
407  AA->val[ii]+=dd;
408  DSDPFunctionReturn(0);
409 }
410 
411 static int DiagMatOnProcessor(void*A,int row,int*flag){
412  *flag=1;
413  return 0;
414 }
415 
416 static int DiagAssemble(void*M){
417  return 0;
418 }
419 
420 static struct DSDPSchurMat_Ops dsdpdiagschurops;
421 
422 #undef __FUNCT__
423 #define __FUNCT__ "DSDPDiagSchurOps"
424 static int DiagSchurOps(struct DSDPSchurMat_Ops *sops){
425  int info;
426  DSDPFunctionBegin;
427  if (!sops) return 0;
428  info=DSDPSchurMatOpsInitialize(sops); DSDPCHKERR(info);
429  sops->matzero=DiagMatZeros;
430  sops->mataddrow=DiagMatAddRow2;
431  sops->mataddelement=DiagMatAddElement;
432  sops->matdestroy=DiagMatDestroy;
433  sops->matfactor=DiagMatCholeskyFactor;
434  sops->matsolve=DiagMatSolve;
435  sops->matadddiagonal=DiagAddDiag;
436  sops->matshiftdiagonal=DiagShiftDiag;
437  sops->mataddelement=DiagAddElement;
438  sops->matscaledmultiply=DiagMultiply;
439  sops->matassemble=DiagAssemble;
440  sops->pmatonprocessor=DiagMatOnProcessor;
441  sops->matrownonzeros=DiagRowNonzeros;
442  sops->id=9;
443  sops->matname=diagmatname;
444  DSDPFunctionReturn(0);
445 }
446 
447 #undef __FUNCT__
448 #define __FUNCT__ "DSDPGetDiagSchurMat"
449 int DSDPGetDiagSchurMat(int n, struct DSDPSchurMat_Ops **sops, void **data){
450  int info=0;
451  diagmat *AA;
452  DSDPFunctionBegin;
453  info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
454  info=DiagSchurOps(&dsdpdiagschurops); DSDPCHKERR(info);
455  if (sops){*sops=&dsdpdiagschurops;}
456  if (data){*data=(void*)AA;}
457  DSDPFunctionReturn(0);
458 }
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Definition: dsdpdualmat.c:423
Error handling, printing, and profiling.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide...
Structure of function pointers that each SDP Delta S matrix type (sparse, dense, diagonal, ...) must implement.
Table of function pointers that operate on the S matrix.
int DSDPDSMatOpsInitialize(struct DSDPDSMat_Ops *aops)
Set pointers to null.
Definition: dsdpdsmat.c:214
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Definition: dsdpschurmat.c:44
Structure of function pointers that each symmetric positive definite matrix type (dense, sparse) must implement.
Symmetric Delta S matrix for one block in the semidefinite cone.