32 #define GETI(a) (int)(sqrt(2*(a)+0.25)-0.5)
33 #define GETJ(b,c) ((b)-((c)*((c)+1))/2)
35 static void getij(
int k,
int *i,
int *j){
51 #define __FUNCT__ "CreateVechMatWData"
52 static int CreateVechMatWdata(
int n,
int ishift,
double alpha,
const int *ind,
const double *vals,
int nnz, vechmat **A){
55 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
56 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
65 static int VechMatAddMultiple(
void* AA,
double scl,
double r[],
int nn,
int n){
66 vechmat* A=(vechmat*)AA;
68 const int* ind=A->ind;
69 const double *val=A->val;
70 double *rr=r-A->ishift;
72 for (k=0; k<nnz; ++k,++ind,++val){
73 *(rr+((*ind))) +=scl*(*val);
78 static int VechMatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
79 vechmat* A=(vechmat*)AA;
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)));
91 static int EigMatVecVec(Eigen*,
double[],
int,
double*);
92 static int VechMatGetRank(
void*,
int*,
int);
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;
99 const double *val=A->val,nnz=A->nnzeros;
102 info=VechMatGetRank(AA,&rank,n);
103 if (nnz>3 && rank<nnz){
104 info=EigMatVecVec(A->Eig,x,n,&vv);
110 for (k=0; k<nnz; ++k,++ind,++val){
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;
129 for (k=0; k<nnz; ++k,++ind){
130 getij((*ind)-ishift,&i,&j);
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;
145 const double *val=A->val;
146 for (k=0; k<nnz; ++k,++ind){
147 getij((*ind)-ishift,&i,&j);
151 fn2+= 2*val[k]*val[k];
154 *fnorm2=fn2*A->alpha*A->alpha;
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;
164 for (k=0; k<nnz; ++k,++ind){
165 getij((*ind)-ishift,&i,&j);
176 static int VechMatCountNonzeros(
void* AA,
int*nnz,
int n){
177 vechmat* A=(vechmat*)AA;
183 #define __FUNCT__ "VechMatDestroy"
184 static int VechMatDestroy(
void* AA){
185 vechmat* A=(vechmat*)AA;
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);
201 DSDPFREE(&A,&info);DSDPCHKERR(info);
208 #define __FUNCT__ "DSDPCreateVechMatEigs"
209 static int CreateEigenLocker(Eigen **EE,
int iptr[],
int neigs,
int n){
213 for (k=0,i=0;i<neigs;i++) k+=iptr[i];
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);
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);
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];}
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;
246 k=0;
if (row>0){ k=A->nnz[row-1];}
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++;
254 for (j=0; j<nsub; j++){
255 if (v[j]==0.0)
continue;
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;
266 *eigenvalue=A->eigval[row];
269 memset((
void*)eigenvector,0,n*
sizeof(
double));
270 if (row==0){ bb=0;}
else {bb=A->nnz[row-1];} ee=A->nnz[row];
272 eigenvector[cols[i]]=an[i];
273 spind[i-bb]=cols[i]; (*nind)++;
276 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
277 for (i=0;i<n;i++)spind[i]=i;
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;
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]];
293 ddd+=dd*dd*eigval[rank];
296 for (rank=0;rank<neigs;rank++){
297 for (dd=0,i=0;i<n;i++){
301 ddd+=dd*dd*eigval[rank];
309 static int VechMatComputeEigs(vechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int,
double[],
int,
double[],
int);
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;
315 double *ss1=0,*ss2=0;
316 const int *ind=A->ind;
318 if (A->factored)
return 0;
319 memset((
void*)iptr,0,3*n*
sizeof(
int));
321 for (isdiag=1,k=0; k<nonzeros; k++){
322 getij(ind[k]-ishift,&i,&j);
324 if (i!=j) {isdiag=0;iptr[j]++;}
327 if (isdiag){ A->factored=1;
return 0;}
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);
336 static int VechMatGetRank(
void *AA,
int *rank,
int n){
337 vechmat* A=(vechmat*)AA;
338 switch (A->factored){
349 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
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;
361 switch (A->factored){
363 memset(vv,0,n*
sizeof(
double));
364 getij(ind[rank]-ishift,&i,&j);
366 *eigenvalue=val[rank]*A->alpha;
371 memset(vv,0,n*
sizeof(
double));
373 getij(ind[k]-ishift,&i,&j);
376 vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
384 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
386 indx[0]=i; indx[1]=j;
388 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
390 indx[0]=i; indx[1]=j;
395 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
396 *eigenvalue=*eigenvalue*A->alpha;
399 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
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]);
415 info=VechMatGetRank(AA,&rank,A->n);DSDPCHKERR(info);
416 printf(
"Detected Rank: %d\n",rank);
423 static const char *datamatname=
"STANDARD VECH MATRIX";
427 if (sops==NULL)
return 0;
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;
442 sops->matname=datamatname;
459 #define __FUNCT__ "DSDPGetVechMat"
461 int info,i,j,k,itmp,nn=n*(n+1)/2;
472 DSDPSETERR3(2,
"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
474 DSDPSETERR1(2,
"Illegal index value: %d. Must be >= 0\n",itmp);
477 for (k=0;k<nnz;++k) dtmp=val[k];
478 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
481 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
482 if (sops){*sops=&vechmatops;}
483 if (smat){*smat=(
void*)AA;}
484 DSDPFunctionReturn(0);
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){
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;
500 iptr=iiptr; perm=iptr+n; invp=perm+n;
512 for (nsub=0,i=0; i<n; i++){
513 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
518 DSDPCALLOC2(&dmatarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
521 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
523 DSDPCALLOC2(&dworkarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
527 if (nsub*nsub*
sizeof(
long int)>nn0*
sizeof(
double)){
528 DSDPCALLOC2(&i2darray,
long int,(nsub*nsub),&info); DSDPCHKERR(info);
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];
536 dmatarray[perm[j]*nsub+perm[i]] += val[k];
540 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
541 W,n,WORK,n1,iiptr+3*n,n2-3*n);
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];}
549 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
550 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
552 for (maxeig=0,i=0;i<nsub;i++){
553 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
555 memset((
void*)iptr,0,nsub*
sizeof(
int));
558 for (neigs=0,k=0; k<nsub; k++){
559 if (fabs(W[k]) > 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;}
572 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
573 DSDPLogInfo(0,49,
" Data Mat has %d eigenvalues: \n",neigs);
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);
582 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
583 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
584 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}