DSDP
vech.c
Go to the documentation of this file.
1 #include "dsdpdatamat_impl.h"
2 #include "dsdpsys.h"
7 typedef struct {
8  int neigs;
9  double *eigval;
10  double *an;
11  int *cols,*nnz;
12 } Eigen;
13 
19 typedef struct {
20  int nnzeros;
21  const int *ind;
22  const double *val;
23  int ishift;
24  double alpha;
25 
26  Eigen *Eig;
27  int factored;
28  int owndata;
29  int n;
30 } vechmat;
31 
32 #define GETI(a) (int)(sqrt(2*(a)+0.25)-0.5)
33 #define GETJ(b,c) ((b)-((c)*((c)+1))/2)
34 
35 static void getij(int k, int *i,int *j){
36  *i=GETI(k);
37  *j=GETJ(k,*i);
38  return;
39 }
40 /*
41 static void geti2(int k, int *i,int *j){
42  int r,c;
43  double rr=sqrt(2*k+0.25)-0.5;
44  r=(int)rr;
45  c=k-(r*(r+1))/2;
46  *i=r;*j=c;
47  return;
48 }
49 */
50 #undef __FUNCT__
51 #define __FUNCT__ "CreateVechMatWData"
52 static int CreateVechMatWdata(int n, int ishift, double alpha,const int *ind, const double *vals, int nnz, vechmat **A){
53  int info;
54  vechmat* V;
55  DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
56  V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
57  V->alpha=alpha;
58  V->owndata=0;
59  V->Eig=0;
60  *A=V;
61  return 0;
62 }
63 
64 
65 static int VechMatAddMultiple(void* AA, double scl, double r[], int nn, int n){
66  vechmat* A=(vechmat*)AA;
67  int k,nnz=A->nnzeros;
68  const int* ind=A->ind;
69  const double *val=A->val;
70  double *rr=r-A->ishift;
71  scl*=A->alpha;
72  for (k=0; k<nnz; ++k,++ind,++val){
73  *(rr+((*ind))) +=scl*(*val);
74  }
75  return 0;
76 }
77 
78 static int VechMatDot(void* AA, double x[], int nn, int n, double *v){
79  vechmat* A=(vechmat*)AA;
80  int k,nnz=A->nnzeros;
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)));
86  }
87  *v=2*vv*A->alpha;
88  return 0;
89 }
90 
91 static int EigMatVecVec(Eigen*, double[], int, double*);
92 static int VechMatGetRank(void*,int*,int);
93 
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;
98  double vv=0,dd;
99  const double *val=A->val,nnz=A->nnzeros;
100 
101  if (A->factored==3){
102  info=VechMatGetRank(AA,&rank,n);
103  if (nnz>3 && rank<nnz){
104  info=EigMatVecVec(A->Eig,x,n,&vv);
105  *v=vv*A->alpha;
106  return 0;
107  }
108  }
109 
110  for (k=0; k<nnz; ++k,++ind,++val){
111  kk=*ind-ishift;
112  i=GETI(kk);
113  j=GETJ(kk,i);
114  dd=x[i]*x[j]*(*val);
115  vv+=2*dd;
116  if (i==j){ vv-=dd; }
117  }
118  *v=vv*A->alpha;
119 
120  return 0;
121 }
122 
123 
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;
128  *nnzz=0;
129  for (k=0; k<nnz; ++k,++ind){
130  getij((*ind)-ishift,&i,&j);
131  if (i==trow){
132  nz[j]++;(*nnzz)++;
133  } else if (j==trow){
134  nz[i]++;(*nnzz)++;
135  }
136  }
137  return 0;
138 }
139 
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;
144  double fn2=0;
145  const double *val=A->val;
146  for (k=0; k<nnz; ++k,++ind){
147  getij((*ind)-ishift,&i,&j);
148  if (i==j){
149  fn2+= val[k]*val[k];
150  } else {
151  fn2+= 2*val[k]*val[k];
152  }
153  }
154  *fnorm2=fn2*A->alpha*A->alpha;
155  return 0;
156 }
157 
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;
163  scl*=A->alpha;
164  for (k=0; k<nnz; ++k,++ind){
165  getij((*ind)-ishift,&i,&j);
166  if (i==trow){
167  /* if (i==j){ r[j]+=scl*val[k];} */
168  r[j]+=scl*val[k];
169  } else if (j==trow){
170  r[i]+=scl*val[k];
171  }
172  }
173  return 0;
174 }
175 
176 static int VechMatCountNonzeros(void* AA, int*nnz, int n){
177  vechmat* A=(vechmat*)AA;
178  *nnz=A->nnzeros;
179  return 0;
180 }
181 
182 #undef __FUNCT__
183 #define __FUNCT__ "VechMatDestroy"
184 static int VechMatDestroy(void* AA){
185  vechmat* A=(vechmat*)AA;
186  int info;
187  if (A->owndata){
188  /* Never happens
189  if (A->ind){ DSDPFREE(&A->ind,&info);DSDPCHKERR(info);}
190  if (A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
191  */
192  return 1;
193  }
194  if (A->Eig){
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);
200  }
201  DSDPFREE(&A,&info);DSDPCHKERR(info);
202  return 0;
203 }
204 
205 
206 
207 #undef __FUNCT__
208 #define __FUNCT__ "DSDPCreateVechMatEigs"
209 static int CreateEigenLocker(Eigen **EE,int iptr[], int neigs, int n){
210  int i,k,info;
211  Eigen *E;
212 
213  for (k=0,i=0;i<neigs;i++) k+=iptr[i];
214  if (k>n*neigs/4){
215 
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);
219  E->neigs=neigs;
220  E->cols=0;
221  E->nnz=0;
222 
223  } else {
224 
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);
230  E->neigs=neigs;
231 
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];}
234  }
235  *EE=E;
236  return 0;
237 }
238 
239 
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;
242  double *an=A->an;
243 
244  A->eigval[row]=eigv;
245  if (cols){
246  k=0; if (row>0){ k=A->nnz[row-1];}
247  cols+=k; an+=k;
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++;
251  }
252  } else {
253  an+=n*row;
254  for (j=0; j<nsub; j++){
255  if (v[j]==0.0) continue;
256  an[idxn[j]]=v[j];
257  }
258  }
259  return 0;
260 }
261 
262 
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;
265  double* an=A->an;
266  *eigenvalue=A->eigval[row];
267  *nind=0;
268  if (cols){
269  memset((void*)eigenvector,0,n*sizeof(double));
270  if (row==0){ bb=0;} else {bb=A->nnz[row-1];} ee=A->nnz[row];
271  for (i=bb;i<ee;i++){
272  eigenvector[cols[i]]=an[i];
273  spind[i-bb]=cols[i]; (*nind)++;
274  }
275  } else {
276  memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
277  for (i=0;i<n;i++)spind[i]=i;
278  *nind=n;
279  }
280  return 0;
281 }
282 
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;
286 
287  if (cols){
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]];
292  }
293  ddd+=dd*dd*eigval[rank];
294  }
295  } else {
296  for (rank=0;rank<neigs;rank++){
297  for (dd=0,i=0;i<n;i++){
298  dd+=an[i]*v[i];
299  }
300  an+=n;
301  ddd+=dd*dd*eigval[rank];
302  }
303  }
304  *vv=ddd;
305  return 0;
306 }
307 
308 
309 static int VechMatComputeEigs(vechmat*,double[],int,double[],int,double[],int,int[],int,double[],int,double[],int);
310 
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;
314  int nn1=0,nn2=0;
315  double *ss1=0,*ss2=0;
316  const int *ind=A->ind;
317 
318  if (A->factored) return 0;
319  memset((void*)iptr,0,3*n*sizeof(int));
320  /* Find number of nonzeros in each row */
321  for (isdiag=1,k=0; k<nonzeros; k++){
322  getij(ind[k]-ishift,&i,&j);
323  iptr[i]++;
324  if (i!=j) {isdiag=0;iptr[j]++;}
325  }
326 
327  if (isdiag){ A->factored=1; return 0;}
328  /* Find most nonzeros per row */
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);
332  A->factored=3;
333  return 0;
334 }
335 
336 static int VechMatGetRank(void *AA,int *rank,int n){
337  vechmat* A=(vechmat*)AA;
338  switch (A->factored){
339  case 1:
340  *rank=A->nnzeros;
341  break;
342  case 2:
343  *rank=2*A->nnzeros;
344  break;
345  case 3:
346  *rank=A->Eig->neigs;
347  break;
348  default:
349  DSDPSETERR(1,"Vech Matrix not factored yet\n");
350  }
351  return 0;
352 }
353 
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;
359 
360  *nind=0;
361  switch (A->factored){
362  case 1:
363  memset(vv,0,n*sizeof(double));
364  getij(ind[rank]-ishift,&i,&j);
365  vv[i]=1.0;
366  *eigenvalue=val[rank]*A->alpha;
367  *nind=1;
368  indx[0]=i;
369  break;
370  case 2:
371  memset(vv,0,n*sizeof(double));
372  k=rank/2;
373  getij(ind[k]-ishift,&i,&j);
374  if (i==j){
375  if (k*2==rank){
376  vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
377  *nind=1;
378  indx[0]=i;
379  } else {
380  *eigenvalue=0;
381  }
382  } else {
383  if (k*2==rank){
384  vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
385  *nind=2;
386  indx[0]=i; indx[1]=j;
387  } else {
388  vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
389  *nind=2;
390  indx[0]=i; indx[1]=j;
391  }
392  }
393  break;
394  case 3:
395  info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
396  *eigenvalue=*eigenvalue*A->alpha;
397  break;
398  default:
399  DSDPSETERR(1,"Vech Matrix not factored yet\n");
400  }
401 
402  return 0;
403 }
404 
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]);
413  }
414  if (A->factored>0){
415  info=VechMatGetRank(AA,&rank,A->n);DSDPCHKERR(info);
416  printf("Detected Rank: %d\n",rank);
417  }
418  return 0;
419 }
420 
421 
422 static struct DSDPDataMat_Ops vechmatops;
423 static const char *datamatname="STANDARD VECH MATRIX";
424 
425 static int VechMatOpsInitialize(struct DSDPDataMat_Ops *sops){
426  int info;
427  if (sops==NULL) return 0;
428  info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
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;
441  sops->id=3;
442  sops->matname=datamatname;
443  return 0;
444 }
445 
458 #undef __FUNCT__
459 #define __FUNCT__ "DSDPGetVechMat"
460 int DSDPGetVechMat(int n,int ishift,double alpha, const int ind[], const double val[],int nnz, struct DSDPDataMat_Ops**sops, void**smat){
461  int info,i,j,k,itmp,nn=n*(n+1)/2;
462  double dtmp;
463  vechmat* AA;
464  DSDPFunctionBegin;
465  for (k=0;k<nnz;++k){
466  itmp=ind[k]-ishift;
467  if (itmp>=nn){
468  getij(itmp,&i,&j);
469  /*
470  DSDPSETERR(2,"Illegal index value: Element %d in array has row %d (>0) or column %d (>0) is greater than %d. \n",k+1,i+1,j+1,n);
471  */
472  DSDPSETERR3(2,"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
473  } else if (itmp<0){
474  DSDPSETERR1(2,"Illegal index value: %d. Must be >= 0\n",itmp);
475  }
476  }
477  for (k=0;k<nnz;++k) dtmp=val[k];
478  info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
479  AA->factored=0;
480  AA->Eig=0;
481  info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
482  if (sops){*sops=&vechmatops;}
483  if (smat){*smat=(void*)AA;}
484  DSDPFunctionReturn(0);
485 }
486 
487 
488 #undef __FUNCT__
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){
491 
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;
499 
500  iptr=iiptr; perm=iptr+n; invp=perm+n;
501  /* These operations were done before calling this routine * /
502  / * Integer arrays corresponding to rows with nonzeros and inverse map * /
503  memset((void*)iiptr,0,3*n*sizeof(int));
504 
505  / * Find number of nonzeros in each row * /
506  for (i=0,k=0; k<nonzeros; k++){
507  getij(ind[k],i,n,&i,&j);
508  iptr[i]++; iptr[j]++;
509  }
510  */
511  /* Number of rows with a nonzero. Order the rows with nonzeros. */
512  for (nsub=0,i=0; i<n; i++){
513  if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
514  }
515 
516  /* create a dense array in which to put numbers */
517  if (nsub*nsub>nn1){
518  DSDPCALLOC2(&dmatarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
519  ownarray1=1;
520  }
521  memset((void*)dmatarray,0,nsub*nsub*sizeof(double));
522  if (nsub*nsub>nn2){
523  DSDPCALLOC2(&dworkarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
524  ownarray2=1;
525  }
526 
527  if (nsub*nsub*sizeof(long int)>nn0*sizeof(double)){
528  DSDPCALLOC2(&i2darray,long int,(nsub*nsub),&info); DSDPCHKERR(info);
529  ownarray3=1;
530  }
531 
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];
535  if (i!=j){
536  dmatarray[perm[j]*nsub+perm[i]] += val[k];
537  }
538  }
539  /* Call LAPACK to compute the eigenvalues */
540  info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
541  W,n,WORK,n1,iiptr+3*n,n2-3*n);
542  if (info){
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];}
548  }
549  info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
550  W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
551  }
552  for (maxeig=0,i=0;i<nsub;i++){
553  if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
554  }
555  memset((void*)iptr,0,nsub*sizeof(int));
556  /* Compute sparsity pattern for eigenvalue and eigenvector structures */
557  /* Count the nonzero eigenvalues */
558  for (neigs=0,k=0; k<nsub; k++){
559  if (fabs(W[k]) /* /maxeig */ > 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;}
563  }
564  neigs++;
565  /*
566  } else if (fabs(W[k])>1.0e-100){
567  printf("SKIPPING EIGENVALUE: %4.4e, max is : %4.4e\n",W[k],maxeig);
568  */
569  }
570  }
571 
572  info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
573  DSDPLogInfo(0,49," Data Mat has %d eigenvalues: \n",neigs);
574  /* Copy into structure */
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);
578  neigs++;
579  }
580  }
581 
582  if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
583  if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
584  if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
585  return 0;
586 }
587 
Error handling, printing, and profiling.
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition: dsdpdatamat.c:47
Structure of function pointers that each SDP data matrix type (sparse, dense, constant, identity, ...) must implement.
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.
Definition: vech.c:460