20static const char* lapackname=
"DENSE,SYMMETRIC,PACKED STORAGE";
22int 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"
50static void dtpuscalevec(
double alpha,
double v1[],
double v2[],
double v3[],
int n){
53 v3[i] = (alpha*v1[i]*v2[i]);
57static 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);
64static 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"
85static 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);
98static 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);
111static 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 );
130static int DTPUMatShiftDiagonal(
void* AA,
double shift){
131 dtpumat* A=(dtpumat*) AA;
143#define __FUNCT__ "DTPUMatAssemble"
144static int DTPUMatAssemble(
void*M){
146 double shift=1.0e-15;
148 info= DTPUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
149 DSDPFunctionReturn(0);
152static 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"
168static 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"
178static 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);
189static 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);
201static 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));
209static int DTPUMatGetSize(
void *AA,
int *n){
210 dtpumat* A=(dtpumat*) AA;
215static 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);}
224static 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]);
241static struct DSDPSchurMat_Ops dsdpmmatops;
243static 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"
266int 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);
282static 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);
294static 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);
319static 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);
331static 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);
341static 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];}
358static int DTPUMatLogDet(
void* AA,
double *dd){
359 dtpumat* A=(dtpumat*) AA;
361 double d=0,*v=A->val,*ss=A->sscale;
371static 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);
389static 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);
398static int DTPUMatScaleDiagonal(
void* AA,
double dd){
399 dtpumat* A=(dtpumat*) AA;
409static 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);
419static 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);
451static int DTPUMatFull(
void*A,
int*full){
457static int DTPUMatGetDenseArray(
void* A,
double *v[],
int*n){
458 dtpumat* ABA=(dtpumat*)A;
460 *n=(ABA->n)*(ABA->n+1)/2;
464static int DTPUMatRestoreDenseArray(
void* A,
double *v[],
int *n){
469static int DDenseSetXMat(
void*A,
double v[],
int nn,
int n){
471 dtpumat* ABA=(dtpumat*)A;
474 memcpy((
void*)vv,(
void*)v,nn*
sizeof(
double));
479static 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];
497static 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"
515int 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"
530int 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);
546static 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"
569int 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"
585int 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"
625int 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);
640static int switchptr(
void* SD,
void *SP){
651#define __FUNCT__ "DSDPLAPACKDualMatCreate2"
652int 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"
678static 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);
692static int DvechmatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int n){
695 for (k=0;k<n;k++) nz[k]++;
699static 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++){
716static 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);
723static 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);
733static int DvechEigVecVec(
void*,
double[],
int,
double*);
734static 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];
757static 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;
774static int DvechmatCountNonzeros(
void* AA,
int *nnz,
int n){
780static 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"
797static 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);
809static 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"
827static 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);
836static int EigMatSetEig(Eigen* A,
int row,
double eigv,
double v[],
int n){
839 memcpy((
void*)(an+n*row),(
void*)v,n*
sizeof(
double));
844static 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));
852static int DvechmatComputeEigs(dvechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int);
854static 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);
863static 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");
873static 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");
888static 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");
912static 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"
937int 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"
955static 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 DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
int DSDPDSMatOpsInitialize(struct DSDPDSMat_Ops *aops)
Set pointers to null.
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.
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.