27 #define GETI(a,b) (int)((int)a/(int)b) 28 #define GETJ(a,b) (int)((int)a%(int)b) 30 static void getij(
int k,
int n,
int *i,
int *j){
37 #define __FUNCT__ "CreateVechMatWData" 38 static int CreateVechMatWdata(
int n,
int ishift,
double alpha,
const int *ind,
const double *vals,
int nnz, vechmat **A){
41 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
42 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
49 static int VechMatAddMultiple(
void* AA,
double scl,
double r[],
int nn,
int n){
50 vechmat* A=(vechmat*)AA;
52 const int *ind=A->ind,nnz=A->nnzeros;
53 const double *val=A->val;
54 double *rr=r-A->ishift;
56 for (k=0; k<nnz; ++k){
57 *(rr+(ind[k])) +=scl*(val[k]);
62 static int VechMatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
63 vechmat* A=(vechmat*)AA;
65 const int *ind=A->ind;
66 double vv=0,*xx=x-A->ishift;
67 const double *val=A->val;
68 for (k=0;k<nnz;++k,++ind,++val){
69 vv+=(*val)*(*(xx+(*ind)));
75 static int EigMatVecVec(Eigen*,
double[],
int,
double*);
76 static int VechMatGetRank(
void*,
int*,
int);
78 static int VechMatVecVec(
void* AA,
double x[],
int n,
double *v){
79 vechmat* A=(vechmat*)AA;
80 int info,rank=n,i=0,j,k,kk;
81 const int *ind=A->ind,ishift=A->ishift,nnz=A->nnzeros;
83 const double *val=A->val;
86 info=VechMatGetRank(AA,&rank,n);
87 if (nnz>3 && rank<nnz){
88 info=EigMatVecVec(A->Eig,x,n,&vv);
94 for (k=0; k<nnz; ++k,++ind,++val){
108 static int VechMatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int nn){
109 vechmat* A=(vechmat*)AA;
111 const int *ind=A->ind, ishift=A->ishift,nnz=A->nnzeros;
113 for (k=0; k<nnz; ++k,ind){
126 static int VechMatFNorm2(
void* AA,
int n,
double *fnorm2){
127 vechmat* A=(vechmat*)AA;
129 const int *ind=A->ind,ishift=A->ishift,nnz=A->nnzeros;
131 const double *val=A->val;
132 for (k=0; k<nnz; ++k){
139 fn2+= 2*val[k]*val[k];
142 *fnorm2=fn2*A->alpha*A->alpha;
146 static int VechMatAddRowMultiple(
void* AA,
int trow,
double scl,
double r[],
int m){
147 vechmat* A=(vechmat*)AA;
148 int i=0,j,k,t,ishift=A->ishift,nnz=A->nnzeros;
149 const int *ind=A->ind;
150 const double *val=A->val;
152 for (k=0; k<nnz; ++k){
165 static int VechMatCountNonzeros(
void* AA,
int*nnz,
int n){
166 vechmat* A=(vechmat*)AA;
172 #define __FUNCT__ "VechMatDestroy" 173 static int VechMatDestroy(
void* AA){
174 vechmat* A=(vechmat*)AA;
184 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
185 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
186 if (A->Eig->cols){DSDPFREE(&A->Eig->cols,&info);DSDPCHKERR(info);}
187 if (A->Eig->nnz){DSDPFREE(&A->Eig->nnz,&info);DSDPCHKERR(info);}
188 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
190 DSDPFREE(&A,&info);DSDPCHKERR(info);
197 #define __FUNCT__ "DSDPCreateVechMatEigs" 198 static int CreateEigenLocker(Eigen **EE,
int iptr[],
int neigs,
int n){
202 for (k=0,i=0;i<neigs;i++) k+=iptr[i];
205 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
206 DSDPCALLOC2(&E->eigval,
double,neigs,&info);DSDPCHKERR(info);
207 DSDPCALLOC2(&E->an,
double,n*neigs,&info);DSDPCHKERR(info);
214 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
215 DSDPCALLOC2(&E->eigval,
double,neigs,&info);DSDPCHKERR(info);
216 DSDPCALLOC2(&E->nnz,
int,neigs,&info);DSDPCHKERR(info);
217 DSDPCALLOC2(&E->an,
double,k,&info);DSDPCHKERR(info);
218 DSDPCALLOC2(&E->cols,
int,k,&info);DSDPCHKERR(info);
221 if (neigs>0) E->nnz[0]=iptr[0];
222 for (i=1;i<neigs;i++){E->nnz[i]=E->nnz[i-1]+iptr[i];}
229 static int EigMatSetEig(Eigen* A,
int row,
double eigv,
int idxn[],
double v[],
int nsub,
int n){
230 int j,k,*cols=A->cols;
234 k=0;
if (row>0){ k=A->nnz[row-1];}
236 for (k=0,j=0; j<nsub; j++){
237 if (v[j]==0.0)
continue;
238 cols[k]=idxn[j]; an[k]=v[j]; k++;
242 for (j=0; j<nsub; j++){
243 if (v[j]==0.0)
continue;
251 static int EigMatGetEig(Eigen* A,
int row,
double *eigenvalue,
double eigenvector[],
int n,
int spind[],
int *nind){
252 int i,*cols=A->cols,bb,ee;
254 *eigenvalue=A->eigval[row];
257 memset((
void*)eigenvector,0,n*
sizeof(
double));
258 if (row==0){ bb=0;}
else {bb=A->nnz[row-1];} ee=A->nnz[row];
260 eigenvector[cols[i]]=an[i];
261 spind[i-bb]=cols[i]; (*nind)++;
264 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
265 for (i=0;i<n;i++)spind[i]=i;
271 static int EigMatVecVec(Eigen* A,
double v[],
int n,
double *vv){
272 int i,rank,*cols=A->cols,neigs=A->neigs,*nnz=A->nnz,bb,ee;
273 double* an=A->an,*eigval=A->eigval,dd,ddd=0;
276 for (rank=0;rank<neigs;rank++){
277 if (rank==0){ bb=0;}
else {bb=nnz[rank-1];} ee=nnz[rank];
278 for (dd=0,i=bb;i<ee;i++){
279 dd+=an[i]*v[cols[i]];
281 ddd+=dd*dd*eigval[rank];
284 for (rank=0;rank<neigs;rank++){
285 for (dd=0,i=0;i<n;i++){
289 ddd+=dd*dd*eigval[rank];
297 static int VechMatComputeEigs(vechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int,
double[],
int,
double[],
int);
299 static int VechMatFactor(
void*AA,
double dmatp[],
int nn0,
double dwork[],
int n,
double ddwork[],
int n1,
int iptr[],
int n2){
301 vechmat* A=(vechmat*)AA;
302 int i,j,k,t,info,isdiag;
303 const int *ind=A->ind,ishift=A->ishift,nonzeros=A->nnzeros;
304 double *ss1=0,*ss2=0;
306 if (A->factored)
return 0;
308 memset((
void*)iptr,0,3*n*
sizeof(
int));
310 for (isdiag=1,k=0; k<nonzeros; k++){
315 if (i!=j) {isdiag=0;iptr[j]++;}
318 if (isdiag){ A->factored=1;
return 0;}
320 for (j=0,i=0; i<n; i++){
if (iptr[i]>j) j=iptr[i]; }
321 if (j<2){ A->factored=2;
return 0; }
323 info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info);
328 static int VechMatGetRank(
void *AA,
int *rank,
int n){
329 vechmat* A=(vechmat*)AA;
330 switch (A->factored){
341 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
346 static int VechMatGetEig(
void* AA,
int rank,
double *eigenvalue,
double vv[],
int n,
int indx[],
int *nind){
347 vechmat* A=(vechmat*)AA;
348 const double *val=A->val,tt=sqrt(0.5);
350 const int *ind=A->ind,ishift=A->ishift;
353 switch (A->factored){
355 memset(vv,0,n*
sizeof(
double));
360 *eigenvalue=val[rank]*A->alpha;
365 memset(vv,0,n*
sizeof(
double));
367 getij(ind[k]-ishift,n,&i,&j);
370 vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
378 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
380 indx[0]=i; indx[1]=j;
382 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
384 indx[0]=i; indx[1]=j;
389 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
390 *eigenvalue=*eigenvalue*A->alpha;
393 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
399 static int VechMatView(
void* AA){
400 vechmat* A=(vechmat*)AA;
401 int info,i=0,j,k,rank=0,ishift=A->ishift,n=A->n,nnz=A->nnzeros;
402 const int *ind=A->ind;
403 const double *val=A->val;
404 for (k=0; k<nnz; k++){
405 getij(ind[k]-ishift,n,&i,&j);
406 printf(
"Row: %d, Column: %d, Value: %10.8f \n",i,j,A->alpha*val[k]);
409 info=VechMatGetRank(AA,&rank,n);DSDPCHKERR(info);
410 printf(
"Detected Rank: %d\n",rank);
417 static const char *datamatname=
"STANDARD VECH MATRIX";
421 if (sops==NULL)
return 0;
423 sops->matvecvec=VechMatVecVec;
424 sops->matdot=VechMatDot;
425 sops->matfnorm2=VechMatFNorm2;
426 sops->mataddrowmultiple=VechMatAddRowMultiple;
427 sops->mataddallmultiple=VechMatAddMultiple;
428 sops->matview=VechMatView;
429 sops->matdestroy=VechMatDestroy;
430 sops->matfactor2=VechMatFactor;
431 sops->matgetrank=VechMatGetRank;
432 sops->matgeteig=VechMatGetEig;
433 sops->matrownz=VechMatGetRowNnz;
434 sops->matnnz=VechMatCountNonzeros;
436 sops->matname=datamatname;
453 #define __FUNCT__ "DSDPGetVecUMat" 455 int info,i,j,k,itmp,nn=n*n;
466 DSDPSETERR3(2,
"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
468 DSDPSETERR1(2,
"Illegal index value: %d. Must be >= 0\n",itmp);
471 for (k=0;k<nnz;++k) dtmp=val[k];
472 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
475 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
476 if (sops){*sops=&vechmatops;}
477 if (smat){*smat=(
void*)AA;}
478 DSDPFunctionReturn(0);
483 #define __FUNCT__ "VechMatComputeEigs" 484 static int VechMatComputeEigs(vechmat* AA,
double DD[],
int nn0,
double W[],
int n,
double WORK[],
int n1,
int iiptr[],
int n2,
double ss1[],
int nn1,
double ss2[],
int nn2){
486 int i,j,k,nsub,neigs,info,*iptr,*perm,*invp;
487 long int *i2darray=(
long int*)DD;
488 int ownarray1=0,ownarray2=0,ownarray3=0;
489 int ishift=AA->ishift,nonzeros=AA->nnzeros;
490 const int *ind=AA->ind;
491 const double *val=AA->val;
492 double *dmatarray=ss1,*dworkarray=ss2,maxeig,eps=1.0e-12,eps2=1.0e-12;
494 iptr=iiptr; perm=iptr+n; invp=perm+n;
506 for (nsub=0,i=0; i<n; i++){
507 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
512 DSDPCALLOC2(&dmatarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
515 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
517 DSDPCALLOC2(&dworkarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
521 if (nsub*nsub*
sizeof(
long int)>nn0*
sizeof(
double)){
522 DSDPCALLOC2(&i2darray,
long int,(nsub*nsub),&info); DSDPCHKERR(info);
527 for (i=0,k=0; k<nonzeros; k++){
528 getij(ind[k]-ishift,n,&i,&j);
529 dmatarray[perm[i]*nsub+perm[j]] += val[k];
531 dmatarray[perm[j]*nsub+perm[i]] += val[k];
535 memset((
void*)W,0,n*
sizeof(
double));
537 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
538 W,nsub,WORK,n1,iiptr+3*n,n2-3*n);
540 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
541 for (i=0,k=0; k<nonzeros; k++){
542 getij(ind[k]-ishift,n,&i,&j);
543 dmatarray[perm[i]*nsub+perm[j]] += val[k];
545 dmatarray[perm[j]*nsub+perm[i]] += val[k];
548 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
549 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
553 for (maxeig=0,i=0;i<nsub;i++){
554 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
556 memset((
void*)iptr,0,nsub*
sizeof(
int));
559 for (neigs=0,k=0; k<nsub; k++){
560 if (fabs(W[k]) > eps){
561 for (j=0;j<nsub;j++){
562 if (fabs(dmatarray[nsub*k+j]) >= eps2){iptr[neigs]++;
563 }
else { dmatarray[nsub*k+j]=0.0;}
573 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
574 DSDPLogInfo(0,49,
" Data Mat has %d eigenvectors.",neigs);
576 for (neigs=0,i=0; i<nsub; i++){
577 if (fabs(W[i]) > eps){
578 info=EigMatSetEig(AA->Eig,neigs,W[i],invp,dmatarray+nsub*i,nsub,n);DSDPCHKERR(info);
583 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
584 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
585 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
Error handling, printing, and profiling.
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,...
Table of function pointers that operate on the data matrix.
int DSDPGetVecUMat(int n, int ishift, double alpha, const int ind[], const double val[], int nnz, struct DSDPDataMat_Ops **sops, void **smat)
Given data in full symmetric format, create a sparse matrix usuable by DSDP.