27 #define GETI(a,b) (int)((int)a/(int)b)
28 #define GETJ(a,b) (int)((int)a%(int)b)
30 static void getij(
int k,
int n,
int *i,
int *j){
37 #define __FUNCT__ "CreateVechMatWData"
38 static int CreateVechMatWdata(
int n,
int ishift,
double alpha,
const int *ind,
const double *vals,
int nnz, vechmat **A){
41 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
42 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
49 static int VechMatAddMultiple(
void* AA,
double scl,
double r[],
int nn,
int n){
50 vechmat* A=(vechmat*)AA;
52 const int *ind=A->ind,nnz=A->nnzeros;
53 const double *val=A->val;
54 double *rr=r-A->ishift;
56 for (k=0; k<nnz; ++k){
57 *(rr+(ind[k])) +=scl*(val[k]);
62 static int VechMatDot(
void* AA,
double x[],
int nn,
int n,
double *v){
63 vechmat* A=(vechmat*)AA;
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)));
75 static int EigMatVecVec(Eigen*,
double[],
int,
double*);
76 static int VechMatGetRank(
void*,
int*,
int);
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;
83 const double *val=A->val;
86 info=VechMatGetRank(AA,&rank,n);
87 if (nnz>3 && rank<nnz){
88 info=EigMatVecVec(A->Eig,x,n,&vv);
94 for (k=0; k<nnz; ++k,++ind,++val){
108 static int VechMatGetRowNnz(
void* AA,
int trow,
int nz[],
int *nnzz,
int nn){
109 vechmat* A=(vechmat*)AA;
111 const int *ind=A->ind, ishift=A->ishift,nnz=A->nnzeros;
113 for (k=0; k<nnz; ++k,ind){
126 static int VechMatFNorm2(
void* AA,
int n,
double *fnorm2){
127 vechmat* A=(vechmat*)AA;
129 const int *ind=A->ind,ishift=A->ishift,nnz=A->nnzeros;
131 const double *val=A->val;
132 for (k=0; k<nnz; ++k){
139 fn2+= 2*val[k]*val[k];
142 *fnorm2=fn2*A->alpha*A->alpha;
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;
152 for (k=0; k<nnz; ++k){
165 static int VechMatCountNonzeros(
void* AA,
int*nnz,
int n){
166 vechmat* A=(vechmat*)AA;
172 #define __FUNCT__ "VechMatDestroy"
173 static int VechMatDestroy(
void* AA){
174 vechmat* A=(vechmat*)AA;
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);
190 DSDPFREE(&A,&info);DSDPCHKERR(info);
197 #define __FUNCT__ "DSDPCreateVechMatEigs"
198 static int CreateEigenLocker(Eigen **EE,
int iptr[],
int neigs,
int n){
202 for (k=0,i=0;i<neigs;i++) k+=iptr[i];
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);
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);
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];}
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;
234 k=0;
if (row>0){ k=A->nnz[row-1];}
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++;
242 for (j=0; j<nsub; j++){
243 if (v[j]==0.0)
continue;
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;
254 *eigenvalue=A->eigval[row];
257 memset((
void*)eigenvector,0,n*
sizeof(
double));
258 if (row==0){ bb=0;}
else {bb=A->nnz[row-1];} ee=A->nnz[row];
260 eigenvector[cols[i]]=an[i];
261 spind[i-bb]=cols[i]; (*nind)++;
264 memcpy((
void*)eigenvector,(
void*)(an+n*row),n*
sizeof(
double));
265 for (i=0;i<n;i++)spind[i]=i;
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;
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]];
281 ddd+=dd*dd*eigval[rank];
284 for (rank=0;rank<neigs;rank++){
285 for (dd=0,i=0;i<n;i++){
289 ddd+=dd*dd*eigval[rank];
297 static int VechMatComputeEigs(vechmat*,
double[],
int,
double[],
int,
double[],
int,
int[],
int,
double[],
int,
double[],
int);
299 static int VechMatFactor(
void*AA,
double dmatp[],
int nn0,
double dwork[],
int n,
double ddwork[],
int n1,
int iptr[],
int n2){
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;
306 if (A->factored)
return 0;
308 memset((
void*)iptr,0,3*n*
sizeof(
int));
310 for (isdiag=1,k=0; k<nonzeros; k++){
315 if (i!=j) {isdiag=0;iptr[j]++;}
318 if (isdiag){ A->factored=1;
return 0;}
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; }
323 info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info);
328 static int VechMatGetRank(
void *AA,
int *rank,
int n){
329 vechmat* A=(vechmat*)AA;
330 switch (A->factored){
341 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
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);
350 const int *ind=A->ind,ishift=A->ishift;
353 switch (A->factored){
355 memset(vv,0,n*
sizeof(
double));
360 *eigenvalue=val[rank]*A->alpha;
365 memset(vv,0,n*
sizeof(
double));
367 getij(ind[k]-ishift,n,&i,&j);
370 vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
378 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
380 indx[0]=i; indx[1]=j;
382 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
384 indx[0]=i; indx[1]=j;
389 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
390 *eigenvalue=*eigenvalue*A->alpha;
393 DSDPSETERR(1,
"Vech Matrix not factored yet\n");
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]);
409 info=VechMatGetRank(AA,&rank,n);DSDPCHKERR(info);
410 printf(
"Detected Rank: %d\n",rank);
417 static const char *datamatname=
"STANDARD VECH MATRIX";
421 if (sops==NULL)
return 0;
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;
436 sops->matname=datamatname;
453 #define __FUNCT__ "DSDPGetVecUMat"
455 int info,i,j,k,itmp,nn=n*n;
466 DSDPSETERR3(2,
"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
468 DSDPSETERR1(2,
"Illegal index value: %d. Must be >= 0\n",itmp);
471 for (k=0;k<nnz;++k) dtmp=val[k];
472 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
475 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
476 if (sops){*sops=&vechmatops;}
477 if (smat){*smat=(
void*)AA;}
478 DSDPFunctionReturn(0);
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){
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;
494 iptr=iiptr; perm=iptr+n; invp=perm+n;
506 for (nsub=0,i=0; i<n; i++){
507 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
512 DSDPCALLOC2(&dmatarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
515 memset((
void*)dmatarray,0,nsub*nsub*
sizeof(
double));
517 DSDPCALLOC2(&dworkarray,
double,(nsub*nsub),&info); DSDPCHKERR(info);
521 if (nsub*nsub*
sizeof(
long int)>nn0*
sizeof(
double)){
522 DSDPCALLOC2(&i2darray,
long int,(nsub*nsub),&info); DSDPCHKERR(info);
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];
531 dmatarray[perm[j]*nsub+perm[i]] += val[k];
535 memset((
void*)W,0,n*
sizeof(
double));
537 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
538 W,nsub,WORK,n1,iiptr+3*n,n2-3*n);
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];
545 dmatarray[perm[j]*nsub+perm[i]] += val[k];
548 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
549 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
553 for (maxeig=0,i=0;i<nsub;i++){
554 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
556 memset((
void*)iptr,0,nsub*
sizeof(
int));
559 for (neigs=0,k=0; k<nsub; k++){
560 if (fabs(W[k]) > 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;}
573 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
574 DSDPLogInfo(0,49,
" Data Mat has %d eigenvectors.",neigs);
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);
583 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
584 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
585 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}