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,...
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
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,...
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