DSDP
vechu.c
Go to the documentation of this file.
1#include "dsdpdatamat_impl.h"
2#include "dsdpsys.h"
7typedef struct {
8 int neigs;
9 double *eigval;
10 double *an;
11 int *cols,*nnz;
12} Eigen;
13
14typedef 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
30static 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"
38static 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
49static 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
62static 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
75static int EigMatVecVec(Eigen*, double[], int, double*);
76static int VechMatGetRank(void*,int*,int);
77
78static 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
108static 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
126static 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
146static 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
165static 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"
173static 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"
198static 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
229static 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
251static 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
271static 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
297static int VechMatComputeEigs(vechmat*,double[],int,double[],int,double[],int,int[],int,double[],int,double[],int);
298
299static 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
328static 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
346static 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
399static 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
416static struct DSDPDataMat_Ops vechmatops;
417static const char *datamatname="STANDARD VECH MATRIX";
418
419static 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"
454int 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"
484static 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
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,...
Error handling, printing, and profiling.
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