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"
12typedef struct {
13 int n;
14 double *val;
15 int owndata;
16} diagmat;
17
18static int DiagMatCreate(int,diagmat**);
19static int DiagMatMult(void*,double[],double[],int);
20static int DiagMatGetSize(void*, int *);
21static int DiagMatAddRow2(void*, int, double, double[], int);
22static int DiagMatDestroy(void*);
23static int DiagMatView(void*);
24static int DiagMatLogDeterminant(void*, double *);
25
26/* static int DiagMatScale(double *, diagmat *); */
27
28static 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
40static 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
49static 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
66static int DiagMatGetSize(void *AA, int *n){
67 diagmat* A=(diagmat*) AA;
68 *n=A->n;
69 return 0;
70}
71
72static 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
81static 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
88static int DiagMatAddElement(void*A, int nrow, double dd){
89 diagmat* AA = (diagmat*)A;
90 AA->val[nrow] += dd;
91 return 0;
92}
93
94static 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
101static 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
111static 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
123static 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
131static int DiagMatInvert(void *A){
132 return 0;
133}
134
135static 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
146static 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
158static 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}
168static 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
178static 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}
188static 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
198static int DiagMatFull(void*A, int* dfull){
199 *dfull=1;
200 return 0;
201}
202
203static struct DSDPDualMat_Ops sdmatopsp;
204static struct DSDPDualMat_Ops sdmatopsu;
205static const char* diagmatname="DIAGONAL";
206
207static 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}
227static 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"
250int 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"
271int 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
290static 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
301static 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}
316static 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
332static struct DSDPDSMat_Ops dsdiagmatopsp;
333static struct DSDPDSMat_Ops dsdiagmatopsu;
334
335#undef __FUNCT__
336#define __FUNCT__ "DSDPDiagDSMatP"
337int 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"
351int 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
364static 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
371static 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
382static 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
393static 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
404static 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
411static int DiagMatOnProcessor(void*A,int row,int*flag){
412 *flag=1;
413 return 0;
414}
415
416static int DiagAssemble(void*M){
417 return 0;
418}
419
420static struct DSDPSchurMat_Ops dsdpdiagschurops;
421
422#undef __FUNCT__
423#define __FUNCT__ "DSDPDiagSchurOps"
424static 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"
449int 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 DSDPDSMatOpsInitialize(struct DSDPDSMat_Ops *aops)
Set pointers to null.
Definition dsdpdsmat.c:214
Structure of function pointers that each SDP Delta S matrix type (sparse, dense, diagonal,...
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Structure of function pointers that each symmetric positive definite matrix type (dense,...
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide.
Error handling, printing, and profiling.
Symmetric Delta S matrix for one block in the semidefinite cone.
Table of function pointers that operate on the S matrix.