1
2 #include "dsdpschurmat_impl.h"
3 #include "dsdpdualmat_impl.h"
4 #include "dsdpdsmat_impl.h"
5 #include "dsdpsys.h"
6 /*!
7 \file diag.c
8 \brief DSDPDualMat, DSDPDSMat, and the DSDPSchurMat implentations for
9 diagonal matrices.
10 */
11
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
DiagMatCreate(int n,diagmat ** M)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
DiagMatDestroy(void * AA)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
DiagMatMult(void * AA,double x[],double y[],int n)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
DiagMatGetSize(void * AA,int * n)66 static int DiagMatGetSize(void *AA, int *n){
67 diagmat* A=(diagmat*) AA;
68 *n=A->n;
69 return 0;
70 }
71
DiagMatView(void * AA)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
DiagMatAddRow2(void * AA,int nrow,double dd,double row[],int n)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
DiagMatAddElement(void * A,int nrow,double dd)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
DiagMatZeros(void * A)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
DiagMatSolve(void * A,double b[],double x[],int n)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
DiagMatSolve2(void * A,int indx[],int nindx,double b[],double x[],int n)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
DiagMatCholeskySolveBackward(void * A,double b[],double x[],int n)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
DiagMatInvert(void * A)131 static int DiagMatInvert(void *A){
132 return 0;
133 }
134
DiagMatCholeskyFactor(void * A,int * flag)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
DiagMatLogDeterminant(void * A,double * dd)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
DiagMatTakeUREntriesP(void * A,double dd[],int nn,int n)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 }
DiagMatTakeUREntriesU(void * A,double dd[],int nn,int n)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
DiagMatInverseAddP(void * A,double alpha,double dd[],int nn,int n)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 }
DiagMatInverseAddU(void * A,double alpha,double dd[],int nn,int n)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
DiagMatFull(void * A,int * dfull)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
DiagDualOpsInitializeP(struct DSDPDualMat_Ops * sops)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 }
DiagDualOpsInitializeU(struct DSDPDualMat_Ops * sops)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"
DSDPDiagDualMatCreateP(int n,struct DSDPDualMat_Ops ** sops1,void ** smat1,struct DSDPDualMat_Ops ** sops2,void ** smat2)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"
DSDPDiagDualMatCreateU(int n,struct DSDPDualMat_Ops ** sops1,void ** smat1,struct DSDPDualMat_Ops ** sops2,void ** smat2)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
DiagMatVecVec(void * A,double x[],int n,double * vv)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
DDiagDSMatOpsInitP(struct DSDPDSMat_Ops * ddiagops)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 }
DDiagDSMatOpsInitU(struct DSDPDSMat_Ops * ddiagops)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"
DSDPCreateDiagDSMatP(int n,struct DSDPDSMat_Ops ** dsmatops,void ** dsmat)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"
DSDPCreateDiagDSMatU(int n,struct DSDPDSMat_Ops ** dsmatops,void ** dsmat)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
DiagRowNonzeros(void * M,int row,double cols[],int * ncols,int nrows)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
DiagAddDiag(void * M,double diag[],int m)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
DiagMultiply(void * M,double xin[],double xout[],int m)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
DiagShiftDiag(void * M,double dd)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
DiagAddElement(void * M,int ii,double dd)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
DiagMatOnProcessor(void * A,int row,int * flag)411 static int DiagMatOnProcessor(void*A,int row,int*flag){
412 *flag=1;
413 return 0;
414 }
415
DiagAssemble(void * M)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"
DiagSchurOps(struct DSDPSchurMat_Ops * sops)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"
DSDPGetDiagSchurMat(int n,struct DSDPSchurMat_Ops ** sops,void ** data)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 }
459