DSDP
vechu.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 
14 typedef struct {
15  int nnzeros;
16  const int *ind;
17  const double *val;
18  int ishift;
19  double alpha;
20 
21  Eigen *Eig;
22  int factored;
23  int owndata;
24  int n;
25 } vechmat;
26 
27 #define GETI(a,b) (int)((int)a/(int)b)
28 #define GETJ(a,b) (int)((int)a%(int)b)
29 
30 static void getij(int k, int n, int *i,int *j){
31  *i=GETI(k,n);
32  *j=GETJ(k,n);
33  return;
34 }
35 
36 #undef __FUNCT__
37 #define __FUNCT__ "CreateVechMatWData"
38 static int CreateVechMatWdata(int n, int ishift, double alpha, const int *ind, const double *vals, int nnz, vechmat **A){
39  int info;
40  vechmat* V;
41  DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
42  V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
43  V->alpha=alpha;
44  V->owndata=0;
45  *A=V;
46  return 0;
47 }
48 
49 static int VechMatAddMultiple(void* AA, double scl, double r[], int nn, int n){
50  vechmat* A=(vechmat*)AA;
51  int k;
52  const int *ind=A->ind,nnz=A->nnzeros;
53  const double *val=A->val;
54  double *rr=r-A->ishift;
55  scl*=A->alpha;
56  for (k=0; k<nnz; ++k){
57  *(rr+(ind[k])) +=scl*(val[k]);
58  }
59  return 0;
60 }
61 
62 static int VechMatDot(void* AA, double x[], int nn, int n, double *v){
63  vechmat* A=(vechmat*)AA;
64  int k,nnz=A->nnzeros;
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)));
70  }
71  *v=2*vv*A->alpha;
72  return 0;
73 }
74 
75 static int EigMatVecVec(Eigen*, double[], int, double*);
76 static int VechMatGetRank(void*,int*,int);
77 
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;
82  double vv=0,dd;
83  const double *val=A->val;
84 
85  if (A->factored==3){
86  info=VechMatGetRank(AA,&rank,n);
87  if (nnz>3 && rank<nnz){
88  info=EigMatVecVec(A->Eig,x,n,&vv);
89  *v=vv*A->alpha;
90  return 0;
91  }
92  }
93 
94  for (k=0; k<nnz; ++k,++ind,++val){
95  kk=*ind-ishift;
96  i=GETI(kk,n);
97  j=GETJ(kk,n);
98  dd=x[i]*x[j]*(*val);
99  vv+=2*dd;
100  if (i==j){ vv-=dd; }
101  }
102  *v=vv*A->alpha;
103 
104  return 0;
105 }
106 
107 
108 static int VechMatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int nn){
109  vechmat* A=(vechmat*)AA;
110  int i=0,j,k,t;
111  const int *ind=A->ind, ishift=A->ishift,nnz=A->nnzeros;
112  *nnzz=0;
113  for (k=0; k<nnz; ++k,ind){
114  t=ind[k]-ishift;
115  i=GETI(t,nn);
116  j=GETJ(t,nn);
117  if (i==trow){
118  nz[j]++;(*nnzz)++;
119  } else if (j==trow){
120  nz[i]++;(*nnzz)++;
121  }
122  }
123  return 0;
124 }
125 
126 static int VechMatFNorm2(void* AA, int n, double *fnorm2){
127  vechmat* A=(vechmat*)AA;
128  int i=0,j,k,t;
129  const int *ind=A->ind,ishift=A->ishift,nnz=A->nnzeros;
130  double fn2=0;
131  const double *val=A->val;
132  for (k=0; k<nnz; ++k){
133  t=ind[k]-ishift;
134  i=GETI(t,n);
135  j=GETJ(t,n);
136  if (i==j){
137  fn2+= val[k]*val[k];
138  } else {
139  fn2+= 2*val[k]*val[k];
140  }
141  }
142  *fnorm2=fn2*A->alpha*A->alpha;
143  return 0;
144 }
145 
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;
151  scl*=A->alpha;
152  for (k=0; k<nnz; ++k){
153  t=ind[k]-ishift;
154  i=GETI(t,m);
155  j=GETJ(t,m);
156  if (i==trow){
157  r[j]+=scl*val[k];
158  } else if (j==trow){
159  r[i]+=scl*val[k];
160  }
161  }
162  return 0;
163 }
164 
165 static int VechMatCountNonzeros(void* AA, int*nnz, int n){
166  vechmat* A=(vechmat*)AA;
167  *nnz=A->nnzeros;
168  return 0;
169 }
170 
171 #undef __FUNCT__
172 #define __FUNCT__ "VechMatDestroy"
173 static int VechMatDestroy(void* AA){
174  vechmat* A=(vechmat*)AA;
175  int info;
176  if (A->owndata){
177  /*
178  if (A->ind){ DSDPFREE(&A->ind,&info);DSDPCHKERR(info);}
179  if (A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
180  */
181  return 1;
182  }
183  if (A->Eig){
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);
189  }
190  DSDPFREE(&A,&info);DSDPCHKERR(info);
191  return 0;
192 }
193 
194 
195 
196 #undef __FUNCT__
197 #define __FUNCT__ "DSDPCreateVechMatEigs"
198 static int CreateEigenLocker(Eigen **EE,int iptr[], int neigs, int n){
199  int i,k,info;
200  Eigen *E;
201 
202  for (k=0,i=0;i<neigs;i++) k+=iptr[i];
203  if (k>n*neigs/4){
204 
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);
208  E->neigs=neigs;
209  E->cols=0;
210  E->nnz=0;
211 
212  } else {
213 
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);
219  E->neigs=neigs;
220 
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];}
223  }
224  *EE=E;
225  return 0;
226 }
227 
228 
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;
231  double *an=A->an;
232  A->eigval[row]=eigv;
233  if (cols){
234  k=0; if (row>0){ k=A->nnz[row-1];}
235  cols+=k; an+=k;
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++;
239  }
240  } else {
241  an+=n*row;
242  for (j=0; j<nsub; j++){
243  if (v[j]==0.0) continue;
244  an[idxn[j]]=v[j];
245  }
246  }
247  return 0;
248 }
249 
250 
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;
253  double* an=A->an;
254  *eigenvalue=A->eigval[row];
255  *nind=0;
256  if (cols){
257  memset((void*)eigenvector,0,n*sizeof(double));
258  if (row==0){ bb=0;} else {bb=A->nnz[row-1];} ee=A->nnz[row];
259  for (i=bb;i<ee;i++){
260  eigenvector[cols[i]]=an[i];
261  spind[i-bb]=cols[i]; (*nind)++;
262  }
263  } else {
264  memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
265  for (i=0;i<n;i++)spind[i]=i;
266  *nind=n;
267  }
268  return 0;
269 }
270 
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;
274 
275  if (cols){
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]];
280  }
281  ddd+=dd*dd*eigval[rank];
282  }
283  } else {
284  for (rank=0;rank<neigs;rank++){
285  for (dd=0,i=0;i<n;i++){
286  dd+=an[i]*v[i];
287  }
288  an+=n;
289  ddd+=dd*dd*eigval[rank];
290  }
291  }
292  *vv=ddd;
293  return 0;
294 }
295 
296 
297 static int VechMatComputeEigs(vechmat*,double[],int,double[],int,double[],int,int[],int,double[],int,double[],int);
298 
299 static int VechMatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){
300 
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;
305  int nn1=0,nn2=0;
306  if (A->factored) return 0;
307 
308  memset((void*)iptr,0,3*n*sizeof(int));
309  /* Find number of nonzeros in each row */
310  for (isdiag=1,k=0; k<nonzeros; k++){
311  t=ind[k]-ishift;
312  i=GETI(t,n);
313  j=GETJ(t,n);
314  iptr[i]++;
315  if (i!=j) {isdiag=0;iptr[j]++;}
316  }
317 
318  if (isdiag){ A->factored=1; return 0;}
319  /* Find most nonzeros per row */
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; }
322 
323  info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info);
324  A->factored=3;
325  return 0;
326 }
327 
328 static int VechMatGetRank(void *AA,int *rank,int n){
329  vechmat* A=(vechmat*)AA;
330  switch (A->factored){
331  case 1:
332  *rank=A->nnzeros;
333  break;
334  case 2:
335  *rank=2*A->nnzeros;
336  break;
337  case 3:
338  *rank=A->Eig->neigs;
339  break;
340  default:
341  DSDPSETERR(1,"Vech Matrix not factored yet\n");
342  }
343  return 0;
344 }
345 
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);
349  int info,i,j,k,t;
350  const int *ind=A->ind,ishift=A->ishift;
351 
352  *nind=0;
353  switch (A->factored){
354  case 1:
355  memset(vv,0,n*sizeof(double));
356  t=ind[rank]-ishift;
357  i=GETI(t,n);
358  j=GETJ(t,n);
359  vv[i]=1.0;
360  *eigenvalue=val[rank]*A->alpha;
361  *nind=1;
362  indx[0]=i;
363  break;
364  case 2:
365  memset(vv,0,n*sizeof(double));
366  k=rank/2;
367  getij(ind[k]-ishift,n,&i,&j);
368  if (i==j){
369  if (k*2==rank){
370  vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
371  *nind=1;
372  indx[0]=i;
373  } else {
374  *eigenvalue=0;
375  }
376  } else {
377  if (k*2==rank){
378  vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
379  *nind=2;
380  indx[0]=i; indx[1]=j;
381  } else {
382  vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
383  *nind=2;
384  indx[0]=i; indx[1]=j;
385  }
386  }
387  break;
388  case 3:
389  info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
390  *eigenvalue=*eigenvalue*A->alpha;
391  break;
392  default:
393  DSDPSETERR(1,"Vech Matrix not factored yet\n");
394  }
395 
396  return 0;
397 }
398 
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]);
407  }
408  if (A->factored>0){
409  info=VechMatGetRank(AA,&rank,n);DSDPCHKERR(info);
410  printf("Detected Rank: %d\n",rank);
411  }
412  return 0;
413 }
414 
415 
416 static struct DSDPDataMat_Ops vechmatops;
417 static const char *datamatname="STANDARD VECH MATRIX";
418 
419 static int VechMatOpsInitialize(struct DSDPDataMat_Ops *sops){
420  int info;
421  if (sops==NULL) return 0;
422  info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
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;
435  sops->id=3;
436  sops->matname=datamatname;
437  return 0;
438 }
439 
452 #undef __FUNCT__
453 #define __FUNCT__ "DSDPGetVecUMat"
454 int DSDPGetVecUMat(int n,int ishift,double alpha,const int ind[], const double val[],int nnz, struct DSDPDataMat_Ops**sops, void**smat){
455  int info,i,j,k,itmp,nn=n*n;
456  double dtmp;
457  vechmat* AA;
458  DSDPFunctionBegin;
459  for (k=0;k<nnz;++k){
460  itmp=ind[k]-ishift;
461  if (itmp>=nn){
462  getij(itmp,n,&i,&j);
463  /*
464  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);
465  */
466  DSDPSETERR3(2,"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
467  } else if (itmp<0){
468  DSDPSETERR1(2,"Illegal index value: %d. Must be >= 0\n",itmp);
469  }
470  }
471  for (k=0;k<nnz;++k) dtmp=val[k];
472  info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
473  AA->factored=0;
474  AA->Eig=0;
475  info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
476  if (sops){*sops=&vechmatops;}
477  if (smat){*smat=(void*)AA;}
478  DSDPFunctionReturn(0);
479 }
480 
481 
482 #undef __FUNCT__
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){
485 
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;
493 
494  iptr=iiptr; perm=iptr+n; invp=perm+n;
495  /* These operations were done before calling this routine * /
496  / * Integer arrays corresponding to rows with nonzeros and inverse map * /
497  memset((void*)iiptr,0,3*n*sizeof(int));
498 
499  / * Find number of nonzeros in each row * /
500  for (i=0,k=0; k<nonzeros; k++){
501  getij(ind[k],i,n,&i,&j);
502  iptr[i]++; iptr[j]++;
503  }
504  */
505  /* Number of rows with a nonzero. Order the rows with nonzeros. */
506  for (nsub=0,i=0; i<n; i++){
507  if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
508  }
509 
510  /* create a dense array in which to put numbers */
511  if (nsub*nsub>nn1){
512  DSDPCALLOC2(&dmatarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
513  ownarray1=1;
514  }
515  memset((void*)dmatarray,0,nsub*nsub*sizeof(double));
516  if (nsub*nsub>nn2){
517  DSDPCALLOC2(&dworkarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
518  ownarray2=1;
519  }
520 
521  if (nsub*nsub*sizeof(long int)>nn0*sizeof(double)){
522  DSDPCALLOC2(&i2darray,long int,(nsub*nsub),&info); DSDPCHKERR(info);
523  ownarray3=1;
524  }
525 
526 
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];
530  if (i!=j){
531  dmatarray[perm[j]*nsub+perm[i]] += val[k];
532  }
533  }
534  /* Call LAPACK to compute the eigenvalues */
535  memset((void*)W,0,n*sizeof(double));
536 
537  info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
538  W,nsub,WORK,n1,iiptr+3*n,n2-3*n);
539  if (info){
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];
544  if (i!=j){
545  dmatarray[perm[j]*nsub+perm[i]] += val[k];
546  }
547  }
548  info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
549  W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
550  }
551  /* dsyev_("V","L",&N,dmatarray,&LDA,W,WORK,&LWORK,&INFO); */
552 
553  for (maxeig=0,i=0;i<nsub;i++){
554  if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
555  }
556  memset((void*)iptr,0,nsub*sizeof(int));
557  /* Compute sparsity pattern for eigenvalue and eigenvector structures */
558  /* Count the nonzero eigenvalues */
559  for (neigs=0,k=0; k<nsub; k++){
560  if (fabs(W[k]) /* /maxeig */ > 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;}
564  }
565  neigs++;
566  /*
567  } else if (fabs(W[k])>1.0e-100){
568  printf("SKIPPING EIGENVALUE: %4.4e, max is : %4.4e\n",W[k],maxeig);
569  */
570  }
571  }
572 
573  info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
574  DSDPLogInfo(0,49," Data Mat has %d eigenvectors.",neigs);
575  /* Copy into structure */
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);
579  neigs++;
580  }
581  }
582 
583  if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
584  if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
585  if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
586  return 0;
587 }
588 
dsdpdatamat_impl.h
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
DSDPGetVecUMat
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.
Definition: vechu.c:454
DSDPDataMatOpsInitialize
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition: dsdpdatamat.c:47
dsdpsys.h
Error handling, printing, and profiling.
DSDPDataMat_Ops
Table of function pointers that operate on the data matrix.
Definition: dsdpdatamat_impl.h:15