32 #define GETI(a) (int)(sqrt(2*(a)+0.25)-0.5) 33 #define GETJ(b,c) ((b)-((c)*((c)+1))/2) 35 static void getij(
int k,
int *i,
int *j){
51 #define __FUNCT__ "CreateVechMatWData" 52 static int CreateVechMatWdata(
int n,
int ishift,
double alpha,
const int *ind,
const double *vals,
int nnz, vechmat **A){
55 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
56 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
65 static int VechMatAddMultiple(
void* AA,
double scl,
double r[],
int nn,
int n){
66 vechmat* A=(vechmat*)AA;
68 const int* ind=A->ind;
69 const double *val=A->val;
70 double *rr=r-A->ishift;
72 for (k=0; k<nnz; ++k,++ind,++val){
73 *(rr+((*ind))) +=scl*(*val);
78 static int VechMatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
79 vechmat* A=(vechmat*)AA;
81 const int *ind=A->ind;
82 double vv=0, *xx=x-A->ishift;
83 const double *val=A->val;
84 for (k=0;k<nnz;++k,++ind,++val){
85 vv+=(*val)*(*(xx+(*ind)));
91 static int EigMatVecVec(Eigen*,
double[],
int,
double*);
92 static int VechMatGetRank(
void*,
int*,
int);
94 static int VechMatVecVec(
void* AA,
double x[],
int n,
double *v){
95 vechmat* A=(vechmat*)AA;
96 int info,rank=n,i=0,j,k,kk;
97 const int *ind=A->ind,ishift=A->ishift;
99 const double *val=A->val,nnz=A->nnzeros;
102 info=VechMatGetRank(AA,&rank,n);
103 if (nnz>3 && rank<nnz){
104 info=EigMatVecVec(A->Eig,x,n,&vv);
110 for (k=0; k<nnz; ++k,++ind,++val){
124 static int VechMatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int nn){
125 vechmat* A=(vechmat*)AA;
126 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
127 const int *ind=A->ind;
129 for (k=0; k<nnz; ++k,++ind){
130 getij((*ind)-ishift,&i,&j);
140 static int VechMatFNorm2(
void* AA,
int n,
double *fnorm2){
141 vechmat* A=(vechmat*)AA;
142 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
143 const int *ind=A->ind;
145 const double *val=A->val;
146 for (k=0; k<nnz; ++k,++ind){
147 getij((*ind)-ishift,&i,&j);
151 fn2+= 2*val[k]*val[k];
154 *fnorm2=fn2*A->alpha*A->alpha;
158 static int VechMatAddRowMultiple(
void* AA,
int trow,
double scl,
double r[],
int m){
159 vechmat* A=(vechmat*)AA;
160 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
161 const int *ind=A->ind;
162 const double *val=A->val;
164 for (k=0; k<nnz; ++k,++ind){
165 getij((*ind)-ishift,&i,&j);
176 static int VechMatCountNonzeros(
void* AA,
int*nnz,
int n){
177 vechmat* A=(vechmat*)AA;
183 #define __FUNCT__ "VechMatDestroy" 184 static int VechMatDestroy(
void* AA){
185 vechmat* A=(vechmat*)AA;
195 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
196 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
197 if (A->Eig->cols){DSDPFREE(&A->Eig->cols,&info);DSDPCHKERR(info);}
198 if (A->Eig->nnz){DSDPFREE(&A->Eig->nnz,&info);DSDPCHKERR(info);}
199 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
201 DSDPFREE(&A,&info);DSDPCHKERR(info);
208 #define __FUNCT__ "DSDPCreateVechMatEigs" 209 static int CreateEigenLocker(Eigen **EE,
int iptr[],
int neigs,
int n){
213 for (k=0,i=0;i<neigs;i++) k+=iptr[i];
216 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
217 DSDPCALLOC2(&E->eigval,
double,neigs,&info);DSDPCHKERR(info);
218 DSDPCALLOC2(&E->an,
double,n*neigs,&info);DSDPCHKERR(info);
225 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
226 DSDPCALLOC2(&E->eigval,
double,neigs,&info);DSDPCHKERR(info);
227 DSDPCALLOC2(&E->nnz,
int,neigs,&info);DSDPCHKERR(info);
228 DSDPCALLOC2(&E->an,
double,k,&info);DSDPCHKERR(info);
229 DSDPCALLOC2(&E->cols,
int,k,&info);DSDPCHKERR(info);
232 if (neigs>0) E->nnz[0]=iptr[0];
233 for (i=1;i<neigs;i++){E->nnz[i]=E->nnz[i-1]+iptr[i];}
240 static int EigMatSetEig(Eigen* A,
int row,
double eigv,
int idxn[],
double v[],
int nsub,
int n){
241 int j,k,*cols=A->cols;
246 k=0;
if (row>0){ k=A->nnz[row-1];}
248 for (k=0,j=0; j<nsub; j++){
249 if (v[j]==0.0)
continue;
250 cols[k]=idxn[j]; an[k]=v[j]; k++;
254 for (j=0; j<nsub; j++){
255 if (v[j]==0.0)
continue;
263 static int EigMatGetEig(Eigen* A,
int row,
double *eigenvalue,
double eigenvector[],
int n,
int spind[],
int *nind){
264 int i,*cols=A->cols,bb,ee;
266 *eigenvalue=A->eigval[row];
269 memset((
void*)eigenvector,0,n*
sizeof(
double));
270 if (row==0){ bb=0;}
else {bb=A->nnz[row-1];} ee=A->nnz[row];
272 eigenvector[cols[i]]=an[i];
273 spind[i-bb]=cols[i]; (*nind)++;
276 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
277 for (i=0;i<n;i++)spind[i]=i;
283 static int EigMatVecVec(Eigen* A,
double v[],
int n,
double *vv){
284 int i,rank,*cols=A->cols,neigs=A->neigs,*nnz=A->nnz,bb,ee;
285 double* an=A->an,*eigval=A->eigval,dd,ddd=0;
288 for (rank=0;rank<neigs;rank++){
289 if (rank==0){ bb=0;}
else {bb=nnz[rank-1];} ee=nnz[rank];
290 for (dd=0,i=bb;i<ee;i++){
291 dd+=an[i]*v[cols[i]];
293 ddd+=dd*dd*eigval[rank];
296 for (rank=0;rank<neigs;rank++){
297 for (dd=0,i=0;i<n;i++){
301 ddd+=dd*dd*eigval[rank];
309 static int VechMatComputeEigs(vechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int,
double[],
int,
double[],
int);
311 static int VechMatFactor(
void*AA,
double dmatp[],
int nn0,
double dwork[],
int n,
double ddwork[],
int n1,
int iptr[],
int n2){
312 vechmat* A=(vechmat*)AA;
313 int i,j,k,ishift=A->ishift,isdiag,nonzeros=A->nnzeros,info;
315 double *ss1=0,*ss2=0;
316 const int *ind=A->ind;
318 if (A->factored)
return 0;
319 memset((
void*)iptr,0,3*n*
sizeof(
int));
321 for (isdiag=1,k=0; k<nonzeros; k++){
322 getij(ind[k]-ishift,&i,&j);
324 if (i!=j) {isdiag=0;iptr[j]++;}
327 if (isdiag){ A->factored=1;
return 0;}
329 for (j=0,i=0; i<n; i++){
if (iptr[i]>j) j=iptr[i]; }
330 if (j<2){ A->factored=2;
return 0; }
331 info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info);
336 static int VechMatGetRank(
void *AA,
int *rank,
int n){
337 vechmat* A=(vechmat*)AA;
338 switch (A->factored){
349 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
354 static int VechMatGetEig(
void* AA,
int rank,
double *eigenvalue,
double vv[],
int n,
int indx[],
int *nind){
355 vechmat* A=(vechmat*)AA;
356 const double *val=A->val,tt=sqrt(0.5);
357 int info,i,j,k,ishift=A->ishift;
358 const int *ind=A->ind;
361 switch (A->factored){
363 memset(vv,0,n*
sizeof(
double));
364 getij(ind[rank]-ishift,&i,&j);
366 *eigenvalue=val[rank]*A->alpha;
371 memset(vv,0,n*
sizeof(
double));
373 getij(ind[k]-ishift,&i,&j);
376 vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
384 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
386 indx[0]=i; indx[1]=j;
388 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
390 indx[0]=i; indx[1]=j;
395 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
396 *eigenvalue=*eigenvalue*A->alpha;
399 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
405 static int VechMatView(
void* AA){
406 vechmat* A=(vechmat*)AA;
407 int info,i=0,j,k,rank=0,ishift=A->ishift,nnz=A->nnzeros;
408 const int *ind=A->ind;
409 const double *val=A->val;
410 for (k=0; k<nnz; k++){
411 getij(ind[k]-ishift,&i,&j);
412 printf(
"Row: %d, Column: %d, Value: %10.8f \n",i,j,A->alpha*val[k]);
415 info=VechMatGetRank(AA,&rank,A->n);DSDPCHKERR(info);
416 printf(
"Detected Rank: %d\n",rank);
423 static const char *datamatname=
"STANDARD VECH MATRIX";
427 if (sops==NULL)
return 0;
429 sops->matvecvec=VechMatVecVec;
430 sops->matdot=VechMatDot;
431 sops->matfnorm2=VechMatFNorm2;
432 sops->mataddrowmultiple=VechMatAddRowMultiple;
433 sops->mataddallmultiple=VechMatAddMultiple;
434 sops->matview=VechMatView;
435 sops->matdestroy=VechMatDestroy;
436 sops->matfactor2=VechMatFactor;
437 sops->matgetrank=VechMatGetRank;
438 sops->matgeteig=VechMatGetEig;
439 sops->matrownz=VechMatGetRowNnz;
440 sops->matnnz=VechMatCountNonzeros;
442 sops->matname=datamatname;
459 #define __FUNCT__ "DSDPGetVechMat" 461 int info,i,j,k,itmp,nn=n*(n+1)/2;
472 DSDPSETERR3(2,
"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
474 DSDPSETERR1(2,
"Illegal index value: %d. Must be >= 0\n",itmp);
477 for (k=0;k<nnz;++k) dtmp=val[k];
478 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
481 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
482 if (sops){*sops=&vechmatops;}
483 if (smat){*smat=(
void*)AA;}
484 DSDPFunctionReturn(0);
489 #define __FUNCT__ "VechMatComputeEigs" 490 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){
492 int i,j,k,nsub,neigs,info,*iptr,*perm,*invp;
493 long int *i2darray=(
long int*)DD;
494 int ownarray1=0,ownarray2=0,ownarray3=0;
495 int ishift=AA->ishift,nonzeros=AA->nnzeros;
496 const int *ind=AA->ind;
497 const double *val=AA->val;
498 double *dmatarray=ss1,*dworkarray=ss2,maxeig,eps=1.0e-12,eps2=1.0e-12;
500 iptr=iiptr; perm=iptr+n; invp=perm+n;
512 for (nsub=0,i=0; i<n; i++){
513 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
518 DSDPCALLOC2(&dmatarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
521 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
523 DSDPCALLOC2(&dworkarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
527 if (nsub*nsub*
sizeof(
long int)>nn0*
sizeof(
double)){
528 DSDPCALLOC2(&i2darray,
long int,(nsub*nsub),&info); DSDPCHKERR(info);
532 for (i=0,k=0; k<nonzeros; k++){
533 getij(ind[k]-ishift,&i,&j);
534 dmatarray[perm[i]*nsub+perm[j]] += val[k];
536 dmatarray[perm[j]*nsub+perm[i]] += val[k];
540 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
541 W,n,WORK,n1,iiptr+3*n,n2-3*n);
543 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
544 for (i=0,k=0; k<nonzeros; k++){
545 getij(ind[k]-ishift,&i,&j);
546 dmatarray[perm[i]*nsub+perm[j]] += val[k];
547 if (i!=j){dmatarray[perm[j]*nsub+perm[i]] += val[k];}
549 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
550 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
552 for (maxeig=0,i=0;i<nsub;i++){
553 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
555 memset((
void*)iptr,0,nsub*
sizeof(
int));
558 for (neigs=0,k=0; k<nsub; k++){
559 if (fabs(W[k]) > eps){
560 for (j=0;j<nsub;j++){
561 if (fabs(dmatarray[nsub*k+j]) >= eps2){iptr[neigs]++;
562 }
else { dmatarray[nsub*k+j]=0.0;}
572 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
573 DSDPLogInfo(0,49,
" Data Mat has %d eigenvalues: \n",neigs);
575 for (neigs=0,i=0; i<nsub; i++){
576 if (fabs(W[i]) > eps){
577 info=EigMatSetEig(AA->Eig,neigs,W[i],invp,dmatarray+nsub*i,nsub,n);DSDPCHKERR(info);
582 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
583 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
584 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 DSDPGetVechMat(int n, int ishift, double alpha, const int ind[], const double val[], int nnz, struct DSDPDataMat_Ops **sops, void **smat)
Given data in packed symmetric format, create a sparse matrix usuable by DSDP.