31static int DTRUMatView(
void*);
35#define __FUNCT__ "DSDPLAPACKROUTINE"
36static void dtruscalevec(
double alpha,
double v1[],
double v2[],
double v3[],
int n){
39 v3[i] = (alpha*v1[i]*v2[i]);
44static void dsydadd(
double x[],
double s[],
double y[],
int n){
47 y[i] += x[i]*(1/(s[i]*s[i])-2);
53static void dtruscalemat(
double vv[],
double ss[],
int n,
int LDA){
55 for (i=0;i<n;i++,vv+=LDA){
56 dtruscalevec(ss[i],vv,ss,vv,i+1);
61static void DTRUMatOwnData(
void* A,
int owndata){
62 dtrumat* AA=(dtrumat*)A;
67static int SUMatGetLDA(
int n){
69 if (n>8 && nlda%2==1){ nlda++;}
71 while (nlda%8!=0){ nlda++;}
79static int DTRUMatCreateWData(
int n,
int LDA,
double nz[],
int nnz, dtrumat**M){
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;
88 for (i=0;i<n;i++)M23->sscale[i]=1.0;
91 if (n<=0){M23->LDA=1;}
99#define __FUNCT__ "DSDPGetEigs"
100int DSDPGetEigs(
double A[],
int n,
double AA[],
int nn0,
long int IA[],
int nn1,
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';
109 if ( n2/2.5 > n || (ni<10*n+1) || (nd<26*n+1) || (nn0 < n*LDA) || (nn1<LDZ*n) ){
114 dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
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;
124 DTRUMatCreateWData(n,n,A,n*n,&M);
125 DTRUMatView((
void*)M);
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];}
149#define __FUNCT__ "DSDPGetEigsSTEP"
150int DSDPGetEigsSTEP(
double A[],
int n,
double AA[],
int nn0,
long int IA[],
int nn1,
152 double WORK[],
int nd,
int LLWORK[],
int ni){
154 info=DSDPGetEigs(A,n,AA,nn0,IA,nn1,W,n2,WORK,nd,LLWORK,ni);
159#define __FUNCT__ "DSDPGetEigs2"
160int DSDPGetEigs2(
double A[],
int n,
double AA[],
int nn0,
long int IA[],
int nn1,
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';
170 DTRUMatCreateWData(n,n,A,n*n,&M);
171 DTRUMatView((
void*)M);
173 dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
179#define __FUNCT__ "DSDP FULL SYMMETRIC U LAPACK ROUTINE"
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';
188 if (A->n != n)
return 1;
189 if (x==0 && n>0)
return 3;
191 dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
193 dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
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;
203 double *AP=A->val,*Y=y,*X=x,*ss=A->sscale,*W=A->workn;
204 char UPLO=A->UPLO,TRANS=
'N',DIAG=
'U';
207 if (A->n != n)
return 1;
208 if (x==0 && n>0)
return 3;
211 memset(y,0,n*
sizeof(
double));
213 memcpy(W,X,n*
sizeof(
double));
215 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
216 daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
218 memcpy(W,X,n*
sizeof(
double));
220 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
221 daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
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);
237static int DTRUMatCholeskyFactor(
void* AA,
int *flag){
238 dtrumat* A=(dtrumat*) AA;
239 ffinteger INFO,LDA=A->LDA,N=A->n;
243 if (A->scaleit){ DTRUMatScale(AA);}
244 dpotrf(&UPLO, &N, AP, &LDA, &INFO );
251int dtrsm2(
char*,
char*,
char*,
char*,ffinteger*,ffinteger*,
double*,
double*,ffinteger*,
double*,ffinteger*);
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';
259 dtruscalevec(1.0,ss,b,x,n);
262 dpotrs(&UPLO, &N, &NRHS, AP, &LDA, x, &LDB, &INFO );
265 ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
267 ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
270 dtruscalevec(1.0,ss,x,x,n);
275static int DTRUMatShiftDiagonal(
void* AA,
double shift){
276 dtrumat* A=(dtrumat*) AA;
277 int i,n=A->n, LDA=A->LDA;
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;
291 daxpy(&nn,&dd,row,&INCX,vv+nrow,&INCY);
293 daxpy(&nn,&dd,row,&ione,vv+nrow*LDA,&ione);
298static int DTRUMatZero(
void* AA){
299 dtrumat* A=(dtrumat*) AA;
300 int mn=A->n*(A->LDA);
302 memset((
void*)vv,0,mn*
sizeof(
double));
307static int DTRUMatGetSize(
void *AA,
int *n){
308 dtrumat* A=(dtrumat*) AA;
313static int DTRUMatDestroy(
void* AA){
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);}
323static int DTRUMatView(
void* AA){
324 dtrumat* M=(dtrumat*) AA;
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]);
332 for (j=i+1; j<M->LDA; j++){
333 printf(
" %9.1e",val[i*LDA+j]);
340static int DTRUMatView2(
void* AA){
341 dtrumat* M=(dtrumat*) AA;
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]);
349 for (j=i+1; j<M->LDA; j++){
350 printf(
" %9.2e",val[i*LDA+j]);
366#define __FUNCT__ "Tassemble"
367static int DTRUMatAssemble(
void*M){
369 double shift=1.0e-15;
371 info= DTRUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
372 DSDPFunctionReturn(0);
375static int DTRUMatRowNonzeros(
void*M,
int row,
double cols[],
int *ncols,
int nrows){
379 for (i=0;i<=row;i++){
382 memset((
void*)(cols+row+1),0,(nrows-row-1)*
sizeof(
int));
383 DSDPFunctionReturn(0);
388#define __FUNCT__ "TAddDiag"
389static int DTRUMatAddDiag(
void*M,
int row,
double dd){
391 dtrumat* ABA=(dtrumat*)M;
392 ffinteger LDA=ABA->LDA;
396 DSDPFunctionReturn(0);
399#define __FUNCT__ "TAddDiag2"
400static int DTRUMatAddDiag2(
void*M,
double diag[],
int m){
402 dtrumat* ABA=(dtrumat*)M;
403 ffinteger LDA=ABA->LDA;
405 for (row=0;row<m;row++){
407 ABA->val[ii] +=diag[row];
409 DSDPFunctionReturn(0);
411static struct DSDPSchurMat_Ops dsdpmmatops;
412static const char* lapackname=
"DENSE,SYMMETRIC U STORAGE";
414static int DSDPInitSchurOps(
struct DSDPSchurMat_Ops* mops){
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;
432 mops->matname=lapackname;
433 DSDPFunctionReturn(0);
438#define __FUNCT__ "DSDPGetLAPACKSUSchurOps"
439int DSDPGetLAPACKSUSchurOps(
int n,
struct DSDPSchurMat_Ops** sops,
void** mdata){
447 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
448 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
450 info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info);
453 DSDPFunctionReturn(0);
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);
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 );
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;
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;
499 if (x==0 && N>0)
return 3;
504 for (i=0;i<N;i++)y[i]=0;
512 for (i=0;i<N;i++){ y[i]=y[i]/ss[i];}
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;
522 if (x==0 && N>0)
return 3;
527 for (i=0;i<N;i++)y[i]=0;
530 y[j]+=AP[j]*x[i]/ss[i];
535 for (i=0;i<-N;i++){ y[i]=y[i]/ss[i];}
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;
544 memcpy((
void*)AP,(
void*)v,nn*
sizeof(
double));
545 dpotri(&UPLO, &N, AP, &LDA, &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 );
552 dtruscalemat(AP,ss,N,LDA);
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;
569 for (j=row+2;j<n;j++){
578 for (row=2*n/2;row<n;row++){
581 for (j=row+1;j<n;j++){
587 A->status=ISymmetric;
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;
597 daxpy(&N,&alpha,v2,&ione,y,&ione);
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;
609 daxpy(&N,&alpha,v2,&ione,y,&ione);
615static void daddrow(
double *v,
double alpha,
int i,
double row[],
int n){
617 ffinteger j,nn=n,ione=1;
618 if (alpha==0){
return;}
620 daxpy(&nn,&alpha,s1,&ione,row,&ione);
622 for (j=i+1;j<n;j++){ row[j]+=alpha*s1[i]; s1+=n; }
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';
637 if (A->status==Inverted){
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);
648 dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
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);
660 dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
666static int DTRUMatSetXMat(
void*A,
double v[],
int nn,
int n){
667 dtrumat* ABA=(dtrumat*)A;
672 memcpy((
void*)vv,(
void*)v,(i+1)*
sizeof(
double));
676 ABA->status=Assemble;
679static int DTRUMatSetXMatP(
void*A,
double v[],
int nn,
int n){
680 dtrumat* ABA=(dtrumat*)A;
685 memcpy((
void*)vv,(
void*)v,(i+1)*
sizeof(
double));
689 ABA->status=Assemble;
693static int DTRUMatFull(
void*A,
int*full){
698static int DTRUMatGetArray(
void*A,
double **v,
int *n){
699 dtrumat* ABA=(dtrumat*)A;
708 if (sops==NULL)
return 0;
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;
732#define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
733static int DSDPLAPACKSUDualMatCreate(
int n,
struct DSDPDualMat_Ops **sops,
void**smat){
740 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
741 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
743 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
746 DSDPFunctionReturn(0);
750static int switchptr(
void *SD,
void *SP){
761#define __FUNCT__ "DSDPLAPACKSUDualMatCreate2"
762int DSDPLAPACKSUDualMatCreate2(
int n,
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);
776 if (sops==NULL)
return 0;
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;
799#define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
800static int DSDPLAPACKSUDualMatCreateP(
int n,
struct DSDPDualMat_Ops **sops,
void**smat){
807 DSDPCALLOC2(&vv,
double,nn,&info);DSDPCHKERR(info);
808 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
810 info=SDualOpsInitializeP(&sdmatopsp);DSDPCHKERR(info);
813 DSDPFunctionReturn(0);
818#define __FUNCT__ "DSDPLAPACKSUDualMatCreate2P"
819int DSDPLAPACKSUDualMatCreate2P(
int n,
824 info=DSDPLAPACKSUDualMatCreateP(n,sops1,smat1);
825 info=DSDPLAPACKSUDualMatCreateP(n,sops2,smat2);
826 info=switchptr(*smat1,*smat2);
827 DSDPFunctionReturn(0);
830static int DTRUMatScaleDiagonal(
void* AA,
double dd){
831 dtrumat* A=(dtrumat*) AA;
832 ffinteger LDA=A->LDA;
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;
847 dsyr(&UPLO,&N,&alpha,x,&ione,v,&LDA);
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;
856 info=DTRUMatScaleDiagonal(AA,tt);
857 dd=dnrm2(&nn,val,&ione);
858 info=DTRUMatScaleDiagonal(AA,1.0/tt);
863static int DTRUMatGetDenseArray(
void* A,
double *v[],
int*n){
864 dtrumat* ABA=(dtrumat*)A;
870static int DTRUMatRestoreDenseArray(
void* A,
double *v[],
int *n){
875static int DDenseSetXMat(
void*A,
double v[],
int nn,
int n){
876 dtrumat* ABA=(dtrumat*)A;
881 memcpy((
void*)vv,(
void*)v,nn*
sizeof(
double));
885 ABA->status=Assemble;
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;
896 dd+=2*x[i]*x[j]*val[j];
899 dd+=x[i]*x[i]*val[i];
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;
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);
919 dsyevx(&JOBZ,&RANGE,&UPLO,&N,AP,&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,&LWORK,IWORK,&IFAIL,&INFO);
924 DSDPFREE(&WORK,&info);
925 DSDPFREE(&IWORK,&info);
933static int DSDPDenseXInitializeOps(
struct DSDPVMat_Ops* densematops){
935 if (!densematops)
return 0;
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;
950 densematops->matname=lapackname;
955#define __FUNCT__ "DSDPXMatUCreateWithData"
956int DSDPXMatUCreateWithData(
int n,
double nz[],
int nnz,
struct DSDPVMat_Ops* *xops,
void * *xmat){
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);
965 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
966 *xops=&turdensematops;
968 DSDPFunctionReturn(0);
972#define __FUNCT__ "DSDPXMatUCreate"
973int DSDPXMatUCreate(
int n,
struct DSDPVMat_Ops* *xops,
void * *xmat){
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);
984static int DSDPDSDenseInitializeOps(
struct DSDPDSMat_Ops* densematops){
986 if (!densematops)
return 0;
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;
996 densematops->matname=lapackname;
1001#define __FUNCT__ "DSDPCreateDSMatWithArray2"
1002int DSDPCreateDSMatWithArray2(
int n,
double vv[],
int nnz,
struct DSDPDSMat_Ops* *dsmatops,
void**dsmat){
1006 info=DTRUMatCreateWData(n,n,vv,nnz,&AA); DSDPCHKERR(info);
1008 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
1009 *dsmatops=&tdsdensematops;
1011 DSDPFunctionReturn(0);
1015#define __FUNCT__ "DSDPCreateXDSMat2"
1016int DSDPCreateXDSMat2(
int n,
struct DSDPDSMat_Ops* *dsmatops,
void**dsmat){
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);
1039#define __FUNCT__ "CreateDvecumatWData"
1040static int CreateDvecumatWdata(
int n,
double vv[], dvecumat **A){
1043 DSDPCALLOC1(&V,dvecumat,&info);DSDPCHKERR(info);
1044 info=DTRUMatCreateWData(n,n,vv,nnz,&V->AA); DSDPCHKERR(info);
1051static int DvecumatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int n){
1054 for (k=0;k<n;k++) nz[k]++;
1058static int DTRUMatGetRowAdd(
void* AA,
int nrow,
double ytmp,
double row[],
int n){
1059 dtrumat* A=(dtrumat*) AA;
1064 for (i=0;i<=nrow;i++){
1065 row[i]+=ytmp*v[nnn+i];
1067 for (i=nrow+1;i<n;i++){
1069 row[i]+=ytmp*v[nrow];
1074static int DvecumatGetRowAdd(
void* AA,
int trow,
double scl,
double r[],
int m){
1076 dvecumat* A=(dvecumat*)AA;
1077 info=DTRUMatGetRowAdd((
void*)A->AA ,trow,scl,r,m);
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);
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;
1096 if (A->Eig && A->Eig->neigs<n/5){
1097 i=DvecuEigVecVec(AA,x,n,v);
1099 for (i=0; i<n; i++){
1101 dd+=2*x[i]*x[j]*val[j];
1103 dd+=x[i]*x[i]*val[i];
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++){
1128static int DvecumatCountNonzeros(
void* AA,
int *nnz,
int n){
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;
1141 d1=ddot(&n2,v1,&ione,v2,&ione);
1156#define __FUNCT__ "DvecumatDestroy"
1157static int DvecumatDestroy(
void* AA){
1158 dvecumat* A=(dvecumat*)AA;
1160 info=DTRUMatDestroy((
void*)(A->AA));
1162 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
1163 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
1165 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
1166 DSDPFREE(&A,&info);DSDPCHKERR(info);
1171static int DvecumatView(
void* AA){
1172 dvecumat* A=(dvecumat*)AA;
1176 for (i=0; i<M->n; i++){
1177 for (j=0; j<M->n; j++){
1178 printf(
" %4.2e",val[j]);
1187#define __FUNCT__ "DSDPCreateDvecumatEigs"
1188static int CreateEigenLocker(Eigen **EE,
int neigs,
int n){
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);
1201static int EigMatSetEig(Eigen* A,
int row,
double eigv,
double v[],
int n){
1203 A->eigval[row]=eigv;
1204 memcpy((
void*)(an+n*row),(
void*)v,n*
sizeof(
double));
1209static int EigMatGetEig(Eigen* A,
int row,
double *eigenvalue,
double eigenvector[],
int n){
1211 *eigenvalue=A->eigval[row];
1212 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
1217static int DvecumatComputeEigs(dvecumat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int);
1219static int DvecumatFactor(
void*AA,
double dmatp[],
int nn0,
double dwork[],
int n,
double ddwork[],
int n1,
int iptr[],
int n2){
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);
1228static int DvecumatGetRank(
void *AA,
int *rank,
int n){
1229 dvecumat* A=(dvecumat*)AA;
1231 *rank=A->Eig->neigs;
1233 DSDPSETERR(1,
"Vecu Matrix not factored yet\n");
1238static int DvecumatGetEig(
void* AA,
int rank,
double *eigenvalue,
double vv[],
int n,
int indz[],
int *nind){
1239 dvecumat* A=(dvecumat*)AA;
1242 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n);DSDPCHKERR(info);
1244 for (i=0;i<n;i++){ indz[i]=i;}
1246 DSDPSETERR(1,
"Vecu Matrix not factored yet\n");
1251static int DvecuEigVecVec(
void* AA,
double v[],
int n,
double *vv){
1252 dvecumat* A=(dvecumat*)AA;
1254 double *an,dd,ddd=0,*eigval;
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++){
1264 ddd+=dd*dd*eigval[rank];
1268 DSDPSETERR(1,
"Vecu Matrix not factored yet\n");
1275static const char *datamatname=
"STANDARD VECU MATRIX";
1279 if (sops==NULL)
return 0;
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;
1294 sops->matname=datamatname;
1299#define __FUNCT__ "DSDPGetDUmat"
1300int DSDPGetDUMat(
int n,
double *val,
struct DSDPDataMat_Ops**sops,
void**smat){
1306 for (k=0;k<n*n;++k) dtmp=val[k];
1307 info=CreateDvecumatWdata(n,val,&A); DSDPCHKERR(info);
1309 info=DvecumatOpsInitialize(&dvecumatops); DSDPCHKERR(info);
1310 if (sops){*sops=&dvecumatops;}
1311 if (smat){*smat=(
void*)A;}
1312 DSDPFunctionReturn(0);
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){
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;
1329 DSDPCALLOC2(&dmatarray,
double,(n*n),&info); DSDPCHKERR(info);
1332 memcpy((
void*)dmatarray,(
void*)val,n*n*
sizeof(
double));
1335 DSDPCALLOC2(&dworkarray,
double,(n*n),&info); DSDPCHKERR(info);
1339 if (n*n*
sizeof(
long int)>nn0*
sizeof(
double)){
1340 DSDPCALLOC2(&i2darray,
long int,(n*n),&info); DSDPCHKERR(info);
1346 info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n,
1347 W,n,WORK,n1,iiptr,n2);
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);
1355 for (neigs=0,i=0;i<n;i++){
1356 if (fabs(W[i])> eps ){ neigs++;}
1359 info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info);
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);
1369 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
1370 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
1371 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.