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,...
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