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 
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,...
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.
Definition: vechu.c:454