31 static int DTRUMatView(
void*);
35 #define __FUNCT__ "DSDPLAPACKROUTINE" 36 static void dtruscalevec(
double alpha,
double v1[],
double v2[],
double v3[],
int n){
39 v3[i] = (alpha*v1[i]*v2[i]);
44 static void dsydadd(
double x[],
double s[],
double y[],
int n){
47 y[i] += x[i]*(1/(s[i]*s[i])-2);
53 static 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);
61 static void DTRUMatOwnData(
void* A,
int owndata){
62 dtrumat* AA=(dtrumat*)A;
67 static int SUMatGetLDA(
int n){
69 if (n>8 && nlda%2==1){ nlda++;}
71 while (nlda%8!=0){ nlda++;}
79 static 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" 100 int 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" 150 int 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" 160 int 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" 181 static int DTRUMatMult(
void* AA,
double x[],
double y[],
int n){
182 dtrumat* A=(dtrumat*) AA;
183 ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
184 double BETA=0.0,ALPHA=1.0;
185 double *AP=A->val,*Y=y,*X=x;
186 char UPLO=A->UPLO,TRANS=
'N';
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);
199 static int DTRUMatMultR(
void* AA,
double x[],
double y[],
int n){
200 dtrumat* A=(dtrumat*) AA;
201 ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
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);
228 static void DTRUMatScale(
void*AA){
229 dtrumat* A=(dtrumat*) AA;
230 int i,N=A->n,LDA=A->LDA;
231 double *ss=A->sscale,*v=A->val;
232 for (i=0;i<N;i++){ ss[i]=*v;v+=(LDA+1);}
233 for (i=0;i<N;i++){
if (ss[i]!=0.0){ss[i]=1.0/sqrt(fabs(ss[i]));}
else {ss[i]=1.0; }}
234 dtruscalemat(A->val,ss,N,LDA);
237 static 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 );
251 int dtrsm2(
char*,
char*,
char*,
char*,ffinteger*,ffinteger*,
double*,
double*,ffinteger*,
double*,ffinteger*);
253 static int DTRUMatSolve(
void* AA,
double b[],
double x[],
int n){
254 dtrumat* A=(dtrumat*) AA;
255 ffinteger ierr,INFO=0,NRHS=1,LDA=A->LDA,LDB=A->LDA,N=A->n;
256 double *AP=A->val,*ss=A->sscale,ONE=1.0;
257 char SIDE=
'L',UPLO=A->UPLO,TRANSA=
'N',DIAG=
'N';
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);
275 static int DTRUMatShiftDiagonal(
void* AA,
double shift){
276 dtrumat* A=(dtrumat*) AA;
277 int i,n=A->n, LDA=A->LDA;
285 static int DTRUMatAddRow(
void* AA,
int nrow,
double dd,
double row[],
int n){
286 dtrumat* A=(dtrumat*) AA;
287 ffinteger ione=1,LDA=A->LDA,nn,INCX=1,INCY=A->LDA;
291 daxpy(&nn,&dd,row,&INCX,vv+nrow,&INCY);
293 daxpy(&nn,&dd,row,&ione,vv+nrow*LDA,&ione);
298 static int DTRUMatZero(
void* AA){
299 dtrumat* A=(dtrumat*) AA;
300 int mn=A->n*(A->LDA);
302 memset((
void*)vv,0,mn*
sizeof(
double));
307 static int DTRUMatGetSize(
void *AA,
int *n){
308 dtrumat* A=(dtrumat*) AA;
313 static 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);}
323 static 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]);
340 static 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" 367 static int DTRUMatAssemble(
void*M){
369 double shift=1.0e-15;
371 info= DTRUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
372 DSDPFunctionReturn(0);
375 static 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" 389 static 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" 400 static 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);
411 static struct DSDPSchurMat_Ops dsdpmmatops;
412 static const char* lapackname=
"DENSE,SYMMETRIC U STORAGE";
414 static 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" 439 int 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);
457 static int DTRUMatCholeskyBackward(
void* AA,
double b[],
double x[],
int n){
458 dtrumat* M=(dtrumat*) AA;
459 ffinteger N=M->n,INCX=1,LDA=M->LDA;
460 double *AP=M->val,*ss=M->sscale;
461 char UPLO=M->UPLO,TRANS=
'N',DIAG=
'N';
462 memcpy(x,b,N*
sizeof(
double));
463 dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
464 dtruscalevec(1.0,ss,x,x,n);
469 static int DTRUMatCholeskyForward(
void* AA,
double b[],
double x[],
int n){
470 dtrumat* M=(dtrumat*) AA;
471 ffinteger N=M->n,INCX=1,LDA=M->LDA;
472 double *AP=M->val,*ss=M->sscale;
473 char UPLO=M->UPLO,TRANS=
'T',DIAG=
'N';
474 dtruscalevec(1.0,ss,b,x,n);
475 dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
479 static int DTRUMatLogDet(
void* AA,
double *dd){
480 dtrumat* A=(dtrumat*) AA;
481 int i,n=A->n,LDA=A->LDA;
482 double d=0,*v=A->val,*ss=A->sscale;
493 static int DTRUMatCholeskyForwardMultiply(
void* AA,
double x[],
double y[],
int n){
494 dtrumat* A=(dtrumat*) AA;
495 ffinteger i,j,N=A->n,LDA=A->LDA;
496 double *AP=A->val,*ss=A->sscale;
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];}
516 static int DTRUMatCholeskyBackwardMultiply(
void* AA,
double x[],
double y[],
int n){
517 dtrumat* A=(dtrumat*) AA;
518 ffinteger i,j,N=A->n,LDA=A->LDA;
519 double *AP=A->val,*ss=A->sscale;
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];}
539 static int DTRUMatInvert(
void* AA){
540 dtrumat* A=(dtrumat*) AA;
541 ffinteger INFO,LDA=A->LDA,N=A->n,nn=LDA*N;
542 double *v=A->val,*AP=A->v2,*ss=A->sscale;
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);
559 static void DTRUSymmetrize(dtrumat* A){
560 int i,j,n=A->n,row,LDA=A->LDA;
561 double *v2=A->v2,*r1=A->v2,*r2=A->v2+LDA,*c1;
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;
591 static int DTRUMatInverseAdd(
void* AA,
double alpha,
double y[],
int nn,
int n){
592 dtrumat* A=(dtrumat*) AA;
593 ffinteger i,LDA=A->LDA,N,ione=1;
597 daxpy(&N,&alpha,v2,&ione,y,&ione);
603 static int DTRUMatInverseAddP(
void* AA,
double alpha,
double y[],
int nn,
int n){
604 dtrumat* A=(dtrumat*) AA;
605 ffinteger N,LDA=A->LDA,i,ione=1;
609 daxpy(&N,&alpha,v2,&ione,y,&ione);
615 static 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; }
629 static int DTRUMatInverseMultiply(
void* AA,
int indx[],
int nind,
double x[],
double y[],
int n){
630 dtrumat* A=(dtrumat*) AA;
631 ffinteger nn=n,LDA=A->LDA,N=A->n,INCX=1,INCY=1;
632 double *AP=A->v2,*s1=A->v2,*s2,*X=x,*Y=y,ALPHA=1.0,BETA=0.0;
633 char UPLO=A->UPLO,TRANS=
'N';
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);
666 static 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;
679 static 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;
693 static int DTRUMatFull(
void*A,
int*full){
698 static 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" 733 static 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);
750 static int switchptr(
void *SD,
void *SP){
761 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate2" 762 int 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" 800 static 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" 819 int DSDPLAPACKSUDualMatCreate2P(
int n,
824 info=DSDPLAPACKSUDualMatCreateP(n,sops1,smat1);
825 info=DSDPLAPACKSUDualMatCreateP(n,sops2,smat2);
826 info=switchptr(*smat1,*smat2);
827 DSDPFunctionReturn(0);
830 static int DTRUMatScaleDiagonal(
void* AA,
double dd){
831 dtrumat* A=(dtrumat*) AA;
832 ffinteger LDA=A->LDA;
842 static int DTRUMatOuterProduct(
void* AA,
double alpha,
double x[],
int n){
843 dtrumat* A=(dtrumat*) AA;
844 ffinteger ione=1,N=n,LDA=A->LDA;
847 dsyr(&UPLO,&N,&alpha,x,&ione,v,&LDA);
851 static int DenseSymPSDNormF2(
void* AA,
int n,
double *dddot){
852 dtrumat* A=(dtrumat*) AA;
853 ffinteger ione=1,nn=A->n*A->n;
854 double dd,tt=sqrt(0.5),*val=A->val;
856 info=DTRUMatScaleDiagonal(AA,tt);
857 dd=dnrm2(&nn,val,&ione);
858 info=DTRUMatScaleDiagonal(AA,1.0/tt);
863 static int DTRUMatGetDenseArray(
void* A,
double *v[],
int*n){
864 dtrumat* ABA=(dtrumat*)A;
870 static int DTRUMatRestoreDenseArray(
void* A,
double *v[],
int *n){
875 static 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;
889 static int DDenseVecVec(
void* A,
double x[],
int n,
double *v){
890 dtrumat* ABA=(dtrumat*)A;
891 int i,j,k=0,LDA=ABA->LDA;
892 double dd=0,*val=ABA->val;
896 dd+=2*x[i]*x[j]*val[j];
899 dd+=x[i]*x[i]*val[i];
908 static int DTRUMatEigs(
void*AA,
double W[],
double IIWORK[],
int nn1,
double *mineig){
909 dtrumat* AAA=(dtrumat*) AA;
910 ffinteger info,INFO=0,M,N=AAA->n;
911 ffinteger IL=1,IU=1,LDA=AAA->LDA,LDZ=LDA,LWORK,IFAIL;
912 ffinteger *IWORK=(ffinteger*)IIWORK;
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);
933 static 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" 956 int 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" 973 int 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);
984 static 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" 1002 int 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" 1016 int 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" 1040 static 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);
1051 static int DvecumatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int n){
1054 for (k=0;k<n;k++) nz[k]++;
1058 static 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];
1074 static 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);
1081 static int DvecumatAddMultiple(
void* AA,
double alpha,
double r[],
int nnn,
int n){
1082 dvecumat* A=(dvecumat*)AA;
1083 ffinteger nn=nnn, ione=1;
1084 double *val=A->AA->val;
1085 daxpy(&nn,&alpha,val,&ione,r,&ione);
1090 static int DvecuEigVecVec(
void*,
double[],
int,
double*);
1091 static int DvecumatVecVec(
void* AA,
double x[],
int n,
double *v){
1092 dvecumat* A=(dvecumat*)AA;
1093 int i,j,k=0,LDA=A->AA->LDA;
1094 double dd=0,*val=A->AA->val;
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];
1112 static int DvecumatFNorm2(
void* AA,
int n,
double *v){
1113 dvecumat* A=(dvecumat*)AA;
1114 long int i,j,k=0,LDA=A->AA->LDA;
1115 double dd=0,*x=A->AA->val;
1116 for (i=0; i<n; i++){
1128 static int DvecumatCountNonzeros(
void* AA,
int *nnz,
int n){
1134 static int DvecumatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
1135 dvecumat* A=(dvecumat*)AA;
1136 double d1,dd=0,*v1=x,*v2=A->AA->val;
1137 ffinteger i,n2,ione=1,LDA=A->AA->LDA;
1141 d1=ddot(&n2,v1,&ione,v2,&ione);
1156 #define __FUNCT__ "DvecumatDestroy" 1157 static 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);
1171 static 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" 1188 static 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);
1201 static 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));
1209 static 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));
1217 static int DvecumatComputeEigs(dvecumat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int);
1219 static 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);
1228 static 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");
1238 static 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");
1251 static 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");
1275 static 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" 1300 int 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" 1318 static 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 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.