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