DSDP
dufull.c
Go to the documentation of this file.
1 #include "dsdpsys.h"
2 #include "dsdpvec.h"
3 #include "dsdplapack.h"
4 #include "dsdpdatamat_impl.h"
5 
11 typedef enum {
12  Init=0,
13  Assemble=1,
14  Factored=2, /* fail to allocate required space */
15  Inverted=3, /* indefinity is detected */
16  ISymmetric=4
17 } MatStatus;
18 
19 typedef struct{
20  char UPLO;
21  int LDA;
22  double *val,*v2;
23  double *sscale;
24  double *workn;
25  int scaleit;
26  int n;
27  int owndata;
28  MatStatus status;
29 } dtrumat;
30 
31 static int DTRUMatView(void*);
32 
33 
34 #undef __FUNCT__
35 #define __FUNCT__ "DSDPLAPACKROUTINE"
36 static void dtruscalevec(double alpha, double v1[], double v2[], double v3[], int n){
37  int i;
38  for (i=0;i<n;i++){
39  v3[i] = (alpha*v1[i]*v2[i]);
40  }
41  return;
42 }
43 
44 static void dsydadd(double x[], double s[], double y[],int n){
45  int i;
46  for (i=0;i<n;i++){
47  y[i] += x[i]*(1/(s[i]*s[i])-2);
48  }
49  return;
50 }
51 
52 
53 static void dtruscalemat(double vv[], double ss[], int n,int LDA){
54  int i;
55  for (i=0;i<n;i++,vv+=LDA){
56  dtruscalevec(ss[i],vv,ss,vv,i+1);
57  }
58  return;
59 }
60 
61 static void DTRUMatOwnData(void* A, int owndata){
62  dtrumat* AA=(dtrumat*)A;
63  AA->owndata=owndata;
64  return;
65 }
66 
67 static int SUMatGetLDA(int n){
68  int nlda=n;
69  if (n>8 && nlda%2==1){ nlda++;}
70  if (n>100){
71  while (nlda%8!=0){ nlda++;}
72  }
73  /*
74  printf("LDA: %d %d %d \n",n,nlda,(int)sizeof(double));
75  */
76  return (nlda);
77 }
78 
79 static int DTRUMatCreateWData(int n,int LDA,double nz[], int nnz, dtrumat**M){
80  int i,info;
81  dtrumat*M23;
82  if (nnz<n*n){DSDPSETERR1(2,"Array must have length of : %d \n",n*n);}
83  DSDPCALLOC1(&M23,dtrumat,&info);DSDPCHKERR(info);
84  DSDPCALLOC2(&M23->sscale,double,n,&info);DSDPCHKERR(info);
85  DSDPCALLOC2(&M23->workn,double,n,&info);DSDPCHKERR(info);
86  M23->owndata=0; M23->val=nz; M23->n=n; M23->UPLO='U';M23->LDA=n;
87  M23->status=Init;
88  for (i=0;i<n;i++)M23->sscale[i]=1.0;
89  M23->scaleit=1;
90  M23->LDA=LDA;
91  if (n<=0){M23->LDA=1;}
92  *M=M23;
93  return 0;
94 }
95 
96 
97 
98 #undef __FUNCT__
99 #define __FUNCT__ "DSDPGetEigs"
100 int DSDPGetEigs(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
101  double W[],int n2,
102  double WORK[],int nd, int LLWORK[], int ni){
103  ffinteger N=n,LDA=n,LDZ=n,LWORK=nd,INFO=0;
104  char UPLO='U',JOBZ='V',RANGE='A';
105  /* Faster, but returns more error codes. ie. INFO>0 sometimes*/
106 
107  LDA=DSDPMax(1,n);
108  LDZ=DSDPMax(1,n);
109  if ( n2/2.5 > n || (ni<10*n+1) || (nd<26*n+1) || (nn0 < n*LDA) || (nn1<LDZ*n) ){
110  /*
111  printf("n: %d, ni: %d, nd: %d\n",n,ni/n,nd/n);
112  printf("SLOW VERSION\n");
113  */
114  dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
115  } else {
116 
117  int i;
118  ffinteger M,IL=1,IU=n,*ISUPPZ=(ffinteger*)IA;
119  ffinteger *IWORK=(ffinteger*)LLWORK,LIWORK=(ffinteger)ni;
120  double *Z=AA,VL=-1e10,VU=1e10,ABSTOL=0;
121  /* ABSTOL=dlamch_("Safe minimum" ); */
122  if (0==1){
123  dtrumat*M;
124  DTRUMatCreateWData(n,n,A,n*n,&M);
125  DTRUMatView((void*)M);
126  }
127  /*
128  printf("N: %d N2: %d , ",n,n2);
129  */
130  /*
131  LWORK=26*n; LIWORK=10*n;
132  */
133  /*
134  printf("JOBZ: %c, RANGE: %c, UPLO: %c, N: %d LDA: %d \n",
135  JOBZ,RANGE,UPLO, N,LDA);
136  printf("VL: %4.4e, VU: %4.4e, IL: %d, IU: %d, ABSTOL: %4.4e, LDZ: %d\n",
137  VL,VU,IL,IU,ABSTOL,LDZ);
138  printf("LWORK: %d, LIWORK: %d\n",LWORK,LIWORK);
139  */
140 
141  dsyevr(&JOBZ,&RANGE,&UPLO,&N,A,&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,Z,&LDZ,ISUPPZ,WORK,&LWORK,IWORK,&LIWORK,&INFO);
142  for (i=0;i<N*N;i++){A[i]=Z[i];}
143 
144  }
145  return INFO;
146 }
147 
148 #undef __FUNCT__
149 #define __FUNCT__ "DSDPGetEigsSTEP"
150 int DSDPGetEigsSTEP(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
151  double W[],int n2,
152  double WORK[],int nd, int LLWORK[], int ni){
153  int info;
154  info=DSDPGetEigs(A,n,AA,nn0,IA,nn1,W,n2,WORK,nd,LLWORK,ni);
155  return info;
156 }
157 
158 #undef __FUNCT__
159 #define __FUNCT__ "DSDPGetEigs2"
160 int DSDPGetEigs2(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
161  double W[],int n2,
162  double WORK[],int nd, int LLWORK[], int ni){
163  ffinteger N=n,LDA=n,LDZ=n,LWORK=nd,INFO=0;
164  char UPLO='U',JOBZ='V';
165  /* Works and uses less memory, but slower by a factor of about 2 or 3 */
166  LDA=DSDPMax(1,n);
167  LDZ=DSDPMax(1,n);
168  if (0==1){
169  dtrumat*M;
170  DTRUMatCreateWData(n,n,A,n*n,&M);
171  DTRUMatView((void*)M);
172  }
173  dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
174  return INFO;
175 }
176 
177 
178 #undef __FUNCT__
179 #define __FUNCT__ "DSDP FULL SYMMETRIC U LAPACK ROUTINE"
180 
181 static int DTRUMatMult(void* AA, double x[], double y[], int n){
182  dtrumat* A=(dtrumat*) AA;
183  ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
184  double BETA=0.0,ALPHA=1.0;
185  double *AP=A->val,*Y=y,*X=x;
186  char UPLO=A->UPLO,TRANS='N';
187 
188  if (A->n != n) return 1;
189  if (x==0 && n>0) return 3;
190  if (0==1){
191  dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
192  } else {
193  dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
194  }
195  return 0;
196 }
197 
198 
199 static int DTRUMatMultR(void* AA, double x[], double y[], int n){
200  dtrumat* A=(dtrumat*) AA;
201  ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
202  double ALPHA=1.0;
203  double *AP=A->val,*Y=y,*X=x,*ss=A->sscale,*W=A->workn;
204  char UPLO=A->UPLO,TRANS='N',DIAG='U';
205 
206  UPLO='L';
207  if (A->n != n) return 1;
208  if (x==0 && n>0) return 3;
209  /* dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY); */
210 
211  memset(y,0,n*sizeof(double));
212 
213  memcpy(W,X,n*sizeof(double));
214  TRANS='N'; UPLO='L';
215  dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
216  daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
217 
218  memcpy(W,X,n*sizeof(double));
219  TRANS='T'; UPLO='L';
220  dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
221  daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
222 
223  dsydadd(x,ss,y,n);
224  return 0;
225 }
226 
227 
228 static void DTRUMatScale(void*AA){
229  dtrumat* A=(dtrumat*) AA;
230  int i,N=A->n,LDA=A->LDA;
231  double *ss=A->sscale,*v=A->val;
232  for (i=0;i<N;i++){ ss[i]=*v;v+=(LDA+1);}
233  for (i=0;i<N;i++){ if (ss[i]!=0.0){ss[i]=1.0/sqrt(fabs(ss[i]));} else {ss[i]=1.0; }}
234  dtruscalemat(A->val,ss,N,LDA);
235 }
236 
237 static int DTRUMatCholeskyFactor(void* AA, int *flag){
238  dtrumat* A=(dtrumat*) AA;
239  ffinteger INFO,LDA=A->LDA,N=A->n;
240  double *AP=A->val;
241  char UPLO=A->UPLO;
242 
243  if (A->scaleit){ DTRUMatScale(AA);}
244  dpotrf(&UPLO, &N, AP, &LDA, &INFO );
245  *flag=INFO;
246  A->status=Factored;
247  return 0;
248 }
249 
250 
251 int dtrsm2(char*,char*,char*,char*,ffinteger*,ffinteger*,double*,double*,ffinteger*,double*,ffinteger*);
252 
253 static int DTRUMatSolve(void* AA, double b[], double x[],int n){
254  dtrumat* A=(dtrumat*) AA;
255  ffinteger ierr,INFO=0,NRHS=1,LDA=A->LDA,LDB=A->LDA,N=A->n;
256  double *AP=A->val,*ss=A->sscale,ONE=1.0;
257  char SIDE='L',UPLO=A->UPLO,TRANSA='N',DIAG='N';
258 
259  dtruscalevec(1.0,ss,b,x,n);
260 
261  if (0==1){
262  dpotrs(&UPLO, &N, &NRHS, AP, &LDA, x, &LDB, &INFO );
263  } else {
264  TRANSA='T';
265  ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
266  TRANSA='N';
267  ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
268  }
269 
270  dtruscalevec(1.0,ss,x,x,n);
271  return INFO;
272 }
273 
274 
275 static int DTRUMatShiftDiagonal(void* AA, double shift){
276  dtrumat* A=(dtrumat*) AA;
277  int i,n=A->n, LDA=A->LDA;
278  double *v=A->val;
279  for (i=0; i<n; i++){
280  v[i*LDA+i]+=shift;
281  }
282  return 0;
283 }
284 
285 static int DTRUMatAddRow(void* AA, int nrow, double dd, double row[], int n){
286  dtrumat* A=(dtrumat*) AA;
287  ffinteger ione=1,LDA=A->LDA,nn,INCX=1,INCY=A->LDA;
288  double *vv=A->val;
289 
290  nn=nrow;
291  daxpy(&nn,&dd,row,&INCX,vv+nrow,&INCY);
292  nn=nrow+1;
293  daxpy(&nn,&dd,row,&ione,vv+nrow*LDA,&ione);
294 
295  return 0;
296 }
297 
298 static int DTRUMatZero(void* AA){
299  dtrumat* A=(dtrumat*) AA;
300  int mn=A->n*(A->LDA);
301  double *vv=A->val;
302  memset((void*)vv,0,mn*sizeof(double));
303  A->status=Assemble;
304  return 0;
305 }
306 
307 static int DTRUMatGetSize(void *AA, int *n){
308  dtrumat* A=(dtrumat*) AA;
309  *n=A->n;
310  return 0;
311 }
312 
313 static int DTRUMatDestroy(void* AA){
314  int info;
315  dtrumat* A=(dtrumat*) AA;
316  if (A && A->owndata){DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
317  if (A && A->sscale) {DSDPFREE(&A->sscale,&info);DSDPCHKERR(info);}
318  if (A && A->workn) {DSDPFREE(&A->workn,&info);DSDPCHKERR(info);}
319  if (A) {DSDPFREE(&A,&info);DSDPCHKERR(info);}
320  return 0;
321 }
322 
323 static int DTRUMatView(void* AA){
324  dtrumat* M=(dtrumat*) AA;
325  int i,j;
326  double *val=M->val;
327  ffinteger LDA=M->LDA;
328  for (i=0; i<M->n; i++){
329  for (j=0; j<=i; j++){
330  printf(" %9.2e",val[i*LDA+j]);
331  }
332  for (j=i+1; j<M->LDA; j++){
333  printf(" %9.1e",val[i*LDA+j]);
334  }
335  printf("\n");
336  }
337  return 0;
338 }
339 
340 static int DTRUMatView2(void* AA){
341  dtrumat* M=(dtrumat*) AA;
342  int i,j;
343  double *val=M->v2;
344  ffinteger LDA=M->LDA;
345  for (i=0; i<M->n; i++){
346  for (j=0; j<=i; j++){
347  printf(" %9.2e",val[i*LDA+j]);
348  }
349  for (j=i+1; j<M->LDA; j++){
350  printf(" %9.2e",val[i*LDA+j]);
351  }
352  printf("\n");
353  }
354  return 0;
355 }
356 
357 
358 #include "dsdpschurmat_impl.h"
359 #include "dsdpdualmat_impl.h"
360 #include "dsdpdatamat_impl.h"
361 #include "dsdpxmat_impl.h"
362 #include "dsdpdsmat_impl.h"
363 #include "dsdpsys.h"
364 
365 #undef __FUNCT__
366 #define __FUNCT__ "Tassemble"
367 static int DTRUMatAssemble(void*M){
368  int info;
369  double shift=1.0e-15;
370  DSDPFunctionBegin;
371  info= DTRUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
372  DSDPFunctionReturn(0);
373 }
374 
375 static int DTRUMatRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
376  int i;
377  DSDPFunctionBegin;
378  *ncols = row+1;
379  for (i=0;i<=row;i++){
380  cols[i]=1.0;
381  }
382  memset((void*)(cols+row+1),0,(nrows-row-1)*sizeof(int));
383  DSDPFunctionReturn(0);
384 }
385 
386 
387 #undef __FUNCT__
388 #define __FUNCT__ "TAddDiag"
389 static int DTRUMatAddDiag(void*M, int row, double dd){
390  int ii;
391  dtrumat* ABA=(dtrumat*)M;
392  ffinteger LDA=ABA->LDA;
393  DSDPFunctionBegin;
394  ii=row*LDA+row;
395  ABA->val[ii] +=dd;
396  DSDPFunctionReturn(0);
397 }
398 #undef __FUNCT__
399 #define __FUNCT__ "TAddDiag2"
400 static int DTRUMatAddDiag2(void*M, double diag[], int m){
401  int row,ii;
402  dtrumat* ABA=(dtrumat*)M;
403  ffinteger LDA=ABA->LDA;
404  DSDPFunctionBegin;
405  for (row=0;row<m;row++){
406  ii=row*LDA+row;
407  ABA->val[ii] +=diag[row];
408  }
409  DSDPFunctionReturn(0);
410 }
411 static struct DSDPSchurMat_Ops dsdpmmatops;
412 static const char* lapackname="DENSE,SYMMETRIC U STORAGE";
413 
414 static int DSDPInitSchurOps(struct DSDPSchurMat_Ops* mops){
415  int info;
416  DSDPFunctionBegin;
417  info=DSDPSchurMatOpsInitialize(mops);DSDPCHKERR(info);
418  mops->matrownonzeros=DTRUMatRowNonzeros;
419  mops->matscaledmultiply=DTRUMatMult;
420  mops->matmultr=DTRUMatMultR;
421  mops->mataddrow=DTRUMatAddRow;
422  mops->mataddelement=DTRUMatAddDiag;
423  mops->matadddiagonal=DTRUMatAddDiag2;
424  mops->matshiftdiagonal=DTRUMatShiftDiagonal;
425  mops->matassemble=DTRUMatAssemble;
426  mops->matfactor=DTRUMatCholeskyFactor;
427  mops->matsolve=DTRUMatSolve;
428  mops->matdestroy=DTRUMatDestroy;
429  mops->matzero=DTRUMatZero;
430  mops->matview=DTRUMatView;
431  mops->id=1;
432  mops->matname=lapackname;
433  DSDPFunctionReturn(0);
434 }
435 
436 
437 #undef __FUNCT__
438 #define __FUNCT__ "DSDPGetLAPACKSUSchurOps"
439 int DSDPGetLAPACKSUSchurOps(int n,struct DSDPSchurMat_Ops** sops, void** mdata){
440  int info,nn,LDA;
441  double *vv;
442  dtrumat *AA;
443  DSDPFunctionBegin;
444 
445  LDA=SUMatGetLDA(n);
446  nn=n*LDA;
447  DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
448  info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
449  AA->owndata=1;
450  info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info);
451  *sops=&dsdpmmatops;
452  *mdata=(void*)AA;
453  DSDPFunctionReturn(0);
454 }
455 
456 
457 static int DTRUMatCholeskyBackward(void* AA, double b[], double x[], int n){
458  dtrumat* M=(dtrumat*) AA;
459  ffinteger N=M->n,INCX=1,LDA=M->LDA;
460  double *AP=M->val,*ss=M->sscale;
461  char UPLO=M->UPLO,TRANS='N',DIAG='N';
462  memcpy(x,b,N*sizeof(double));
463  dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
464  dtruscalevec(1.0,ss,x,x,n);
465  return 0;
466 }
467 
468 
469 static int DTRUMatCholeskyForward(void* AA, double b[], double x[], int n){
470  dtrumat* M=(dtrumat*) AA;
471  ffinteger N=M->n,INCX=1,LDA=M->LDA;
472  double *AP=M->val,*ss=M->sscale;
473  char UPLO=M->UPLO,TRANS='T',DIAG='N';
474  dtruscalevec(1.0,ss,b,x,n);
475  dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
476  return 0;
477 }
478 
479 static int DTRUMatLogDet(void* AA, double *dd){
480  dtrumat* A=(dtrumat*) AA;
481  int i,n=A->n,LDA=A->LDA;
482  double d=0,*v=A->val,*ss=A->sscale;
483  for (i=0; i<n; i++){
484  if (*v<=0) return 1;
485  d+=2*log(*v/ss[i]);
486  v+=LDA+1;
487  }
488  *dd=d;
489  return 0;
490 }
491 
492 
493 static int DTRUMatCholeskyForwardMultiply(void* AA, double x[], double y[], int n){
494  dtrumat* A=(dtrumat*) AA;
495  ffinteger i,j,N=A->n,LDA=A->LDA;
496  double *AP=A->val,*ss=A->sscale;
497  /* char UPLO=A->UPLO,TRANS='N',DIAG='N'; */
498 
499  if (x==0 && N>0) return 3;
500  /*
501  memcpy(y,x,N*sizeof(double));
502  dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,Y,&INCY);
503  */
504  for (i=0;i<N;i++)y[i]=0;
505  for (i=0;i<N;i++){
506  for (j=0;j<=i;j++){
507  y[i]+=AP[j]*x[j];
508  }
509  AP=AP+LDA;
510  }
511 
512  for (i=0;i<N;i++){ y[i]=y[i]/ss[i];}
513  return 0;
514 }
515 
516 static int DTRUMatCholeskyBackwardMultiply(void* AA, double x[], double y[], int n){
517  dtrumat* A=(dtrumat*) AA;
518  ffinteger i,j,N=A->n,LDA=A->LDA;
519  double *AP=A->val,*ss=A->sscale;
520  /* char UPLO=A->UPLO,TRANS='N',DIAG='N'; */
521 
522  if (x==0 && N>0) return 3;
523  /*
524  memcpy(y,x,N*sizeof(double));
525  dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,Y,&INCY);
526  */
527  for (i=0;i<N;i++)y[i]=0;
528  for (i=0;i<N;i++){
529  for (j=0;j<=i;j++){
530  y[j]+=AP[j]*x[i]/ss[i];
531  }
532  AP=AP+LDA;
533  }
534 
535  for (i=0;i<-N;i++){ y[i]=y[i]/ss[i];}
536  return 0;
537 }
538 
539 static int DTRUMatInvert(void* AA){
540  dtrumat* A=(dtrumat*) AA;
541  ffinteger INFO,LDA=A->LDA,N=A->n,nn=LDA*N;
542  double *v=A->val,*AP=A->v2,*ss=A->sscale;
543  char UPLO=A->UPLO;
544  memcpy((void*)AP,(void*)v,nn*sizeof(double));
545  dpotri(&UPLO, &N, AP, &LDA, &INFO );
546  if (INFO){
547  INFO=DTRUMatShiftDiagonal(AA,1e-11); INFO=0;
548  memcpy((void*)AP,(void*)v,nn*sizeof(double));
549  dpotri(&UPLO, &N, AP, &LDA, &INFO );
550  }
551  if (A->scaleit){
552  dtruscalemat(AP,ss,N,LDA);
553  }
554  A->status=Inverted;
555  return INFO;
556 
557 }
558 
559 static void DTRUSymmetrize(dtrumat* A){
560  int i,j,n=A->n,row,LDA=A->LDA;
561  double *v2=A->v2,*r1=A->v2,*r2=A->v2+LDA,*c1;
562  for (i=0;i<n/2;i++){
563  row=2*i;
564  r1=v2+LDA*(row);
565  r2=v2+LDA*(row+1);
566  c1=v2+LDA*(row+1);
567  r1[row+1]=c1[row];
568  c1+=LDA;
569  for (j=row+2;j<n;j++){
570  r1[j]=c1[row];
571  r2[j]=c1[row+1];
572  c1+=LDA;
573  }
574  r1+=LDA*2;
575  r2+=LDA*2;
576  }
577 
578  for (row=2*n/2;row<n;row++){
579  r1=v2+LDA*(row);
580  c1=v2+LDA*(row+1);
581  for (j=row+1;j<n;j++){
582  r1[j]=c1[row];
583  c1+=LDA;
584  }
585  r1+=LDA;
586  }
587  A->status=ISymmetric;
588  return;
589 }
590 
591 static int DTRUMatInverseAdd(void* AA, double alpha, double y[], int nn, int n){
592  dtrumat* A=(dtrumat*) AA;
593  ffinteger i,LDA=A->LDA,N,ione=1;
594  double *v2=A->v2;
595  for (i=0;i<n;i++){
596  N=i+1;
597  daxpy(&N,&alpha,v2,&ione,y,&ione);
598  v2+=LDA; y+=n;
599  }
600  return 0;
601 }
602 
603 static int DTRUMatInverseAddP(void* AA, double alpha, double y[], int nn, int n){
604  dtrumat* A=(dtrumat*) AA;
605  ffinteger N,LDA=A->LDA,i,ione=1;
606  double *v2=A->v2;
607  for (i=0;i<n;i++){
608  N=i+1;
609  daxpy(&N,&alpha,v2,&ione,y,&ione);
610  v2+=LDA; y+=i+1;
611  }
612  return 0;
613 }
614 
615 static void daddrow(double *v, double alpha, int i, double row[], int n){
616  double *s1;
617  ffinteger j,nn=n,ione=1;
618  if (alpha==0){return;}
619  nn=i+1; s1=v+i*n;
620  daxpy(&nn,&alpha,s1,&ione,row,&ione);
621  s1=v+i*n+n;
622  for (j=i+1;j<n;j++){ row[j]+=alpha*s1[i]; s1+=n; }
623  return;
624 }
625 /*
626 static void printrow(double r[], int n){int i;
627  for (i=0;i<n;i++){printf(" %4.2e",r[i]);} printf("\n"); }
628 */
629 static int DTRUMatInverseMultiply(void* AA, int indx[], int nind, double x[], double y[],int n){
630  dtrumat* A=(dtrumat*) AA;
631  ffinteger nn=n,LDA=A->LDA,N=A->n,INCX=1,INCY=1;
632  double *AP=A->v2,*s1=A->v2,*s2,*X=x,*Y=y,ALPHA=1.0,BETA=0.0;
633  char UPLO=A->UPLO,TRANS='N';
634  int i,ii,usefull=1;
635 
636  if (usefull){
637  if (A->status==Inverted){
638  DTRUSymmetrize(A);
639  }
640  if (nind < n/4){
641  memset((void*)y,0,n*sizeof(double));
642  for (ii=0;ii<nind;ii++){
643  i=indx[ii]; nn=n; ALPHA=x[i];s2=s1+i*LDA;
644  daxpy(&nn,&ALPHA,s2,&INCY,y,&INCX);
645  }
646  } else{
647  ALPHA=1.0;
648  dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
649  }
650 
651  } else {
652  if (nind<n/4 ){
653  memset((void*)y,0,n*sizeof(double));
654  for (ii=0;ii<nind;ii++){
655  i=indx[ii]; ALPHA=x[i];
656  daddrow(s1,ALPHA,i,y,n);
657  }
658  } else {
659  ALPHA=1.0;
660  dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
661  }
662  }
663  return 0;
664 }
665 
666 static int DTRUMatSetXMat(void*A, double v[], int nn, int n){
667  dtrumat* ABA=(dtrumat*)A;
668  int i,LDA=ABA->LDA;
669  double *vv=ABA->val;
670  if (vv!=v){
671  for (i=0;i<n;i++){
672  memcpy((void*)vv,(void*)v,(i+1)*sizeof(double));
673  vv+=LDA; v+=n;
674  }
675  }
676  ABA->status=Assemble;
677  return 0;
678 }
679 static int DTRUMatSetXMatP(void*A, double v[], int nn, int n){
680  dtrumat* ABA=(dtrumat*)A;
681  int i,LDA=ABA->LDA;
682  double *vv=ABA->val;
683  if (vv!=v){
684  for (i=0;i<n;i++){
685  memcpy((void*)vv,(void*)v,(i+1)*sizeof(double));
686  v+=(i+1); vv+=LDA;
687  }
688  }
689  ABA->status=Assemble;
690  return 0;
691 }
692 
693 static int DTRUMatFull(void*A, int*full){
694  *full=1;
695  return 0;
696 }
697 
698 static int DTRUMatGetArray(void*A,double **v,int *n){
699  dtrumat* ABA=(dtrumat*)A;
700  *n=ABA->n*ABA->LDA;
701  *v=ABA->val;
702  return 0;
703 }
704 
705 static struct DSDPDualMat_Ops sdmatops;
706 static int SDualOpsInitialize(struct DSDPDualMat_Ops* sops){
707  int info;
708  if (sops==NULL) return 0;
709  info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info);
710  sops->matseturmat=DTRUMatSetXMat;
711  sops->matgetarray=DTRUMatGetArray;
712  sops->matcholesky=DTRUMatCholeskyFactor;
713  sops->matsolveforward=DTRUMatCholeskyForward;
714  sops->matsolvebackward=DTRUMatCholeskyBackward;
715  sops->matinvert=DTRUMatInvert;
716  sops->matinverseadd=DTRUMatInverseAdd;
717  sops->matinversemultiply=DTRUMatInverseMultiply;
718  sops->matforwardmultiply=DTRUMatCholeskyForwardMultiply;
719  sops->matbackwardmultiply=DTRUMatCholeskyBackwardMultiply;
720  sops->matfull=DTRUMatFull;
721  sops->matdestroy=DTRUMatDestroy;
722  sops->matgetsize=DTRUMatGetSize;
723  sops->matview=DTRUMatView;
724  sops->matlogdet=DTRUMatLogDet;
725  sops->matname=lapackname;
726  sops->id=1;
727  return 0;
728 }
729 
730 
731 #undef __FUNCT__
732 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
733 static int DSDPLAPACKSUDualMatCreate(int n,struct DSDPDualMat_Ops **sops, void**smat){
734  dtrumat *AA;
735  int info,nn,LDA=n;
736  double *vv;
737  DSDPFunctionBegin;
738  LDA=SUMatGetLDA(n);
739  nn=n*LDA;
740  DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
741  info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
742  AA->owndata=1;
743  info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
744  *sops=&sdmatops;
745  *smat=(void*)AA;
746  DSDPFunctionReturn(0);
747 }
748 
749 
750 static int switchptr(void *SD,void *SP){
751  dtrumat *s1,*s2;
752  s1=(dtrumat*)(SD);
753  s2=(dtrumat*)(SP);
754  s1->v2=s2->val;
755  s2->v2=s1->val;
756  return 0;
757 }
758 
759 
760 #undef __FUNCT__
761 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate2"
762 int DSDPLAPACKSUDualMatCreate2(int n,
763  struct DSDPDualMat_Ops **sops1, void**smat1,
764  struct DSDPDualMat_Ops **sops2, void**smat2){
765  int info;
766  DSDPFunctionBegin;
767  info=DSDPLAPACKSUDualMatCreate(n,sops1,smat1);DSDPCHKERR(info);
768  info=DSDPLAPACKSUDualMatCreate(n,sops2,smat2);DSDPCHKERR(info);
769  info=switchptr(*smat1,*smat2);DSDPCHKERR(info);
770  DSDPFunctionReturn(0);
771 }
772 
773 static struct DSDPDualMat_Ops sdmatopsp;
774 static int SDualOpsInitializeP(struct DSDPDualMat_Ops* sops){
775  int info;
776  if (sops==NULL) return 0;
777  info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info);
778  sops->matseturmat=DTRUMatSetXMatP;
779  sops->matgetarray=DTRUMatGetArray;
780  sops->matcholesky=DTRUMatCholeskyFactor;
781  sops->matsolveforward=DTRUMatCholeskyForward;
782  sops->matsolvebackward=DTRUMatCholeskyBackward;
783  sops->matinvert=DTRUMatInvert;
784  sops->matinverseadd=DTRUMatInverseAddP;
785  sops->matinversemultiply=DTRUMatInverseMultiply;
786  sops->matforwardmultiply=DTRUMatCholeskyForwardMultiply;
787  sops->matbackwardmultiply=DTRUMatCholeskyBackwardMultiply;
788  sops->matfull=DTRUMatFull;
789  sops->matdestroy=DTRUMatDestroy;
790  sops->matgetsize=DTRUMatGetSize;
791  sops->matview=DTRUMatView;
792  sops->matlogdet=DTRUMatLogDet;
793  sops->matname=lapackname;
794  sops->id=1;
795  return 0;
796 }
797 
798 #undef __FUNCT__
799 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
800 static int DSDPLAPACKSUDualMatCreateP(int n, struct DSDPDualMat_Ops **sops, void**smat){
801  dtrumat *AA;
802  int info,nn,LDA;
803  double *vv;
804  DSDPFunctionBegin;
805  LDA=SUMatGetLDA(n);
806  nn=LDA*n;
807  DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
808  info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
809  AA->owndata=1;
810  info=SDualOpsInitializeP(&sdmatopsp);DSDPCHKERR(info);
811  *sops=&sdmatopsp;
812  *smat=(void*)AA;
813  DSDPFunctionReturn(0);
814 }
815 
816 
817 #undef __FUNCT__
818 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate2P"
819 int DSDPLAPACKSUDualMatCreate2P(int n,
820  struct DSDPDualMat_Ops* *sops1, void**smat1,
821  struct DSDPDualMat_Ops* *sops2, void**smat2){
822  int info;
823  DSDPFunctionBegin;
824  info=DSDPLAPACKSUDualMatCreateP(n,sops1,smat1);
825  info=DSDPLAPACKSUDualMatCreateP(n,sops2,smat2);
826  info=switchptr(*smat1,*smat2);
827  DSDPFunctionReturn(0);
828 }
829 
830 static int DTRUMatScaleDiagonal(void* AA, double dd){
831  dtrumat* A=(dtrumat*) AA;
832  ffinteger LDA=A->LDA;
833  int i,n=A->n;
834  double *v=A->val;
835  for (i=0; i<n; i++){
836  *v*=dd;
837  v+=LDA+1;
838  }
839  return 0;
840 }
841 
842 static int DTRUMatOuterProduct(void* AA, double alpha, double x[], int n){
843  dtrumat* A=(dtrumat*) AA;
844  ffinteger ione=1,N=n,LDA=A->LDA;
845  double *v=A->val;
846  char UPLO=A->UPLO;
847  dsyr(&UPLO,&N,&alpha,x,&ione,v,&LDA);
848  return 0;
849 }
850 
851 static int DenseSymPSDNormF2(void* AA, int n, double *dddot){
852  dtrumat* A=(dtrumat*) AA;
853  ffinteger ione=1,nn=A->n*A->n;
854  double dd,tt=sqrt(0.5),*val=A->val;
855  int info;
856  info=DTRUMatScaleDiagonal(AA,tt);
857  dd=dnrm2(&nn,val,&ione);
858  info=DTRUMatScaleDiagonal(AA,1.0/tt);
859  *dddot=dd*dd*2;
860  return 0;
861 }
862 
863 static int DTRUMatGetDenseArray(void* A, double *v[], int*n){
864  dtrumat* ABA=(dtrumat*)A;
865  *v=ABA->val;
866  *n=ABA->n*ABA->LDA;
867  return 0;
868 }
869 
870 static int DTRUMatRestoreDenseArray(void* A, double *v[], int *n){
871  *v=0;*n=0;
872  return 0;
873 }
874 
875 static int DDenseSetXMat(void*A, double v[], int nn, int n){
876  dtrumat* ABA=(dtrumat*)A;
877  int i,LDA=ABA->LDA;
878  double *vv=ABA->val;
879  if (v!=vv){
880  for (i=0;i<n;i++){
881  memcpy((void*)vv,(void*)v,nn*sizeof(double));
882  v+=n;vv+=LDA;
883  }
884  }
885  ABA->status=Assemble;
886  return 0;
887 }
888 
889 static int DDenseVecVec(void* A, double x[], int n, double *v){
890  dtrumat* ABA=(dtrumat*)A;
891  int i,j,k=0,LDA=ABA->LDA;
892  double dd=0,*val=ABA->val;
893  *v=0.0;
894  for (i=0; i<n; i++){
895  for (j=0;j<i;j++){
896  dd+=2*x[i]*x[j]*val[j];
897  k++;
898  }
899  dd+=x[i]*x[i]*val[i];
900  k+=LDA;
901  }
902  *v=dd;
903  return 0;
904 }
905 
906 
907 
908 static int DTRUMatEigs(void*AA,double W[],double IIWORK[], int nn1, double *mineig){
909  dtrumat* AAA=(dtrumat*) AA;
910  ffinteger info,INFO=0,M,N=AAA->n;
911  ffinteger IL=1,IU=1,LDA=AAA->LDA,LDZ=LDA,LWORK,IFAIL;
912  ffinteger *IWORK=(ffinteger*)IIWORK;
913  double *AP=AAA->val;
914  double Z=0,VL=-1e10,VU=1e10,*WORK,ABSTOL=1e-13;
915  char UPLO=AAA->UPLO,JOBZ='N',RANGE='I';
916  DSDPCALLOC2(&WORK,double,8*N,&info);
917  DSDPCALLOC2(&IWORK,ffinteger,5*N,&info);
918  LWORK=8*N;
919  dsyevx(&JOBZ,&RANGE,&UPLO,&N,AP,&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,&LWORK,IWORK,&IFAIL,&INFO);
920  /*
921  ffinteger LIWORK=nn1;
922  dsyevd(&JOBZ,&UPLO,&N,AP,&LDA,W,WORK,&LWORK,IWORK,&LIWORK,&INFO);
923  */
924  DSDPFREE(&WORK,&info);
925  DSDPFREE(&IWORK,&info);
926  *mineig=W[0];
927  return INFO;
928 }
929 
930 
931 static struct DSDPVMat_Ops turdensematops;
932 
933 static int DSDPDenseXInitializeOps(struct DSDPVMat_Ops* densematops){
934  int info;
935  if (!densematops) return 0;
936  info=DSDPVMatOpsInitialize(densematops); DSDPCHKERR(info);
937  densematops->matview=DTRUMatView;
938  densematops->matscalediagonal=DTRUMatScaleDiagonal;
939  densematops->matshiftdiagonal=DTRUMatShiftDiagonal;
940  densematops->mataddouterproduct=DTRUMatOuterProduct;
941  densematops->matmult=DTRUMatMult;
942  densematops->matdestroy=DTRUMatDestroy;
943  densematops->matfnorm2=DenseSymPSDNormF2;
944  densematops->matgetsize=DTRUMatGetSize;
945  densematops->matzeroentries=DTRUMatZero;
946  densematops->matgeturarray=DTRUMatGetDenseArray;
947  densematops->matrestoreurarray=DTRUMatRestoreDenseArray;
948  densematops->matmineig=DTRUMatEigs;
949  densematops->id=1;
950  densematops->matname=lapackname;
951  return 0;
952 }
953 
954 #undef __FUNCT__
955 #define __FUNCT__ "DSDPXMatUCreateWithData"
956 int DSDPXMatUCreateWithData(int n,double nz[],int nnz,struct DSDPVMat_Ops* *xops, void * *xmat){
957  int i,info;
958  double dtmp;
959  dtrumat*AA;
960  DSDPFunctionBegin;
961  if (nnz<n*n){DSDPSETERR1(2,"Array must have length of : %d \n",n*n);}
962  for (i=0;i<n*n;i++) dtmp=nz[i];
963  info=DTRUMatCreateWData(n,n,nz,nnz,&AA); DSDPCHKERR(info);
964  AA->owndata=0;
965  info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
966  *xops=&turdensematops;
967  *xmat=(void*)AA;
968  DSDPFunctionReturn(0);
969 }
970 
971 #undef __FUNCT__
972 #define __FUNCT__ "DSDPXMatUCreate"
973 int DSDPXMatUCreate(int n,struct DSDPVMat_Ops* *xops, void * *xmat){
974  int info,nn=n*n;
975  double *vv;
976  DSDPFunctionBegin;
977  DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
978  info=DSDPXMatUCreateWithData(n,vv,nn,xops,xmat);DSDPCHKERR(info);
979  DTRUMatOwnData((dtrumat*)(*xmat),1);
980  DSDPFunctionReturn(0);
981 }
982 
983 static struct DSDPDSMat_Ops tdsdensematops;
984 static int DSDPDSDenseInitializeOps(struct DSDPDSMat_Ops* densematops){
985  int info;
986  if (!densematops) return 0;
987  info=DSDPDSMatOpsInitialize(densematops); DSDPCHKERR(info);
988  densematops->matseturmat=DDenseSetXMat;
989  densematops->matview=DTRUMatView;
990  densematops->matdestroy=DTRUMatDestroy;
991  densematops->matgetsize=DTRUMatGetSize;
992  densematops->matzeroentries=DTRUMatZero;
993  densematops->matmult=DTRUMatMult;
994  densematops->matvecvec=DDenseVecVec;
995  densematops->id=1;
996  densematops->matname=lapackname;
997  return 0;
998 }
999 
1000 #undef __FUNCT__
1001 #define __FUNCT__ "DSDPCreateDSMatWithArray2"
1002 int DSDPCreateDSMatWithArray2(int n,double vv[],int nnz,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
1003  int info;
1004  dtrumat*AA;
1005  DSDPFunctionBegin;
1006  info=DTRUMatCreateWData(n,n,vv,nnz,&AA); DSDPCHKERR(info);
1007  AA->owndata=0;
1008  info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
1009  *dsmatops=&tdsdensematops;
1010  *dsmat=(void*)AA;
1011  DSDPFunctionReturn(0);
1012 }
1013 
1014 #undef __FUNCT__
1015 #define __FUNCT__ "DSDPCreateXDSMat2"
1016 int DSDPCreateXDSMat2(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
1017  int info,nn=n*n;
1018  double *vv;
1019  DSDPFunctionBegin;
1020  DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
1021  info=DSDPCreateDSMatWithArray2(n,vv,nn,dsmatops,dsmat);DSDPCHKERR(info);
1022  DTRUMatOwnData((dtrumat*)(*dsmat),1);
1023  DSDPFunctionReturn(0);
1024 }
1025 
1026 
1027 typedef struct {
1028  int neigs;
1029  double *eigval;
1030  double *an;
1031 } Eigen;
1032 
1033 typedef struct {
1034  dtrumat* AA;
1035  Eigen *Eig;
1036 } dvecumat;
1037 
1038 #undef __FUNCT__
1039 #define __FUNCT__ "CreateDvecumatWData"
1040 static int CreateDvecumatWdata(int n, double vv[], dvecumat **A){
1041  int info,nnz=n*n;
1042  dvecumat* V;
1043  DSDPCALLOC1(&V,dvecumat,&info);DSDPCHKERR(info);
1044  info=DTRUMatCreateWData(n,n,vv,nnz,&V->AA); DSDPCHKERR(info);
1045  V->Eig=0;
1046  *A=V;
1047  return 0;
1048 }
1049 
1050 
1051 static int DvecumatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int n){
1052  int k;
1053  *nnzz=n;
1054  for (k=0;k<n;k++) nz[k]++;
1055  return 0;
1056 }
1057 
1058 static int DTRUMatGetRowAdd(void* AA, int nrow, double ytmp, double row[], int n){
1059  dtrumat* A=(dtrumat*) AA;
1060  ffinteger i,nnn=n;
1061  double *v=A->val;
1062 
1063  nnn=nrow*n;
1064  for (i=0;i<=nrow;i++){
1065  row[i]+=ytmp*v[nnn+i];
1066  }
1067  for (i=nrow+1;i<n;i++){
1068  nnn+=nrow;
1069  row[i]+=ytmp*v[nrow];
1070  }
1071  return 0;
1072 }
1073 
1074 static int DvecumatGetRowAdd(void* AA, int trow, double scl, double r[], int m){
1075  int info;
1076  dvecumat* A=(dvecumat*)AA;
1077  info=DTRUMatGetRowAdd((void*)A->AA ,trow,scl,r,m);
1078  return 0;
1079 }
1080 
1081 static int DvecumatAddMultiple(void* AA, double alpha, double r[], int nnn, int n){
1082  dvecumat* A=(dvecumat*)AA;
1083  ffinteger nn=nnn, ione=1;
1084  double *val=A->AA->val;
1085  daxpy(&nn,&alpha,val,&ione,r,&ione);
1086  return 0;
1087 }
1088 
1089 
1090 static int DvecuEigVecVec(void*, double[], int, double*);
1091 static int DvecumatVecVec(void* AA, double x[], int n, double *v){
1092  dvecumat* A=(dvecumat*)AA;
1093  int i,j,k=0,LDA=A->AA->LDA;
1094  double dd=0,*val=A->AA->val;
1095  *v=0.0;
1096  if (A->Eig && A->Eig->neigs<n/5){
1097  i=DvecuEigVecVec(AA,x,n,v);
1098  } else {
1099  for (i=0; i<n; i++){
1100  for (j=0;j<i;j++){
1101  dd+=2*x[i]*x[j]*val[j];
1102  }
1103  dd+=x[i]*x[i]*val[i];
1104  k+=LDA;
1105  }
1106  *v=dd;
1107  }
1108  return 0;
1109 }
1110 
1111 
1112 static int DvecumatFNorm2(void* AA, int n, double *v){
1113  dvecumat* A=(dvecumat*)AA;
1114  long int i,j,k=0,LDA=A->AA->LDA;
1115  double dd=0,*x=A->AA->val;
1116  for (i=0; i<n; i++){
1117  for (j=0;j<i;j++){
1118  dd+=2*x[j]*x[j];
1119  }
1120  dd+=x[i]*x[i];
1121  k+=LDA;
1122  }
1123  *v=dd;
1124  return 0;
1125 }
1126 
1127 
1128 static int DvecumatCountNonzeros(void* AA, int *nnz, int n){
1129  *nnz=n*(n+1)/2;
1130  return 0;
1131 }
1132 
1133 
1134 static int DvecumatDot(void* AA, double x[], int nn, int n, double *v){
1135  dvecumat* A=(dvecumat*)AA;
1136  double d1,dd=0,*v1=x,*v2=A->AA->val;
1137  ffinteger i,n2,ione=1,LDA=A->AA->LDA;
1138 
1139  for (i=0;i<n;i++){
1140  n2=i+1;
1141  d1=ddot(&n2,v1,&ione,v2,&ione);
1142  v1+=n; v2+=LDA;
1143  dd+=d1;
1144  }
1145  *v=2*dd;
1146  return 0;
1147 }
1148 
1149 /*
1150 static int DvecumatNormF2(void* AA, int n, double *v){
1151  dvecumat* A=(dvecumat*)AA;
1152  return(DTRUMatNormF2((void*)(A->AA), n,v));
1153 }
1154 */
1155 #undef __FUNCT__
1156 #define __FUNCT__ "DvecumatDestroy"
1157 static int DvecumatDestroy(void* AA){
1158  dvecumat* A=(dvecumat*)AA;
1159  int info;
1160  info=DTRUMatDestroy((void*)(A->AA));
1161  if (A->Eig){
1162  DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
1163  DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
1164  }
1165  DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
1166  DSDPFREE(&A,&info);DSDPCHKERR(info);
1167  return 0;
1168 }
1169 
1170 
1171 static int DvecumatView(void* AA){
1172  dvecumat* A=(dvecumat*)AA;
1173  dtrumat* M=A->AA;
1174  int i,j,LDA=M->LDA;
1175  double *val=M->val;
1176  for (i=0; i<M->n; i++){
1177  for (j=0; j<M->n; j++){
1178  printf(" %4.2e",val[j]);
1179  }
1180  val+=LDA;
1181  }
1182  return 0;
1183 }
1184 
1185 
1186 #undef __FUNCT__
1187 #define __FUNCT__ "DSDPCreateDvecumatEigs"
1188 static int CreateEigenLocker(Eigen **EE,int neigs, int n){
1189  int info;
1190  Eigen *E;
1191 
1192  DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
1193  DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info);
1194  DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info);
1195  E->neigs=neigs;
1196  *EE=E;
1197  return 0;
1198 }
1199 
1200 
1201 static int EigMatSetEig(Eigen* A,int row, double eigv, double v[], int n){
1202  double *an=A->an;
1203  A->eigval[row]=eigv;
1204  memcpy((void*)(an+n*row),(void*)v,n*sizeof(double));
1205  return 0;
1206 }
1207 
1208 
1209 static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n){
1210  double* an=A->an;
1211  *eigenvalue=A->eigval[row];
1212  memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
1213  return 0;
1214 }
1215 
1216 
1217 static int DvecumatComputeEigs(dvecumat*,double[],int,double[],int,double[],int,int[],int);
1218 
1219 static int DvecumatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){
1220 
1221  int info;
1222  dvecumat* A=(dvecumat*)AA;
1223  if (A->Eig) return 0;
1224  info=DvecumatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2);DSDPCHKERR(info);
1225  return 0;
1226 }
1227 
1228 static int DvecumatGetRank(void *AA,int *rank, int n){
1229  dvecumat* A=(dvecumat*)AA;
1230  if (A->Eig){
1231  *rank=A->Eig->neigs;
1232  } else {
1233  DSDPSETERR(1,"Vecu Matrix not factored yet\n");
1234  }
1235  return 0;
1236 }
1237 
1238 static int DvecumatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indz[], int *nind){
1239  dvecumat* A=(dvecumat*)AA;
1240  int i,info;
1241  if (A->Eig){
1242  info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n);DSDPCHKERR(info);
1243  *nind=n;
1244  for (i=0;i<n;i++){ indz[i]=i;}
1245  } else {
1246  DSDPSETERR(1,"Vecu Matrix not factored yet\n");
1247  }
1248  return 0;
1249 }
1250 
1251 static int DvecuEigVecVec(void* AA, double v[], int n, double *vv){
1252  dvecumat* A=(dvecumat*)AA;
1253  int i,rank,neigs;
1254  double *an,dd,ddd=0,*eigval;
1255  if (A->Eig){
1256  an=A->Eig->an;
1257  neigs=A->Eig->neigs;
1258  eigval=A->Eig->eigval;
1259  for (rank=0;rank<neigs;rank++){
1260  for (dd=0,i=0;i<n;i++){
1261  dd+=v[i]*an[i];
1262  }
1263  an+=n;
1264  ddd+=dd*dd*eigval[rank];
1265  }
1266  *vv=ddd;
1267  } else {
1268  DSDPSETERR(1,"Vecu Matrix not factored yet\n");
1269  }
1270  return 0;
1271 }
1272 
1273 
1274 static struct DSDPDataMat_Ops dvecumatops;
1275 static const char *datamatname="STANDARD VECU MATRIX";
1276 
1277 static int DvecumatOpsInitialize(struct DSDPDataMat_Ops *sops){
1278  int info;
1279  if (sops==NULL) return 0;
1280  info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
1281  sops->matvecvec=DvecumatVecVec;
1282  sops->matdot=DvecumatDot;
1283  sops->mataddrowmultiple=DvecumatGetRowAdd;
1284  sops->mataddallmultiple=DvecumatAddMultiple;
1285  sops->matview=DvecumatView;
1286  sops->matdestroy=DvecumatDestroy;
1287  sops->matfactor2=DvecumatFactor;
1288  sops->matgetrank=DvecumatGetRank;
1289  sops->matgeteig=DvecumatGetEig;
1290  sops->matrownz=DvecumatGetRowNnz;
1291  sops->matfnorm2=DvecumatFNorm2;
1292  sops->matnnz=DvecumatCountNonzeros;
1293  sops->id=1;
1294  sops->matname=datamatname;
1295  return 0;
1296 }
1297 
1298 #undef __FUNCT__
1299 #define __FUNCT__ "DSDPGetDUmat"
1300 int DSDPGetDUMat(int n,double *val, struct DSDPDataMat_Ops**sops, void**smat){
1301  int info,k;
1302  double dtmp;
1303  dvecumat* A;
1304  DSDPFunctionBegin;
1305 
1306  for (k=0;k<n*n;++k) dtmp=val[k];
1307  info=CreateDvecumatWdata(n,val,&A); DSDPCHKERR(info);
1308  A->Eig=0;
1309  info=DvecumatOpsInitialize(&dvecumatops); DSDPCHKERR(info);
1310  if (sops){*sops=&dvecumatops;}
1311  if (smat){*smat=(void*)A;}
1312  DSDPFunctionReturn(0);
1313 }
1314 
1315 
1316 #undef __FUNCT__
1317 #define __FUNCT__ "DvecumatComputeEigs"
1318 static int DvecumatComputeEigs(dvecumat* AA,double DD[], int nn0, double W[], int n, double WORK[], int n1, int iiptr[], int n2){
1319 
1320  int i,neigs,info;
1321  long int *i2darray=(long int*)DD;
1322  int ownarray1=0,ownarray2=0,ownarray3=0;
1323  double *val=AA->AA->val;
1324  double *dmatarray=0,*dworkarray=0,eps=1.0e-12;
1325  int nn1=0,nn2=0;
1326 
1327  /* create a dense array in which to put numbers */
1328  if (n*n>nn1){
1329  DSDPCALLOC2(&dmatarray,double,(n*n),&info); DSDPCHKERR(info);
1330  ownarray1=1;
1331  }
1332  memcpy((void*)dmatarray,(void*)val,n*n*sizeof(double));
1333 
1334  if (n*n>nn2){
1335  DSDPCALLOC2(&dworkarray,double,(n*n),&info); DSDPCHKERR(info);
1336  ownarray2=1;
1337  }
1338 
1339  if (n*n*sizeof(long int)>nn0*sizeof(double)){
1340  DSDPCALLOC2(&i2darray,long int,(n*n),&info); DSDPCHKERR(info);
1341  ownarray3=1;
1342  }
1343 
1344 
1345  /* Call LAPACK to compute the eigenvalues */
1346  info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n,
1347  W,n,WORK,n1,iiptr,n2);
1348  if (info){
1349  memcpy((void*)dmatarray,(void*)val,n*n*sizeof(double));
1350  info=DSDPGetEigs2(dmatarray,n,dworkarray,n*n,i2darray,n*n,
1351  W,n,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
1352  }
1353 
1354  /* Count the nonzero eigenvalues */
1355  for (neigs=0,i=0;i<n;i++){
1356  if (fabs(W[i])> eps ){ neigs++;}
1357  }
1358 
1359  info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info);
1360 
1361  /* Copy into structure */
1362  for (neigs=0,i=0; i<n; i++){
1363  if (fabs(W[i]) > eps){
1364  info=EigMatSetEig(AA->Eig,neigs,W[i],dmatarray+n*i,n);DSDPCHKERR(info);
1365  neigs++;
1366  }
1367  }
1368 
1369  if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
1370  if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
1371  if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
1372  return 0;
1373 }
1374 
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