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
11typedef 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
19typedef 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
31static int DTRUMatView(void*);
32
33
34#undef __FUNCT__
35#define __FUNCT__ "DSDPLAPACKROUTINE"
36static 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
44static 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
53static 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
61static void DTRUMatOwnData(void* A, int owndata){
62 dtrumat* AA=(dtrumat*)A;
63 AA->owndata=owndata;
64 return;
65}
66
67static 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
79static 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"
100int 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"
150int 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"
160int 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
181static 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
199static 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
228static 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
237static 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
251int dtrsm2(char*,char*,char*,char*,ffinteger*,ffinteger*,double*,double*,ffinteger*,double*,ffinteger*);
252
253static 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
275static 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
285static 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
298static 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
307static int DTRUMatGetSize(void *AA, int *n){
308 dtrumat* A=(dtrumat*) AA;
309 *n=A->n;
310 return 0;
311}
312
313static 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
323static 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
340static 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"
367static 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
375static 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"
389static 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"
400static 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}
411static struct DSDPSchurMat_Ops dsdpmmatops;
412static const char* lapackname="DENSE,SYMMETRIC U STORAGE";
413
414static 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"
439int 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
457static 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
469static 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
479static 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
493static 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
516static 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
539static 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
559static 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
591static 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
603static 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
615static 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/*
626static void printrow(double r[], int n){int i;
627 for (i=0;i<n;i++){printf(" %4.2e",r[i]);} printf("\n"); }
628*/
629static 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
666static 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}
679static 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
693static int DTRUMatFull(void*A, int*full){
694 *full=1;
695 return 0;
696}
697
698static 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
705static struct DSDPDualMat_Ops sdmatops;
706static 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"
733static 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
750static 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"
762int 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
773static struct DSDPDualMat_Ops sdmatopsp;
774static 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"
800static 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"
819int 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
830static 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
842static 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
851static 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
863static 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
870static int DTRUMatRestoreDenseArray(void* A, double *v[], int *n){
871 *v=0;*n=0;
872 return 0;
873}
874
875static 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
889static 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
908static 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
931static struct DSDPVMat_Ops turdensematops;
932
933static 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"
956int 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"
973int 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
983static struct DSDPDSMat_Ops tdsdensematops;
984static 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"
1002int 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"
1016int 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
1027typedef struct {
1028 int neigs;
1029 double *eigval;
1030 double *an;
1031} Eigen;
1032
1033typedef struct {
1034 dtrumat* AA;
1035 Eigen *Eig;
1036} dvecumat;
1037
1038#undef __FUNCT__
1039#define __FUNCT__ "CreateDvecumatWData"
1040static 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
1051static 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
1058static 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
1074static 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
1081static 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
1090static int DvecuEigVecVec(void*, double[], int, double*);
1091static 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
1112static 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
1128static int DvecumatCountNonzeros(void* AA, int *nnz, int n){
1129 *nnz=n*(n+1)/2;
1130 return 0;
1131}
1132
1133
1134static 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/*
1150static 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"
1157static 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
1171static 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"
1188static 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
1201static 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
1209static 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
1217static int DvecumatComputeEigs(dvecumat*,double[],int,double[],int,double[],int,int[],int);
1218
1219static 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
1228static 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
1238static 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
1251static 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
1274static struct DSDPDataMat_Ops dvecumatops;
1275static const char *datamatname="STANDARD VECU MATRIX";
1276
1277static 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"
1300int 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"
1318static 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 DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition dsdpdatamat.c:47
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
int DSDPDSMatOpsInitialize(struct DSDPDSMat_Ops *aops)
Set pointers to null.
Definition dsdpdsmat.c:214
Structure of function pointers that each SDP Delta S matrix type (sparse, dense, diagonal,...
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Structure of function pointers that each symmetric positive definite matrix type (dense,...
DSDP uses BLAS and LAPACK for many of its operations.
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide.
Error handling, printing, and profiling.
Vector operations used by the solver.
int DSDPVMatOpsInitialize(struct DSDPVMat_Ops *aops)
Set function pointers to null.
Definition dsdpxmat.c:377
Structure of function pointers that each dense matrix array type (upper full, packed symmetric,...
Symmetric Delta S matrix for one block in the semidefinite cone.
Table of function pointers that operate on the data matrix.
Table of function pointers that operate on the S matrix.
Table of function pointers that operate on the dense matrix.