DSDP
dlpack.c
Go to the documentation of this file.
1 #include "dsdpsys.h"
2 #include "dsdpvec.h"
3 #include "dsdplapack.h"
4 
10 typedef 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 
20 static const char* lapackname="DENSE,SYMMETRIC,PACKED STORAGE";
21 
22 int 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"
50 static 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 
57 static 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 
64 static 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 
85 static 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 
98 static 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 
111 static 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 
130 static 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"
144 static 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 
152 static 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"
168 static 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"
178 static 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 
189 static 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 
201 static 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 
209 static int DTPUMatGetSize(void *AA, int *n){
210  dtpumat* A=(dtpumat*) AA;
211  *n=A->n;
212  return 0;
213 }
214 
215 static 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 
224 static 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"
241 static struct DSDPSchurMat_Ops dsdpmmatops;
242 
243 static 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"
266 int 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 
282 static 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 
294 static 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 
319 static 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 
331 static 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 
341 static 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 
358 static 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 
371 static 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 
389 static 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 
398 static 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 
409 static 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 
419 static 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 /*
433 static 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 
451 static int DTPUMatFull(void*A, int*full){
452  *full=1;
453  return 0;
454 }
455 
456 
457 static 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 
464 static int DTPUMatRestoreDenseArray(void* A, double *v[], int *n){
465  *v=0;*n=0;
466  return 0;
467 }
468 
469 static 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 
479 static 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 
496 static struct DSDPDSMat_Ops tdsdensematops;
497 static 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"
515 int 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"
530 int 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 
544 static struct DSDPVMat_Ops turdensematops;
545 
546 static 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"
569 int 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"
585 int 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 
599 static struct DSDPDualMat_Ops sdmatops;
600 static 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"
625 int 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 
640 static 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"
652 int 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 
664 typedef struct {
665  int neigs;
666  double *eigval;
667  double *an;
668 } Eigen;
669 
670 typedef struct {
671  dtpumat* AA;
672  double alpha;
673  Eigen Eig;
674 } dvechmat;
675 
676 #undef __FUNCT__
677 #define __FUNCT__ "CreateDvechmatWData"
678 static 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 
692 static 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 
699 static 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 
716 static 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 
723 static 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 
733 static int DvechEigVecVec(void*, double[], int, double*);
734 static 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 
757 static 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 
774 static int DvechmatCountNonzeros(void* AA, int *nnz, int n){
775  *nnz=n*(n+1)/2;
776  return 0;
777 }
778 
779 
780 static 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 /*
790 static 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"
797 static 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 
809 static 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"
827 static 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 
836 static 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 
844 static 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 
852 static int DvechmatComputeEigs(dvechmat*,double[],int,double[],int,double[],int,int[],int);
853 
854 static 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 
863 static 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 
873 static 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 
888 static 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 
911 static struct DSDPDataMat_Ops dvechmatops;
912 static const char *datamatname="DENSE VECH MATRIX";
913 
914 static 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"
937 int 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"
955 static 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 DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Definition: dsdpdualmat.c:423
Structure of function pointers that each dense matrix array type (upper full, packed symmetric...
Table of function pointers that operate on the dense matrix.
Definition: dsdpxmat_impl.h:13
Error handling, printing, and profiling.
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition: dsdpdatamat.c:47
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.
Structure of function pointers that each SDP data matrix type (sparse, dense, constant, identity, ...) 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
Vector operations used by the solver.
Table of function pointers that operate on the data matrix.
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.
DSDP uses BLAS and LAPACK for many of its operations.
Symmetric Delta S matrix for one block in the semidefinite cone.
int DSDPVMatOpsInitialize(struct DSDPVMat_Ops *aops)
Set function pointers to null.
Definition: dsdpxmat.c:377