DSDP
dlpack.c
Go to the documentation of this file.
1#include "dsdpsys.h"
2#include "dsdpvec.h"
3#include "dsdplapack.h"
4
10typedef struct{
11 char UPLO;
12 double *val;
13 double *v2;
14 double *sscale;
15 int scaleit;
16 int n;
17 int owndata;
18} dtpumat;
19
20static const char* lapackname="DENSE,SYMMETRIC,PACKED STORAGE";
21
22int DTPUMatEigs(void*AA,double W[],double IIWORK[], int nn1, double *mineig){
23 dtpumat* AAA=(dtpumat*) AA;
24 ffinteger info,INFO=0,M,N=AAA->n;
25 ffinteger IL=1,IU=1,LDZ=1,IFAIL;
26 ffinteger *IWORK=(ffinteger*)IIWORK;
27 double *AP=AAA->val,ABSTOL=1e-13;
28 double Z=0,VL=-1e10,VU=1;
29 double *WORK;
30 char UPLO=AAA->UPLO,JOBZ='N',RANGE='I';
31
32 DSDPCALLOC2(&WORK,double,7*N,&info);DSDPCHKERR(info);
33 DSDPCALLOC2(&IWORK,ffinteger,5*N,&info);DSDPCHKERR(info);
34 dspevx(&JOBZ,&RANGE,&UPLO,&N,AP,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,IWORK,&IFAIL,&INFO);
35
36 /*
37 DSDPCALLOC2(&WORK,double,2*N,&info);
38 LWORK=2*N;
39 dspevd(&JOBZ,&UPLO,&N,AP,W,&Z,&LDZ,WORK,&LWORK,IWORK,&LIWORK,&INFO);
40 */
41 *mineig=W[0];
42 DSDPFREE(&WORK,&info);DSDPCHKERR(info);
43 DSDPFREE(&IWORK,&info);DSDPCHKERR(info);
44 return INFO;
45}
46
47
48#undef __FUNCT__
49#define __FUNCT__ "DSDPLAPACKROUTINE"
50static void dtpuscalevec(double alpha, double v1[], double v2[], double v3[], int n){
51 int i;
52 for (i=0;i<n;i++){
53 v3[i] = (alpha*v1[i]*v2[i]);
54 }
55}
56
57static void dtpuscalemat(double vv[], double ss[], int n){
58 int i;
59 for (i=0;i<n;i++,vv+=i){
60 dtpuscalevec(ss[i],vv,ss,vv,i+1);
61 }
62}
63
64static int DTPUMatCreateWData(int n,double nz[], int nnz, dtpumat**M){
65 int i,nn=(n*n+n)/2,info;
66 double dtmp;
67 dtpumat*M23;
68 if (nnz<nn){DSDPSETERR1(2,"Array must have length of : %d \n",nn);}
69 for (i=0;i<nnz;i++) dtmp=nz[i];
70 DSDPCALLOC1(&M23,dtpumat,&info);DSDPCHKERR(info);
71 DSDPCALLOC2(&M23->sscale,double,n,&info);DSDPCHKERR(info);
72 M23->owndata=0; M23->val=nz; M23->n=n; M23->UPLO='U';
73 for (i=0;i<n;i++)M23->sscale[i]=1.0;
74 M23->scaleit=0;
75 *M=M23;
76 return 0;
77}
78
79
80
81#undef __FUNCT__
82#define __FUNCT__ "DSDPLAPACK ROUTINE"
83
84
85static int DTPUMatMult(void* AA, double x[], double y[], int n){
86 dtpumat* A=(dtpumat*) AA;
87 ffinteger ione=1,N=n;
88 double BETA=0.0,ALPHA=1.0;
89 double *AP=A->val,*Y=y,*X=x;
90 char UPLO=A->UPLO;
91
92 if (A->n != n) return 1;
93 if (x==0 && n>0) return 3;
94 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione);
95 return 0;
96}
97
98static int DTPUMatSolve(void* AA, double b[], double x[], int n){
99 dtpumat* A=(dtpumat*) AA;
100 ffinteger INFO,NRHS=1,LDB=A->n,N=A->n;
101 double *AP=A->val,*ss=A->sscale;
102 char UPLO=A->UPLO;
103
104 if (N>0) LDB=N;
105 dtpuscalevec(1.0,ss,b,x,n);
106 dpptrs(&UPLO, &N, &NRHS, AP, x, &LDB, &INFO );
107 dtpuscalevec(1.0,ss,x,x,n);
108 return INFO;
109}
110
111static int DTPUMatCholeskyFactor(void* AA, int *flag){
112 dtpumat* A=(dtpumat*) AA;
113 int i;
114 ffinteger INFO,LDA=1,N=A->n;
115 double *AP=A->val,*ss=A->sscale,*v;
116 char UPLO=A->UPLO;
117
118 if (N<=0) LDA=1;
119 else LDA=N;
120 if (A->scaleit){
121 for (v=AP,i=0;i<N;i++){ ss[i]=*v;v+=(i+2);}
122 for (i=0;i<N;i++){ ss[i]=1.0/sqrt(fabs(ss[i])+1.0e-8); }
123 dtpuscalemat(AP,ss,N);
124 }
125 dpptrf(&UPLO, &N, AP, &INFO );
126 *flag=INFO;
127 return 0;
128}
129
130static int DTPUMatShiftDiagonal(void* AA, double shift){
131 dtpumat* A=(dtpumat*) AA;
132 int i,n=A->n;
133 double *v=A->val;
134 for (i=0; i<n; i++){
135 *v+=shift;
136 v+=i+2;
137 }
138 return 0;
139}
140
141
142#undef __FUNCT__
143#define __FUNCT__ "DTPUMatAssemble"
144static int DTPUMatAssemble(void*M){
145 int info;
146 double shift=1.0e-15;
147 DSDPFunctionBegin;
148 info= DTPUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
149 DSDPFunctionReturn(0);
150}
151
152static int DTPUMatRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
153 int i;
154 DSDPFunctionBegin;
155 *ncols = row+1;
156 for (i=0;i<=row;i++){
157 cols[i]=1.0;
158 }
159 for (i=row+1;i<nrows;i++){
160 cols[i]=0.0;
161 }
162 DSDPFunctionReturn(0);
163}
164
165
166#undef __FUNCT__
167#define __FUNCT__ "DTPUMatDiag"
168static int DTPUMatDiag(void*M, int row, double dd){
169 int ii;
170 dtpumat* ABA=(dtpumat*)M;
171 DSDPFunctionBegin;
172 ii=row*(row+1)/2+row;
173 ABA->val[ii] +=dd;
174 DSDPFunctionReturn(0);
175}
176#undef __FUNCT__
177#define __FUNCT__ "DTPUMatDiag2"
178static int DTPUMatDiag2(void*M, double diag[], int m){
179 int row,ii;
180 dtpumat* ABA=(dtpumat*)M;
181 DSDPFunctionBegin;
182 for (row=0;row<m;row++){
183 ii=row*(row+1)/2+row;
184 ABA->val[ii] +=diag[row];
185 }
186 DSDPFunctionReturn(0);
187}
188
189static int DTPUMatAddRow(void* AA, int nrow, double dd, double row[], int n){
190 dtpumat* A=(dtpumat*) AA;
191 ffinteger ione=1,nn,nnn;
192 double *vv=A->val;
193
194 nnn=nrow*(nrow+1)/2;
195 nn=nrow+1;
196 daxpy(&nn,&dd,row,&ione,vv+nnn,&ione);
197
198 return 0;
199}
200
201static int DTPUMatZero(void* AA){
202 dtpumat* A=(dtpumat*) AA;
203 int mn=A->n*(A->n+1)/2;
204 double *vv=A->val;
205 memset((void*)vv,0,mn*sizeof(double));
206 return 0;
207}
208
209static int DTPUMatGetSize(void *AA, int *n){
210 dtpumat* A=(dtpumat*) AA;
211 *n=A->n;
212 return 0;
213}
214
215static int DTPUMatDestroy(void* AA){
216 int info;
217 dtpumat* A=(dtpumat*) AA;
218 if (A && A->owndata){DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
219 if (A && A->sscale) {DSDPFREE(&A->sscale,&info);DSDPCHKERR(info);}
220 if (A) {DSDPFREE(&A,&info);DSDPCHKERR(info);}
221 return 0;
222}
223
224static int DTPUMatView(void* AA){
225 dtpumat* M=(dtpumat*) AA;
226 int i,j,kk=0;
227 double *val=M->val;
228 for (i=0; i<M->n; i++){
229 for (j=0; j<=i; j++){
230 printf(" %9.2e",val[kk]);
231 kk++;
232 }
233 printf("\n");
234 }
235 return 0;
236}
237
238
239#include "dsdpschurmat_impl.h"
240#include "dsdpsys.h"
241static struct DSDPSchurMat_Ops dsdpmmatops;
242
243static int DSDPInitSchurOps(struct DSDPSchurMat_Ops* mops){
244 int info;
245 DSDPFunctionBegin;
246 info=DSDPSchurMatOpsInitialize(mops);DSDPCHKERR(info);
247 mops->matrownonzeros=DTPUMatRowNonzeros;
248 mops->matscaledmultiply=DTPUMatMult;
249 mops->mataddrow=DTPUMatAddRow;
250 mops->mataddelement=DTPUMatDiag;
251 mops->matadddiagonal=DTPUMatDiag2;
252 mops->matshiftdiagonal=DTPUMatShiftDiagonal;
253 mops->matassemble=DTPUMatAssemble;
254 mops->matfactor=DTPUMatCholeskyFactor;
255 mops->matsolve=DTPUMatSolve;
256 mops->matdestroy=DTPUMatDestroy;
257 mops->matzero=DTPUMatZero;
258 mops->matview=DTPUMatView;
259 mops->id=1;
260 mops->matname=lapackname;
261 DSDPFunctionReturn(0);
262}
263
264#undef __FUNCT__
265#define __FUNCT__ "DSDPGetLAPACKPUSchurOps"
266int DSDPGetLAPACKPUSchurOps(int n,struct DSDPSchurMat_Ops** sops, void** mdata){
267 int info,nn=n*(n+1)/2;
268 double *vv;
269 dtpumat*AA;
270 DSDPFunctionBegin;
271 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
272 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
273 AA->owndata=1;
274 AA->scaleit=1;
275 info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info);
276 *sops=&dsdpmmatops;
277 *mdata=(void*)AA;
278 DSDPFunctionReturn(0);
279}
280
281
282static void daddrow(double *v, double alpha, int i, double row[], int n){
283 double *s1;
284 ffinteger j,nn=n,ione=1;
285 nn=i+1; s1=v+i*(i+1)/2;
286 daxpy(&nn,&alpha,s1,&ione,row,&ione);
287 for (j=i+1;j<n;j++){
288 s1+=j;
289 row[j]+=alpha*s1[i];
290 }
291 return;
292}
293
294static int DTPUMatInverseMult(void* AA, int indx[], int nind, double x[], double y[], int n){
295 dtpumat* A=(dtpumat*) AA;
296 ffinteger ione=1,N=n;
297 double BETA=0.0,ALPHA=1.0;
298 double *AP=A->v2,*Y=y,*X=x;
299 int i,ii;
300 char UPLO=A->UPLO;
301
302 if (A->n != n) return 1;
303 if (x==0 && n>0) return 3;
304
305 if (nind<n/4 ){
306 memset((void*)y,0,n*sizeof(double));
307 for (ii=0;ii<nind;ii++){
308 i=indx[ii]; ALPHA=x[i];
309 daddrow(AP,ALPHA,i,y,n);
310 }
311 } else {
312 ALPHA=1.0;
313 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione);
314 }
315 return 0;
316}
317
318
319static int DTPUMatCholeskyBackward(void* AA, double b[], double x[], int n){
320 dtpumat* M=(dtpumat*) AA;
321 ffinteger N=M->n,INCX=1;
322 double *AP=M->val,*ss=M->sscale;
323 char UPLO=M->UPLO,TRANS='N',DIAG='N';
324 memcpy(x,b,N*sizeof(double));
325 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX );
326 dtpuscalevec(1.0,ss,x,x,n);
327 return 0;
328}
329
330
331static int DTPUMatCholeskyForward(void* AA, double b[], double x[], int n){
332 dtpumat* M=(dtpumat*) AA;
333 ffinteger N=M->n,INCX=1;
334 double *AP=M->val,*ss=M->sscale;
335 char UPLO=M->UPLO,TRANS='T',DIAG='N';
336 dtpuscalevec(1.0,ss,b,x,n);
337 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX);
338 return 0;
339}
340
341static int DenseSymPSDCholeskyForwardMultiply(void* AA, double x[], double y[], int n){
342 dtpumat* A=(dtpumat*) AA;
343 int i,j,k=0;
344 ffinteger N=A->n;
345 double *AP=A->val,*ss=A->sscale;
346
347 if (x==0 && N>0) return 3;
348 for (i=0;i<N;i++){
349 for (j=0;j<=i;j++){
350 y[i]+=AP[k]*x[j];
351 k++;
352 }
353 }
354 for (i=0;i<N;i++){ y[i]=y[i]/ss[i];}
355 return 0;
356}
357
358static int DTPUMatLogDet(void* AA, double *dd){
359 dtpumat* A=(dtpumat*) AA;
360 int i,n=A->n;
361 double d=0,*v=A->val,*ss=A->sscale;
362 for (i=0; i<n; i++){
363 if (*v<=0) return 1;
364 d+=2*log(*v/ss[i]);
365 v+=i+2;
366 }
367 *dd=d;
368 return 0;
369}
370
371static int DTPUMatInvert(void* AA){
372 dtpumat* A=(dtpumat*) AA;
373 ffinteger INFO,N=A->n,nn=N*(N+1)/2;
374 double *v=A->val,*AP=A->v2,*ss=A->sscale;
375 char UPLO=A->UPLO;
376 memcpy((void*)AP,(void*)v,nn*sizeof(double));
377 dpptri(&UPLO, &N, AP, &INFO );
378 if (INFO){
379 INFO=DTPUMatShiftDiagonal(AA,1e-11);
380 memcpy((void*)AP,(void*)v,nn*sizeof(double));
381 dpptri(&UPLO, &N, AP, &INFO );
382 }
383 if (A->scaleit){
384 dtpuscalemat(AP,ss,N);
385 }
386 return INFO;
387}
388
389static int DTPUMatInverseAdd(void* AA, double alpha, double y[], int nn, int n){
390 dtpumat* A=(dtpumat*) AA;
391 ffinteger N=nn,ione=1;
392 double *v2=A->v2;
393 daxpy(&N,&alpha,v2,&ione,y,&ione);
394 return 0;
395}
396
397
398static int DTPUMatScaleDiagonal(void* AA, double dd){
399 dtpumat* A=(dtpumat*) AA;
400 int i,n=A->n;
401 double *v=A->val;
402 for (i=0; i<n; i++){
403 *v*=dd;
404 v+=i+2;
405 }
406 return 0;
407}
408
409static int DTPUMatOuterProduct(void* AA, double alpha, double x[], int n){
410 dtpumat* A=(dtpumat*) AA;
411 ffinteger ione=1,N=n;
412 double *v=A->val;
413 char UPLO=A->UPLO;
414 dspr(&UPLO,&N,&alpha,x,&ione,v);
415 return 0;
416}
417
418
419static int DenseSymPSDNormF2(void* AA, int n, double *dddot){
420 dtpumat* A=(dtpumat*) AA;
421 ffinteger ione=1,nn=A->n*(A->n+1)/2;
422 double dd,tt=sqrt(0.5),*val=A->val;
423 int info;
424 info=DTPUMatScaleDiagonal(AA,tt);
425 dd=dnrm2(&nn,val,&ione);
426 info=DTPUMatScaleDiagonal(AA,1.0/tt);
427 *dddot=dd*dd*2;
428 return 0;
429}
430
431
432/*
433static int DTPUMatFNorm2(void* AA, double *mnorm){
434 dtpumat* A=(dtpumat*) AA;
435 ffinteger ione=1,n;
436 double vv=0,*AP=A->val;
437 n=A->n*(A->n+1)/2;
438 vv=dnrm2(&n,AP,&ione);
439 *mnorm=2.0*vv;
440 return 1;
441}
442*/
443
444#include "dsdpdualmat_impl.h"
445#include "dsdpdatamat_impl.h"
446#include "dsdpxmat_impl.h"
447#include "dsdpdsmat_impl.h"
448
449
450
451static int DTPUMatFull(void*A, int*full){
452 *full=1;
453 return 0;
454}
455
456
457static int DTPUMatGetDenseArray(void* A, double *v[], int*n){
458 dtpumat* ABA=(dtpumat*)A;
459 *v=ABA->val;
460 *n=(ABA->n)*(ABA->n+1)/2;
461 return 0;
462}
463
464static int DTPUMatRestoreDenseArray(void* A, double *v[], int *n){
465 *v=0;*n=0;
466 return 0;
467}
468
469static int DDenseSetXMat(void*A, double v[], int nn, int n){
470 double *vv;
471 dtpumat* ABA=(dtpumat*)A;
472 vv=ABA->val;
473 if (v!=vv){
474 memcpy((void*)vv,(void*)v,nn*sizeof(double));
475 }
476 return 0;
477}
478
479static int DDenseVecVec(void* A, double x[], int n, double *v){
480 dtpumat* ABA=(dtpumat*)A;
481 int i,j,k=0;
482 double dd=0,*val=ABA->val;
483 *v=0.0;
484 for (i=0; i<n; i++){
485 for (j=0;j<i;j++){
486 dd+=2*x[i]*x[j]*val[k];
487 k++;
488 }
489 dd+=x[i]*x[i]*val[k];
490 k++;
491 }
492 *v=dd;
493 return 0;
494}
495
496static struct DSDPDSMat_Ops tdsdensematops;
497static int DSDPDSDenseInitializeOps(struct DSDPDSMat_Ops* densematops){
498 int info;
499 if (!densematops) return 0;
500 info=DSDPDSMatOpsInitialize(densematops); DSDPCHKERR(info);
501 densematops->matseturmat=DDenseSetXMat;
502 densematops->matview=DTPUMatView;
503 densematops->matdestroy=DTPUMatDestroy;
504 densematops->matgetsize=DTPUMatGetSize;
505 densematops->matzeroentries=DTPUMatZero;
506 densematops->matmult=DTPUMatMult;
507 densematops->matvecvec=DDenseVecVec;
508 densematops->id=1;
509 densematops->matname=lapackname;
510 return 0;
511}
512
513#undef __FUNCT__
514#define __FUNCT__ "DSDPCreateDSMatWithArray"
515int DSDPCreateDSMatWithArray(int n,double vv[], int nnz, struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
516 int info;
517 dtpumat*AA;
518 DSDPFunctionBegin;
519 info=DTPUMatCreateWData(n,vv,nnz,&AA); DSDPCHKERR(info);
520 AA->owndata=0;
521 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
522 *dsmatops=&tdsdensematops;
523 *dsmat=(void*)AA;
524 DSDPFunctionReturn(0);
525}
526
527
528#undef __FUNCT__
529#define __FUNCT__ "DSDPCreateDSMat"
530int DSDPCreateDSMat(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
531 int info,nn=n*(n+1)/2;
532 double *vv;
533 dtpumat*AA;
534 DSDPFunctionBegin;
535 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
536 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
537 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
538 *dsmatops=&tdsdensematops;
539 *dsmat=(void*)AA;
540 AA->owndata=1;
541 DSDPFunctionReturn(0);
542}
543
544static struct DSDPVMat_Ops turdensematops;
545
546static int DSDPDenseXInitializeOps(struct DSDPVMat_Ops* densematops){
547 int info;
548 if (!densematops) return 0;
549 info=DSDPVMatOpsInitialize(densematops); DSDPCHKERR(info);
550 densematops->matview=DTPUMatView;
551 densematops->matscalediagonal=DTPUMatScaleDiagonal;
552 densematops->matshiftdiagonal=DTPUMatShiftDiagonal;
553 densematops->mataddouterproduct=DTPUMatOuterProduct;
554 densematops->matdestroy=DTPUMatDestroy;
555 densematops->matfnorm2=DenseSymPSDNormF2;
556 densematops->matgetsize=DTPUMatGetSize;
557 densematops->matzeroentries=DTPUMatZero;
558 densematops->matgeturarray=DTPUMatGetDenseArray;
559 densematops->matrestoreurarray=DTPUMatRestoreDenseArray;
560 densematops->matmineig=DTPUMatEigs;
561 densematops->matmult=DTPUMatMult;
562 densematops->id=1;
563 densematops->matname=lapackname;
564 return 0;
565}
566
567#undef __FUNCT__
568#define __FUNCT__ "DSDPXMatCreate"
569int DSDPXMatPCreate(int n,struct DSDPVMat_Ops* *xops, void * *xmat){
570 int info,nn=n*(n+1)/2;
571 double *vv;
572 dtpumat*AA;
573 DSDPFunctionBegin;
574 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
575 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
576 AA->owndata=1;
577 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
578 *xops=&turdensematops;
579 *xmat=(void*)AA;
580 DSDPFunctionReturn(0);
581}
582
583#undef __FUNCT__
584#define __FUNCT__ "DSDPXMatCreate"
585int DSDPXMatPCreateWithData(int n,double nz[],int nnz,struct DSDPVMat_Ops* *xops, void * *xmat){
586 int i,info;
587 double dtmp;
588 dtpumat*AA;
589 DSDPFunctionBegin;
590 for (i=0;i<n*(n+1)/2;i++) dtmp=nz[i];
591 info=DTPUMatCreateWData(n,nz,nnz,&AA); DSDPCHKERR(info);
592 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
593 *xops=&turdensematops;
594 *xmat=(void*)AA;
595 DSDPFunctionReturn(0);
596}
597
598
599static struct DSDPDualMat_Ops sdmatops;
600static int SDualOpsInitialize(struct DSDPDualMat_Ops* sops){
601 int info;
602 if (sops==NULL) return 0;
603 info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info);
604 sops->matseturmat=DDenseSetXMat;
605 sops->matcholesky=DTPUMatCholeskyFactor;
606 sops->matsolveforward=DTPUMatCholeskyForward;
607 sops->matsolvebackward=DTPUMatCholeskyBackward;
608 sops->matinvert=DTPUMatInvert;
609 sops->matinverseadd=DTPUMatInverseAdd;
610 sops->matinversemultiply=DTPUMatInverseMult;
611 sops->matforwardmultiply=DenseSymPSDCholeskyForwardMultiply;
612 sops->matfull=DTPUMatFull;
613 sops->matdestroy=DTPUMatDestroy;
614 sops->matgetsize=DTPUMatGetSize;
615 sops->matview=DTPUMatView;
616 sops->matlogdet=DTPUMatLogDet;
617 sops->matname=lapackname;
618 sops->id=1;
619 return 0;
620}
621
622
623#undef __FUNCT__
624#define __FUNCT__ "DSDPLAPACKDualMatCreate"
625int DSDPLAPACKPUDualMatCreate(int n,struct DSDPDualMat_Ops **sops, void**smat){
626 int info,nn=n*(n+1)/2;
627 double *vv;
628 dtpumat*AA;
629 DSDPFunctionBegin;
630 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
631 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
632 AA->owndata=1;;
633 AA->scaleit=1;
634 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
635 *sops=&sdmatops;
636 *smat=(void*)AA;
637 DSDPFunctionReturn(0);
638}
639
640static int switchptr(void* SD,void *SP){
641 dtpumat *s1,*s2;
642 s1=(dtpumat*)(SD);
643 s2=(dtpumat*)(SP);
644 s1->v2=s2->val;
645 s2->v2=s1->val;
646 return 0;
647}
648
649
650#undef __FUNCT__
651#define __FUNCT__ "DSDPLAPACKDualMatCreate2"
652int DSDPLAPACKPUDualMatCreate2(int n,
653 struct DSDPDualMat_Ops **sops1, void**smat1,
654 struct DSDPDualMat_Ops **sops2, void**smat2){
655 int info;
656 DSDPFunctionBegin;
657 info=DSDPLAPACKPUDualMatCreate(n,sops1,smat1);DSDPCHKERR(info);
658 info=DSDPLAPACKPUDualMatCreate(n,sops2,smat2);DSDPCHKERR(info);
659 info=switchptr(*smat1,*smat2);
660 DSDPFunctionReturn(0);
661}
662
663
664typedef struct {
665 int neigs;
666 double *eigval;
667 double *an;
668} Eigen;
669
670typedef struct {
671 dtpumat* AA;
672 double alpha;
673 Eigen Eig;
674} dvechmat;
675
676#undef __FUNCT__
677#define __FUNCT__ "CreateDvechmatWData"
678static int CreateDvechmatWdata(int n, double alpha, double vv[], dvechmat **A){
679 int info,nn=(n*n+n)/2;
680 dvechmat* V;
681 DSDPCALLOC1(&V,dvechmat,&info);DSDPCHKERR(info);
682 info=DTPUMatCreateWData(n,vv,nn,&V->AA); DSDPCHKERR(info);
683 V->Eig.neigs=-1;
684 V->Eig.eigval=0;
685 V->Eig.an=0;
686 V->alpha=alpha;
687 *A=V;
688 return 0;
689}
690
691
692static int DvechmatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int n){
693 int k;
694 *nnzz=n;
695 for (k=0;k<n;k++) nz[k]++;
696 return 0;
697}
698
699static int DTPUMatGetRowAdd(void* AA, int nrow, double ytmp, double row[], int n){
700 dtpumat* A=(dtpumat*) AA;
701 ffinteger i,k,nnn=n;
702 double *v=A->val;
703
704 nnn=nrow*(nrow+1)/2;
705 for (i=0;i<nrow;i++){
706 row[i]+=ytmp*v[nnn+i];
707 }
708 row[nrow]+=ytmp*v[nnn+nrow];
709 for (i=nrow+1;i<n;i++){
710 k=i*(i+1)/2+nrow;
711 row[i]+=ytmp*v[k];
712 }
713 return 0;
714}
715
716static int DvechmatGetRowAdd(void* AA, int trow, double scl, double r[], int m){
717 int info;
718 dvechmat* A=(dvechmat*)AA;
719 info=DTPUMatGetRowAdd((void*)A->AA ,trow,scl*A->alpha,r,m);
720 return 0;
721}
722
723static int DvechmatAddMultiple(void* AA, double alpha, double r[], int nnn, int n){
724 dvechmat* A=(dvechmat*)AA;
725 ffinteger nn=nnn, ione=1;
726 double *val=A->AA->val;
727 alpha*=A->alpha;
728 daxpy(&nn,&alpha,val,&ione,r,&ione);
729 return 0;
730}
731
732
733static int DvechEigVecVec(void*, double[], int, double*);
734static int DvechmatVecVec(void* AA, double x[], int n, double *v){
735 dvechmat* A=(dvechmat*)AA;
736 int i,j,k=0;
737 double dd=0,*val=A->AA->val;
738 *v=0.0;
739 if (A->Eig.neigs<n/5){
740 i=DvechEigVecVec(AA,x,n,&dd);
741 *v=dd*A->alpha;
742 } else {
743 for (i=0; i<n; i++){
744 for (j=0;j<i;j++){
745 dd+=2*x[i]*x[j]*val[k];
746 k++;
747 }
748 dd+=x[i]*x[i]*val[k];
749 k++;
750 }
751 *v=dd*A->alpha;
752 }
753 return 0;
754}
755
756
757static int DvechmatFNorm2(void* AA, int n, double *v){
758 dvechmat* A=(dvechmat*)AA;
759 long int i,j,k=0;
760 double dd=0,*x=A->AA->val;
761 for (i=0; i<n; i++){
762 for (j=0;j<i;j++){
763 dd+=2*x[k]*x[k];
764 k++;
765 }
766 dd+=x[k]*x[k];
767 k++;
768 }
769 *v=dd*A->alpha*A->alpha;
770 return 0;
771}
772
773
774static int DvechmatCountNonzeros(void* AA, int *nnz, int n){
775 *nnz=n*(n+1)/2;
776 return 0;
777}
778
779
780static int DvechmatDot(void* AA, double x[], int nn, int n, double *v){
781 dvechmat* A=(dvechmat*)AA;
782 ffinteger ione=1,nnn=nn;
783 double dd,*val=A->AA->val;
784 dd=ddot(&nnn,val,&ione,x,&ione);
785 *v=2*dd*A->alpha;
786 return 0;
787}
788
789/*
790static int DvechmatNormF2(void* AA, int n, double *v){
791 dvechmat* A=(dvechmat*)AA;
792 return(DTPUMatNormF2((void*)(A->AA), n,v));
793}
794*/
795#undef __FUNCT__
796#define __FUNCT__ "DvechmatDestroy"
797static int DvechmatDestroy(void* AA){
798 dvechmat* A=(dvechmat*)AA;
799 int info;
800 info=DTPUMatDestroy((void*)(A->AA));
801 if (A->Eig.an){DSDPFREE(&A->Eig.an,&info);DSDPCHKERR(info);}
802 if (A->Eig.eigval){DSDPFREE(&A->Eig.eigval,&info);DSDPCHKERR(info);}
803 A->Eig.neigs=-1;
804 DSDPFREE(&A,&info);DSDPCHKERR(info);
805 return 0;
806}
807
808
809static int DvechmatView(void* AA){
810 dvechmat* A=(dvechmat*)AA;
811 dtpumat* M=A->AA;
812 int i,j,kk=0;
813 double *val=M->val;
814 for (i=0; i<M->n; i++){
815 for (j=0; j<=i; j++){
816 printf(" %4.2e",A->alpha*val[kk]);
817 kk++;
818 }
819 printf(" \n");
820 }
821 return 0;
822}
823
824
825#undef __FUNCT__
826#define __FUNCT__ "DSDPCreateDvechmatEigs"
827static int CreateEigenLocker(Eigen *E,int neigs, int n){
828 int info;
829 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info);
830 DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info);
831 E->neigs=neigs;
832 return 0;
833}
834
835
836static int EigMatSetEig(Eigen* A,int row, double eigv, double v[], int n){
837 double *an=A->an;
838 A->eigval[row]=eigv;
839 memcpy((void*)(an+n*row),(void*)v,n*sizeof(double));
840 return 0;
841}
842
843
844static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n){
845 double* an=A->an;
846 *eigenvalue=A->eigval[row];
847 memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
848 return 0;
849}
850
851
852static int DvechmatComputeEigs(dvechmat*,double[],int,double[],int,double[],int,int[],int);
853
854static int DvechmatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){
855
856 int info;
857 dvechmat* A=(dvechmat*)AA;
858 if (A->Eig.neigs>=0) return 0;
859 info=DvechmatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2);DSDPCHKERR(info);
860 return 0;
861}
862
863static int DvechmatGetRank(void *AA,int *rank, int n){
864 dvechmat* A=(dvechmat*)AA;
865 if (A->Eig.neigs>=0){
866 *rank=A->Eig.neigs;
867 } else {
868 DSDPSETERR(1,"Vech Matrix not factored yet\n");
869 }
870 return 0;
871}
872
873static int DvechmatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indz[], int *nind){
874 dvechmat* A=(dvechmat*)AA;
875 int i,info;
876 double dd;
877 if (A->Eig.neigs>0){
878 info=EigMatGetEig(&A->Eig,rank,&dd,vv,n);DSDPCHKERR(info);
879 *nind=n;
880 *eigenvalue=dd*A->alpha;
881 for (i=0;i<n;i++){ indz[i]=i;}
882 } else {
883 DSDPSETERR(1,"Vech Matrix not factored yet\n");
884 }
885 return 0;
886}
887
888static int DvechEigVecVec(void* AA, double v[], int n, double *vv){
889 dvechmat* A=(dvechmat*)AA;
890 int i,rank,neigs;
891 double *an,dd,ddd=0,*eigval;
892 if (A->Eig.neigs>=0){
893 an=A->Eig.an;
894 neigs=A->Eig.neigs;
895 eigval=A->Eig.eigval;
896 for (rank=0;rank<neigs;rank++){
897 for (dd=0,i=0;i<n;i++){
898 dd+=v[i]*an[i];
899 }
900 an+=n;
901 ddd+=dd*dd*eigval[rank];
902 }
903 *vv=ddd*A->alpha;
904 } else {
905 DSDPSETERR(1,"Vech Matrix not factored yet\n");
906 }
907 return 0;
908}
909
910
911static struct DSDPDataMat_Ops dvechmatops;
912static const char *datamatname="DENSE VECH MATRIX";
913
914static int DvechmatOpsInitialize(struct DSDPDataMat_Ops *sops){
915 int info;
916 if (sops==NULL) return 0;
917 info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
918 sops->matvecvec=DvechmatVecVec;
919 sops->matdot=DvechmatDot;
920 sops->mataddrowmultiple=DvechmatGetRowAdd;
921 sops->mataddallmultiple=DvechmatAddMultiple;
922 sops->matview=DvechmatView;
923 sops->matdestroy=DvechmatDestroy;
924 sops->matfactor2=DvechmatFactor;
925 sops->matgetrank=DvechmatGetRank;
926 sops->matgeteig=DvechmatGetEig;
927 sops->matrownz=DvechmatGetRowNnz;
928 sops->matfnorm2=DvechmatFNorm2;
929 sops->matnnz=DvechmatCountNonzeros;
930 sops->id=1;
931 sops->matname=datamatname;
932 return 0;
933}
934
935#undef __FUNCT__
936#define __FUNCT__ "DSDPGetDmat"
937int DSDPGetDMat(int n,double alpha, double *val, struct DSDPDataMat_Ops**sops, void**smat){
938 int info,k;
939 double dtmp;
940 dvechmat* A;
941 DSDPFunctionBegin;
942
943 for (k=0;k<n*(n+1)/2;++k) dtmp=val[k];
944 info=CreateDvechmatWdata(n,alpha,val,&A); DSDPCHKERR(info);
945 A->Eig.neigs=-1;
946 info=DvechmatOpsInitialize(&dvechmatops); DSDPCHKERR(info);
947 if (sops){*sops=&dvechmatops;}
948 if (smat){*smat=(void*)A;}
949 DSDPFunctionReturn(0);
950}
951
952
953#undef __FUNCT__
954#define __FUNCT__ "DvechmatComputeEigs"
955static int DvechmatComputeEigs(dvechmat* AA,double DD[], int nn0, double W[], int n, double WORK[], int n1, int iiptr[], int n2){
956
957 int i,j,k,neigs,info;
958 long int *i2darray=(long int*)DD;
959 int ownarray1=0,ownarray2=0,ownarray3=0;
960 double *val=AA->AA->val;
961 double *dmatarray=0,*dworkarray=0,eps=1.0e-12;
962 int nn1=0,nn2=0;
963
964 /* create a dense array in which to put numbers */
965 if (n*n>nn1){
966 DSDPCALLOC2(&dmatarray,double,(n*n),&info); DSDPCHKERR(info);
967 ownarray1=1;
968 }
969 memset((void*)dmatarray,0,n*n*sizeof(double));
970
971 if (n*n>nn2){
972 DSDPCALLOC2(&dworkarray,double,(n*n),&info); DSDPCHKERR(info);
973 ownarray2=1;
974 }
975
976 if (n*n*sizeof(long int)>nn0*sizeof(double)){
977 DSDPCALLOC2(&i2darray,long int,(n*n),&info); DSDPCHKERR(info);
978 ownarray3=1;
979 }
980
981
982 for (k=0,i=0; i<n; i++){
983 for (j=0; j<=i; j++){
984 dmatarray[i*n+j] += val[k];
985 if (i!=j){
986 dmatarray[j*n+i] += val[k];
987 }
988 k++;
989 }
990 }
991 /* Call LAPACK to compute the eigenvalues */
992 info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n,
993 W,n,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
994
995 /* Count the nonzero eigenvalues */
996 for (neigs=0,i=0;i<n;i++){
997 if (fabs(W[i])> eps ){ neigs++;}
998 }
999
1000 info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info);
1001
1002 /* Copy into structure */
1003 for (neigs=0,i=0; i<n; i++){
1004 if (fabs(W[i]) > eps){
1005 info=EigMatSetEig(&AA->Eig,neigs,W[i],dmatarray+n*i,n);DSDPCHKERR(info);
1006 neigs++;
1007 }
1008 }
1009
1010 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
1011 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
1012 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
1013 return 0;
1014}
1015
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition dsdpdatamat.c:47
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
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,...
DSDP uses BLAS and LAPACK for many of its operations.
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.
Vector operations used by the solver.
int DSDPVMatOpsInitialize(struct DSDPVMat_Ops *aops)
Set function pointers to null.
Definition dsdpxmat.c:377
Structure of function pointers that each dense matrix array type (upper full, packed symmetric,...
Symmetric Delta S matrix for one block in the semidefinite cone.
Table of function pointers that operate on the data matrix.
Table of function pointers that operate on the S matrix.
Table of function pointers that operate on the dense matrix.