20 static const char* lapackname=
"DENSE,SYMMETRIC,PACKED STORAGE";
22 int DTPUMatEigs(
void*AA,
double W[],
double IIWORK[],
int nn1,
double *mineig){
23 dtpumat* AAA=(dtpumat*) AA;
24 ffinteger info,INFO=0,M,N=AAA->n;
25 ffinteger IL=1,IU=1,LDZ=1,IFAIL;
26 ffinteger *IWORK=(ffinteger*)IIWORK;
27 double *AP=AAA->val,ABSTOL=1e-13;
28 double Z=0,VL=-1e10,VU=1;
30 char UPLO=AAA->UPLO,JOBZ=
'N',RANGE=
'I';
32 DSDPCALLOC2(&WORK,
double,7*N,&info);DSDPCHKERR(info);
33 DSDPCALLOC2(&IWORK,ffinteger,5*N,&info);DSDPCHKERR(info);
34 dspevx(&JOBZ,&RANGE,&UPLO,&N,AP,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,IWORK,&IFAIL,&INFO);
42 DSDPFREE(&WORK,&info);DSDPCHKERR(info);
43 DSDPFREE(&IWORK,&info);DSDPCHKERR(info);
49 #define __FUNCT__ "DSDPLAPACKROUTINE" 50 static void dtpuscalevec(
double alpha,
double v1[],
double v2[],
double v3[],
int n){
53 v3[i] = (alpha*v1[i]*v2[i]);
57 static void dtpuscalemat(
double vv[],
double ss[],
int n){
59 for (i=0;i<n;i++,vv+=i){
60 dtpuscalevec(ss[i],vv,ss,vv,i+1);
64 static int DTPUMatCreateWData(
int n,
double nz[],
int nnz, dtpumat**M){
65 int i,nn=(n*n+n)/2,info;
68 if (nnz<nn){DSDPSETERR1(2,
"Array must have length of : %d \n",nn);}
69 for (i=0;i<nnz;i++) dtmp=nz[i];
70 DSDPCALLOC1(&M23,dtpumat,&info);DSDPCHKERR(info);
71 DSDPCALLOC2(&M23->sscale,
double,n,&info);DSDPCHKERR(info);
72 M23->owndata=0; M23->val=nz; M23->n=n; M23->UPLO=
'U';
73 for (i=0;i<n;i++)M23->sscale[i]=1.0;
82 #define __FUNCT__ "DSDPLAPACK ROUTINE" 85 static int DTPUMatMult(
void* AA,
double x[],
double y[],
int n){
86 dtpumat* A=(dtpumat*) AA;
88 double BETA=0.0,ALPHA=1.0;
89 double *AP=A->val,*Y=y,*X=x;
92 if (A->n != n)
return 1;
93 if (x==0 && n>0)
return 3;
94 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione);
98 static int DTPUMatSolve(
void* AA,
double b[],
double x[],
int n){
99 dtpumat* A=(dtpumat*) AA;
100 ffinteger INFO,NRHS=1,LDB=A->n,N=A->n;
101 double *AP=A->val,*ss=A->sscale;
105 dtpuscalevec(1.0,ss,b,x,n);
106 dpptrs(&UPLO, &N, &NRHS, AP, x, &LDB, &INFO );
107 dtpuscalevec(1.0,ss,x,x,n);
111 static int DTPUMatCholeskyFactor(
void* AA,
int *flag){
112 dtpumat* A=(dtpumat*) AA;
114 ffinteger INFO,LDA=1,N=A->n;
115 double *AP=A->val,*ss=A->sscale,*v;
121 for (v=AP,i=0;i<N;i++){ ss[i]=*v;v+=(i+2);}
122 for (i=0;i<N;i++){ ss[i]=1.0/sqrt(fabs(ss[i])+1.0e-8); }
123 dtpuscalemat(AP,ss,N);
125 dpptrf(&UPLO, &N, AP, &INFO );
130 static int DTPUMatShiftDiagonal(
void* AA,
double shift){
131 dtpumat* A=(dtpumat*) AA;
143 #define __FUNCT__ "DTPUMatAssemble" 144 static int DTPUMatAssemble(
void*M){
146 double shift=1.0e-15;
148 info= DTPUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
149 DSDPFunctionReturn(0);
152 static int DTPUMatRowNonzeros(
void*M,
int row,
double cols[],
int *ncols,
int nrows){
156 for (i=0;i<=row;i++){
159 for (i=row+1;i<nrows;i++){
162 DSDPFunctionReturn(0);
167 #define __FUNCT__ "DTPUMatDiag" 168 static int DTPUMatDiag(
void*M,
int row,
double dd){
170 dtpumat* ABA=(dtpumat*)M;
172 ii=row*(row+1)/2+row;
174 DSDPFunctionReturn(0);
177 #define __FUNCT__ "DTPUMatDiag2" 178 static int DTPUMatDiag2(
void*M,
double diag[],
int m){
180 dtpumat* ABA=(dtpumat*)M;
182 for (row=0;row<m;row++){
183 ii=row*(row+1)/2+row;
184 ABA->val[ii] +=diag[row];
186 DSDPFunctionReturn(0);
189 static int DTPUMatAddRow(
void* AA,
int nrow,
double dd,
double row[],
int n){
190 dtpumat* A=(dtpumat*) AA;
191 ffinteger ione=1,nn,nnn;
196 daxpy(&nn,&dd,row,&ione,vv+nnn,&ione);
201 static int DTPUMatZero(
void* AA){
202 dtpumat* A=(dtpumat*) AA;
203 int mn=A->n*(A->n+1)/2;
205 memset((
void*)vv,0,mn*
sizeof(
double));
209 static int DTPUMatGetSize(
void *AA,
int *n){
210 dtpumat* A=(dtpumat*) AA;
215 static int DTPUMatDestroy(
void* AA){
217 dtpumat* A=(dtpumat*) AA;
218 if (A && A->owndata){DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
219 if (A && A->sscale) {DSDPFREE(&A->sscale,&info);DSDPCHKERR(info);}
220 if (A) {DSDPFREE(&A,&info);DSDPCHKERR(info);}
224 static int DTPUMatView(
void* AA){
225 dtpumat* M=(dtpumat*) AA;
228 for (i=0; i<M->n; i++){
229 for (j=0; j<=i; j++){
230 printf(
" %9.2e",val[kk]);
241 static struct DSDPSchurMat_Ops dsdpmmatops;
243 static int DSDPInitSchurOps(
struct DSDPSchurMat_Ops* mops){
247 mops->matrownonzeros=DTPUMatRowNonzeros;
248 mops->matscaledmultiply=DTPUMatMult;
249 mops->mataddrow=DTPUMatAddRow;
250 mops->mataddelement=DTPUMatDiag;
251 mops->matadddiagonal=DTPUMatDiag2;
252 mops->matshiftdiagonal=DTPUMatShiftDiagonal;
253 mops->matassemble=DTPUMatAssemble;
254 mops->matfactor=DTPUMatCholeskyFactor;
255 mops->matsolve=DTPUMatSolve;
256 mops->matdestroy=DTPUMatDestroy;
257 mops->matzero=DTPUMatZero;
258 mops->matview=DTPUMatView;
260 mops->matname=lapackname;
261 DSDPFunctionReturn(0);
265 #define __FUNCT__ "DSDPGetLAPACKPUSchurOps" 266 int DSDPGetLAPACKPUSchurOps(
int n,
struct DSDPSchurMat_Ops** sops,
void** mdata){
267 int info,nn=n*(n+1)/2;
271 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
272 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
275 info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info);
278 DSDPFunctionReturn(0);
282 static void daddrow(
double *v,
double alpha,
int i,
double row[],
int n){
284 ffinteger j,nn=n,ione=1;
285 nn=i+1; s1=v+i*(i+1)/2;
286 daxpy(&nn,&alpha,s1,&ione,row,&ione);
294 static int DTPUMatInverseMult(
void* AA,
int indx[],
int nind,
double x[],
double y[],
int n){
295 dtpumat* A=(dtpumat*) AA;
296 ffinteger ione=1,N=n;
297 double BETA=0.0,ALPHA=1.0;
298 double *AP=A->v2,*Y=y,*X=x;
302 if (A->n != n)
return 1;
303 if (x==0 && n>0)
return 3;
306 memset((
void*)y,0,n*
sizeof(
double));
307 for (ii=0;ii<nind;ii++){
308 i=indx[ii]; ALPHA=x[i];
309 daddrow(AP,ALPHA,i,y,n);
313 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione);
319 static int DTPUMatCholeskyBackward(
void* AA,
double b[],
double x[],
int n){
320 dtpumat* M=(dtpumat*) AA;
321 ffinteger N=M->n,INCX=1;
322 double *AP=M->val,*ss=M->sscale;
323 char UPLO=M->UPLO,TRANS=
'N',DIAG=
'N';
324 memcpy(x,b,N*
sizeof(
double));
325 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX );
326 dtpuscalevec(1.0,ss,x,x,n);
331 static int DTPUMatCholeskyForward(
void* AA,
double b[],
double x[],
int n){
332 dtpumat* M=(dtpumat*) AA;
333 ffinteger N=M->n,INCX=1;
334 double *AP=M->val,*ss=M->sscale;
335 char UPLO=M->UPLO,TRANS=
'T',DIAG=
'N';
336 dtpuscalevec(1.0,ss,b,x,n);
337 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX);
341 static int DenseSymPSDCholeskyForwardMultiply(
void* AA,
double x[],
double y[],
int n){
342 dtpumat* A=(dtpumat*) AA;
345 double *AP=A->val,*ss=A->sscale;
347 if (x==0 && N>0)
return 3;
354 for (i=0;i<N;i++){ y[i]=y[i]/ss[i];}
358 static int DTPUMatLogDet(
void* AA,
double *dd){
359 dtpumat* A=(dtpumat*) AA;
361 double d=0,*v=A->val,*ss=A->sscale;
371 static int DTPUMatInvert(
void* AA){
372 dtpumat* A=(dtpumat*) AA;
373 ffinteger INFO,N=A->n,nn=N*(N+1)/2;
374 double *v=A->val,*AP=A->v2,*ss=A->sscale;
376 memcpy((
void*)AP,(
void*)v,nn*
sizeof(
double));
377 dpptri(&UPLO, &N, AP, &INFO );
379 INFO=DTPUMatShiftDiagonal(AA,1e-11);
380 memcpy((
void*)AP,(
void*)v,nn*
sizeof(
double));
381 dpptri(&UPLO, &N, AP, &INFO );
384 dtpuscalemat(AP,ss,N);
389 static int DTPUMatInverseAdd(
void* AA,
double alpha,
double y[],
int nn,
int n){
390 dtpumat* A=(dtpumat*) AA;
391 ffinteger N=nn,ione=1;
393 daxpy(&N,&alpha,v2,&ione,y,&ione);
398 static int DTPUMatScaleDiagonal(
void* AA,
double dd){
399 dtpumat* A=(dtpumat*) AA;
409 static int DTPUMatOuterProduct(
void* AA,
double alpha,
double x[],
int n){
410 dtpumat* A=(dtpumat*) AA;
411 ffinteger ione=1,N=n;
414 dspr(&UPLO,&N,&alpha,x,&ione,v);
419 static int DenseSymPSDNormF2(
void* AA,
int n,
double *dddot){
420 dtpumat* A=(dtpumat*) AA;
421 ffinteger ione=1,nn=A->n*(A->n+1)/2;
422 double dd,tt=sqrt(0.5),*val=A->val;
424 info=DTPUMatScaleDiagonal(AA,tt);
425 dd=dnrm2(&nn,val,&ione);
426 info=DTPUMatScaleDiagonal(AA,1.0/tt);
451 static int DTPUMatFull(
void*A,
int*full){
457 static int DTPUMatGetDenseArray(
void* A,
double *v[],
int*n){
458 dtpumat* ABA=(dtpumat*)A;
460 *n=(ABA->n)*(ABA->n+1)/2;
464 static int DTPUMatRestoreDenseArray(
void* A,
double *v[],
int *n){
469 static int DDenseSetXMat(
void*A,
double v[],
int nn,
int n){
471 dtpumat* ABA=(dtpumat*)A;
474 memcpy((
void*)vv,(
void*)v,nn*
sizeof(
double));
479 static int DDenseVecVec(
void* A,
double x[],
int n,
double *v){
480 dtpumat* ABA=(dtpumat*)A;
482 double dd=0,*val=ABA->val;
486 dd+=2*x[i]*x[j]*val[k];
489 dd+=x[i]*x[i]*val[k];
497 static int DSDPDSDenseInitializeOps(
struct DSDPDSMat_Ops* densematops){
499 if (!densematops)
return 0;
501 densematops->matseturmat=DDenseSetXMat;
502 densematops->matview=DTPUMatView;
503 densematops->matdestroy=DTPUMatDestroy;
504 densematops->matgetsize=DTPUMatGetSize;
505 densematops->matzeroentries=DTPUMatZero;
506 densematops->matmult=DTPUMatMult;
507 densematops->matvecvec=DDenseVecVec;
509 densematops->matname=lapackname;
514 #define __FUNCT__ "DSDPCreateDSMatWithArray" 515 int DSDPCreateDSMatWithArray(
int n,
double vv[],
int nnz,
struct DSDPDSMat_Ops* *dsmatops,
void**dsmat){
519 info=DTPUMatCreateWData(n,vv,nnz,&AA); DSDPCHKERR(info);
521 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
522 *dsmatops=&tdsdensematops;
524 DSDPFunctionReturn(0);
529 #define __FUNCT__ "DSDPCreateDSMat" 530 int DSDPCreateDSMat(
int n,
struct DSDPDSMat_Ops* *dsmatops,
void**dsmat){
531 int info,nn=n*(n+1)/2;
535 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
536 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
537 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
538 *dsmatops=&tdsdensematops;
541 DSDPFunctionReturn(0);
546 static int DSDPDenseXInitializeOps(
struct DSDPVMat_Ops* densematops){
548 if (!densematops)
return 0;
550 densematops->matview=DTPUMatView;
551 densematops->matscalediagonal=DTPUMatScaleDiagonal;
552 densematops->matshiftdiagonal=DTPUMatShiftDiagonal;
553 densematops->mataddouterproduct=DTPUMatOuterProduct;
554 densematops->matdestroy=DTPUMatDestroy;
555 densematops->matfnorm2=DenseSymPSDNormF2;
556 densematops->matgetsize=DTPUMatGetSize;
557 densematops->matzeroentries=DTPUMatZero;
558 densematops->matgeturarray=DTPUMatGetDenseArray;
559 densematops->matrestoreurarray=DTPUMatRestoreDenseArray;
560 densematops->matmineig=DTPUMatEigs;
561 densematops->matmult=DTPUMatMult;
563 densematops->matname=lapackname;
568 #define __FUNCT__ "DSDPXMatCreate" 569 int DSDPXMatPCreate(
int n,
struct DSDPVMat_Ops* *xops,
void * *xmat){
570 int info,nn=n*(n+1)/2;
574 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
575 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
577 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
578 *xops=&turdensematops;
580 DSDPFunctionReturn(0);
584 #define __FUNCT__ "DSDPXMatCreate" 585 int DSDPXMatPCreateWithData(
int n,
double nz[],
int nnz,
struct DSDPVMat_Ops* *xops,
void * *xmat){
590 for (i=0;i<n*(n+1)/2;i++) dtmp=nz[i];
591 info=DTPUMatCreateWData(n,nz,nnz,&AA); DSDPCHKERR(info);
592 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
593 *xops=&turdensematops;
595 DSDPFunctionReturn(0);
602 if (sops==NULL)
return 0;
604 sops->matseturmat=DDenseSetXMat;
605 sops->matcholesky=DTPUMatCholeskyFactor;
606 sops->matsolveforward=DTPUMatCholeskyForward;
607 sops->matsolvebackward=DTPUMatCholeskyBackward;
608 sops->matinvert=DTPUMatInvert;
609 sops->matinverseadd=DTPUMatInverseAdd;
610 sops->matinversemultiply=DTPUMatInverseMult;
611 sops->matforwardmultiply=DenseSymPSDCholeskyForwardMultiply;
612 sops->matfull=DTPUMatFull;
613 sops->matdestroy=DTPUMatDestroy;
614 sops->matgetsize=DTPUMatGetSize;
615 sops->matview=DTPUMatView;
616 sops->matlogdet=DTPUMatLogDet;
617 sops->matname=lapackname;
624 #define __FUNCT__ "DSDPLAPACKDualMatCreate" 625 int DSDPLAPACKPUDualMatCreate(
int n,
struct DSDPDualMat_Ops **sops,
void**smat){
626 int info,nn=n*(n+1)/2;
630 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
631 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info);
634 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
637 DSDPFunctionReturn(0);
640 static int switchptr(
void* SD,
void *SP){
651 #define __FUNCT__ "DSDPLAPACKDualMatCreate2" 652 int DSDPLAPACKPUDualMatCreate2(
int n,
657 info=DSDPLAPACKPUDualMatCreate(n,sops1,smat1);DSDPCHKERR(info);
658 info=DSDPLAPACKPUDualMatCreate(n,sops2,smat2);DSDPCHKERR(info);
659 info=switchptr(*smat1,*smat2);
660 DSDPFunctionReturn(0);
677 #define __FUNCT__ "CreateDvechmatWData" 678 static int CreateDvechmatWdata(
int n,
double alpha,
double vv[], dvechmat **A){
679 int info,nn=(n*n+n)/2;
681 DSDPCALLOC1(&V,dvechmat,&info);DSDPCHKERR(info);
682 info=DTPUMatCreateWData(n,vv,nn,&V->AA); DSDPCHKERR(info);
692 static int DvechmatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int n){
695 for (k=0;k<n;k++) nz[k]++;
699 static int DTPUMatGetRowAdd(
void* AA,
int nrow,
double ytmp,
double row[],
int n){
700 dtpumat* A=(dtpumat*) AA;
705 for (i=0;i<nrow;i++){
706 row[i]+=ytmp*v[nnn+i];
708 row[nrow]+=ytmp*v[nnn+nrow];
709 for (i=nrow+1;i<n;i++){
716 static int DvechmatGetRowAdd(
void* AA,
int trow,
double scl,
double r[],
int m){
718 dvechmat* A=(dvechmat*)AA;
719 info=DTPUMatGetRowAdd((
void*)A->AA ,trow,scl*A->alpha,r,m);
723 static int DvechmatAddMultiple(
void* AA,
double alpha,
double r[],
int nnn,
int n){
724 dvechmat* A=(dvechmat*)AA;
725 ffinteger nn=nnn, ione=1;
726 double *val=A->AA->val;
728 daxpy(&nn,&alpha,val,&ione,r,&ione);
733 static int DvechEigVecVec(
void*,
double[],
int,
double*);
734 static int DvechmatVecVec(
void* AA,
double x[],
int n,
double *v){
735 dvechmat* A=(dvechmat*)AA;
737 double dd=0,*val=A->AA->val;
739 if (A->Eig.neigs<n/5){
740 i=DvechEigVecVec(AA,x,n,&dd);
745 dd+=2*x[i]*x[j]*val[k];
748 dd+=x[i]*x[i]*val[k];
757 static int DvechmatFNorm2(
void* AA,
int n,
double *v){
758 dvechmat* A=(dvechmat*)AA;
760 double dd=0,*x=A->AA->val;
769 *v=dd*A->alpha*A->alpha;
774 static int DvechmatCountNonzeros(
void* AA,
int *nnz,
int n){
780 static int DvechmatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
781 dvechmat* A=(dvechmat*)AA;
782 ffinteger ione=1,nnn=nn;
783 double dd,*val=A->AA->val;
784 dd=ddot(&nnn,val,&ione,x,&ione);
796 #define __FUNCT__ "DvechmatDestroy" 797 static int DvechmatDestroy(
void* AA){
798 dvechmat* A=(dvechmat*)AA;
800 info=DTPUMatDestroy((
void*)(A->AA));
801 if (A->Eig.an){DSDPFREE(&A->Eig.an,&info);DSDPCHKERR(info);}
802 if (A->Eig.eigval){DSDPFREE(&A->Eig.eigval,&info);DSDPCHKERR(info);}
804 DSDPFREE(&A,&info);DSDPCHKERR(info);
809 static int DvechmatView(
void* AA){
810 dvechmat* A=(dvechmat*)AA;
814 for (i=0; i<M->n; i++){
815 for (j=0; j<=i; j++){
816 printf(
" %4.2e",A->alpha*val[kk]);
826 #define __FUNCT__ "DSDPCreateDvechmatEigs" 827 static int CreateEigenLocker(Eigen *E,
int neigs,
int n){
829 DSDPCALLOC2(&E->eigval,
double,neigs,&info);DSDPCHKERR(info);
830 DSDPCALLOC2(&E->an,
double,n*neigs,&info);DSDPCHKERR(info);
836 static int EigMatSetEig(Eigen* A,
int row,
double eigv,
double v[],
int n){
839 memcpy((
void*)(an+n*row),(
void*)v,n*
sizeof(
double));
844 static int EigMatGetEig(Eigen* A,
int row,
double *eigenvalue,
double eigenvector[],
int n){
846 *eigenvalue=A->eigval[row];
847 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
852 static int DvechmatComputeEigs(dvechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int);
854 static int DvechmatFactor(
void*AA,
double dmatp[],
int nn0,
double dwork[],
int n,
double ddwork[],
int n1,
int iptr[],
int n2){
857 dvechmat* A=(dvechmat*)AA;
858 if (A->Eig.neigs>=0)
return 0;
859 info=DvechmatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2);DSDPCHKERR(info);
863 static int DvechmatGetRank(
void *AA,
int *rank,
int n){
864 dvechmat* A=(dvechmat*)AA;
865 if (A->Eig.neigs>=0){
868 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
873 static int DvechmatGetEig(
void* AA,
int rank,
double *eigenvalue,
double vv[],
int n,
int indz[],
int *nind){
874 dvechmat* A=(dvechmat*)AA;
878 info=EigMatGetEig(&A->Eig,rank,&dd,vv,n);DSDPCHKERR(info);
880 *eigenvalue=dd*A->alpha;
881 for (i=0;i<n;i++){ indz[i]=i;}
883 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
888 static int DvechEigVecVec(
void* AA,
double v[],
int n,
double *vv){
889 dvechmat* A=(dvechmat*)AA;
891 double *an,dd,ddd=0,*eigval;
892 if (A->Eig.neigs>=0){
895 eigval=A->Eig.eigval;
896 for (rank=0;rank<neigs;rank++){
897 for (dd=0,i=0;i<n;i++){
901 ddd+=dd*dd*eigval[rank];
905 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
912 static const char *datamatname=
"DENSE VECH MATRIX";
916 if (sops==NULL)
return 0;
918 sops->matvecvec=DvechmatVecVec;
919 sops->matdot=DvechmatDot;
920 sops->mataddrowmultiple=DvechmatGetRowAdd;
921 sops->mataddallmultiple=DvechmatAddMultiple;
922 sops->matview=DvechmatView;
923 sops->matdestroy=DvechmatDestroy;
924 sops->matfactor2=DvechmatFactor;
925 sops->matgetrank=DvechmatGetRank;
926 sops->matgeteig=DvechmatGetEig;
927 sops->matrownz=DvechmatGetRowNnz;
928 sops->matfnorm2=DvechmatFNorm2;
929 sops->matnnz=DvechmatCountNonzeros;
931 sops->matname=datamatname;
936 #define __FUNCT__ "DSDPGetDmat" 937 int DSDPGetDMat(
int n,
double alpha,
double *val,
struct DSDPDataMat_Ops**sops,
void**smat){
943 for (k=0;k<n*(n+1)/2;++k) dtmp=val[k];
944 info=CreateDvechmatWdata(n,alpha,val,&A); DSDPCHKERR(info);
946 info=DvechmatOpsInitialize(&dvechmatops); DSDPCHKERR(info);
947 if (sops){*sops=&dvechmatops;}
948 if (smat){*smat=(
void*)A;}
949 DSDPFunctionReturn(0);
954 #define __FUNCT__ "DvechmatComputeEigs" 955 static int DvechmatComputeEigs(dvechmat* AA,
double DD[],
int nn0,
double W[],
int n,
double WORK[],
int n1,
int iiptr[],
int n2){
957 int i,j,k,neigs,info;
958 long int *i2darray=(
long int*)DD;
959 int ownarray1=0,ownarray2=0,ownarray3=0;
960 double *val=AA->AA->val;
961 double *dmatarray=0,*dworkarray=0,eps=1.0e-12;
966 DSDPCALLOC2(&dmatarray,
double,(n*n),&info); DSDPCHKERR(info);
969 memset((
void*)dmatarray,0,n*n*
sizeof(
double));
972 DSDPCALLOC2(&dworkarray,
double,(n*n),&info); DSDPCHKERR(info);
976 if (n*n*
sizeof(
long int)>nn0*
sizeof(
double)){
977 DSDPCALLOC2(&i2darray,
long int,(n*n),&info); DSDPCHKERR(info);
982 for (k=0,i=0; i<n; i++){
983 for (j=0; j<=i; j++){
984 dmatarray[i*n+j] += val[k];
986 dmatarray[j*n+i] += val[k];
992 info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n,
993 W,n,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
996 for (neigs=0,i=0;i<n;i++){
997 if (fabs(W[i])> eps ){ neigs++;}
1000 info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info);
1003 for (neigs=0,i=0; i<n; i++){
1004 if (fabs(W[i]) > eps){
1005 info=EigMatSetEig(&AA->Eig,neigs,W[i],dmatarray+n*i,n);DSDPCHKERR(info);
1010 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
1011 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
1012 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Structure of function pointers that each dense matrix array type (upper full, packed symmetric,...
Table of function pointers that operate on the dense matrix.
Error handling, printing, and profiling.
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
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.
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.
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.