DSDP
|
00001 #include "dsdpdatamat_impl.h" 00002 #include "dsdpsys.h" 00007 typedef struct { 00008 int neigs; 00009 double *eigval; 00010 double *an; 00011 int *cols,*nnz; 00012 } Eigen; 00013 00019 typedef struct { 00020 int nnzeros; 00021 const int *ind; 00022 const double *val; 00023 int ishift; 00024 double alpha; 00025 00026 Eigen *Eig; 00027 int factored; 00028 int owndata; 00029 int n; 00030 } vechmat; 00031 00032 #define GETI(a) (int)(sqrt(2*(a)+0.25)-0.5) 00033 #define GETJ(b,c) ((b)-((c)*((c)+1))/2) 00034 00035 static void getij(int k, int *i,int *j){ 00036 *i=GETI(k); 00037 *j=GETJ(k,*i); 00038 return; 00039 } 00040 /* 00041 static void geti2(int k, int *i,int *j){ 00042 int r,c; 00043 double rr=sqrt(2*k+0.25)-0.5; 00044 r=(int)rr; 00045 c=k-(r*(r+1))/2; 00046 *i=r;*j=c; 00047 return; 00048 } 00049 */ 00050 #undef __FUNCT__ 00051 #define __FUNCT__ "CreateVechMatWData" 00052 static int CreateVechMatWdata(int n, int ishift, double alpha,const int *ind, const double *vals, int nnz, vechmat **A){ 00053 int info; 00054 vechmat* V; 00055 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info); 00056 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz; 00057 V->alpha=alpha; 00058 V->owndata=0; 00059 V->Eig=0; 00060 *A=V; 00061 return 0; 00062 } 00063 00064 00065 static int VechMatAddMultiple(void* AA, double scl, double r[], int nn, int n){ 00066 vechmat* A=(vechmat*)AA; 00067 int k,nnz=A->nnzeros; 00068 const int* ind=A->ind; 00069 const double *val=A->val; 00070 double *rr=r-A->ishift; 00071 scl*=A->alpha; 00072 for (k=0; k<nnz; ++k,++ind,++val){ 00073 *(rr+((*ind))) +=scl*(*val); 00074 } 00075 return 0; 00076 } 00077 00078 static int VechMatDot(void* AA, double x[], int nn, int n, double *v){ 00079 vechmat* A=(vechmat*)AA; 00080 int k,nnz=A->nnzeros; 00081 const int *ind=A->ind; 00082 double vv=0, *xx=x-A->ishift; 00083 const double *val=A->val; 00084 for (k=0;k<nnz;++k,++ind,++val){ 00085 vv+=(*val)*(*(xx+(*ind))); 00086 } 00087 *v=2*vv*A->alpha; 00088 return 0; 00089 } 00090 00091 static int EigMatVecVec(Eigen*, double[], int, double*); 00092 static int VechMatGetRank(void*,int*,int); 00093 00094 static int VechMatVecVec(void* AA, double x[], int n, double *v){ 00095 vechmat* A=(vechmat*)AA; 00096 int info,rank=n,i=0,j,k,kk; 00097 const int *ind=A->ind,ishift=A->ishift; 00098 double vv=0,dd; 00099 const double *val=A->val,nnz=A->nnzeros; 00100 00101 if (A->factored==3){ 00102 info=VechMatGetRank(AA,&rank,n); 00103 if (nnz>3 && rank<nnz){ 00104 info=EigMatVecVec(A->Eig,x,n,&vv); 00105 *v=vv*A->alpha; 00106 return 0; 00107 } 00108 } 00109 00110 for (k=0; k<nnz; ++k,++ind,++val){ 00111 kk=*ind-ishift; 00112 i=GETI(kk); 00113 j=GETJ(kk,i); 00114 dd=x[i]*x[j]*(*val); 00115 vv+=2*dd; 00116 if (i==j){ vv-=dd; } 00117 } 00118 *v=vv*A->alpha; 00119 00120 return 0; 00121 } 00122 00123 00124 static int VechMatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int nn){ 00125 vechmat* A=(vechmat*)AA; 00126 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros; 00127 const int *ind=A->ind; 00128 *nnzz=0; 00129 for (k=0; k<nnz; ++k,++ind){ 00130 getij((*ind)-ishift,&i,&j); 00131 if (i==trow){ 00132 nz[j]++;(*nnzz)++; 00133 } else if (j==trow){ 00134 nz[i]++;(*nnzz)++; 00135 } 00136 } 00137 return 0; 00138 } 00139 00140 static int VechMatFNorm2(void* AA, int n, double *fnorm2){ 00141 vechmat* A=(vechmat*)AA; 00142 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros; 00143 const int *ind=A->ind; 00144 double fn2=0; 00145 const double *val=A->val; 00146 for (k=0; k<nnz; ++k,++ind){ 00147 getij((*ind)-ishift,&i,&j); 00148 if (i==j){ 00149 fn2+= val[k]*val[k]; 00150 } else { 00151 fn2+= 2*val[k]*val[k]; 00152 } 00153 } 00154 *fnorm2=fn2*A->alpha*A->alpha; 00155 return 0; 00156 } 00157 00158 static int VechMatAddRowMultiple(void* AA, int trow, double scl, double r[], int m){ 00159 vechmat* A=(vechmat*)AA; 00160 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros; 00161 const int *ind=A->ind; 00162 const double *val=A->val; 00163 scl*=A->alpha; 00164 for (k=0; k<nnz; ++k,++ind){ 00165 getij((*ind)-ishift,&i,&j); 00166 if (i==trow){ 00167 /* if (i==j){ r[j]+=scl*val[k];} */ 00168 r[j]+=scl*val[k]; 00169 } else if (j==trow){ 00170 r[i]+=scl*val[k]; 00171 } 00172 } 00173 return 0; 00174 } 00175 00176 static int VechMatCountNonzeros(void* AA, int*nnz, int n){ 00177 vechmat* A=(vechmat*)AA; 00178 *nnz=A->nnzeros; 00179 return 0; 00180 } 00181 00182 #undef __FUNCT__ 00183 #define __FUNCT__ "VechMatDestroy" 00184 static int VechMatDestroy(void* AA){ 00185 vechmat* A=(vechmat*)AA; 00186 int info; 00187 if (A->owndata){ 00188 /* Never happens 00189 if (A->ind){ DSDPFREE(&A->ind,&info);DSDPCHKERR(info);} 00190 if (A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);} 00191 */ 00192 return 1; 00193 } 00194 if (A->Eig){ 00195 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info); 00196 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info); 00197 if (A->Eig->cols){DSDPFREE(&A->Eig->cols,&info);DSDPCHKERR(info);} 00198 if (A->Eig->nnz){DSDPFREE(&A->Eig->nnz,&info);DSDPCHKERR(info);} 00199 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info); 00200 } 00201 DSDPFREE(&A,&info);DSDPCHKERR(info); 00202 return 0; 00203 } 00204 00205 00206 00207 #undef __FUNCT__ 00208 #define __FUNCT__ "DSDPCreateVechMatEigs" 00209 static int CreateEigenLocker(Eigen **EE,int iptr[], int neigs, int n){ 00210 int i,k,info; 00211 Eigen *E; 00212 00213 for (k=0,i=0;i<neigs;i++) k+=iptr[i]; 00214 if (k>n*neigs/4){ 00215 00216 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info); 00217 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info); 00218 DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info); 00219 E->neigs=neigs; 00220 E->cols=0; 00221 E->nnz=0; 00222 00223 } else { 00224 00225 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info); 00226 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info); 00227 DSDPCALLOC2(&E->nnz,int,neigs,&info);DSDPCHKERR(info); 00228 DSDPCALLOC2(&E->an,double,k,&info);DSDPCHKERR(info); 00229 DSDPCALLOC2(&E->cols,int,k,&info);DSDPCHKERR(info); 00230 E->neigs=neigs; 00231 00232 if (neigs>0) E->nnz[0]=iptr[0]; 00233 for (i=1;i<neigs;i++){E->nnz[i]=E->nnz[i-1]+iptr[i];} 00234 } 00235 *EE=E; 00236 return 0; 00237 } 00238 00239 00240 static int EigMatSetEig(Eigen* A,int row, double eigv, int idxn[], double v[], int nsub,int n){ 00241 int j,k,*cols=A->cols; 00242 double *an=A->an; 00243 00244 A->eigval[row]=eigv; 00245 if (cols){ 00246 k=0; if (row>0){ k=A->nnz[row-1];} 00247 cols+=k; an+=k; 00248 for (k=0,j=0; j<nsub; j++){ 00249 if (v[j]==0.0) continue; 00250 cols[k]=idxn[j]; an[k]=v[j]; k++; 00251 } 00252 } else { 00253 an+=n*row; 00254 for (j=0; j<nsub; j++){ 00255 if (v[j]==0.0) continue; 00256 an[idxn[j]]=v[j]; 00257 } 00258 } 00259 return 0; 00260 } 00261 00262 00263 static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n, int spind[], int *nind){ 00264 int i,*cols=A->cols,bb,ee; 00265 double* an=A->an; 00266 *eigenvalue=A->eigval[row]; 00267 *nind=0; 00268 if (cols){ 00269 memset((void*)eigenvector,0,n*sizeof(double)); 00270 if (row==0){ bb=0;} else {bb=A->nnz[row-1];} ee=A->nnz[row]; 00271 for (i=bb;i<ee;i++){ 00272 eigenvector[cols[i]]=an[i]; 00273 spind[i-bb]=cols[i]; (*nind)++; 00274 } 00275 } else { 00276 memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double)); 00277 for (i=0;i<n;i++)spind[i]=i; 00278 *nind=n; 00279 } 00280 return 0; 00281 } 00282 00283 static int EigMatVecVec(Eigen* A, double v[], int n, double *vv){ 00284 int i,rank,*cols=A->cols,neigs=A->neigs,*nnz=A->nnz,bb,ee; 00285 double* an=A->an,*eigval=A->eigval,dd,ddd=0; 00286 00287 if (cols){ 00288 for (rank=0;rank<neigs;rank++){ 00289 if (rank==0){ bb=0;} else {bb=nnz[rank-1];} ee=nnz[rank]; 00290 for (dd=0,i=bb;i<ee;i++){ 00291 dd+=an[i]*v[cols[i]]; 00292 } 00293 ddd+=dd*dd*eigval[rank]; 00294 } 00295 } else { 00296 for (rank=0;rank<neigs;rank++){ 00297 for (dd=0,i=0;i<n;i++){ 00298 dd+=an[i]*v[i]; 00299 } 00300 an+=n; 00301 ddd+=dd*dd*eigval[rank]; 00302 } 00303 } 00304 *vv=ddd; 00305 return 0; 00306 } 00307 00308 00309 static int VechMatComputeEigs(vechmat*,double[],int,double[],int,double[],int,int[],int,double[],int,double[],int); 00310 00311 static int VechMatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){ 00312 vechmat* A=(vechmat*)AA; 00313 int i,j,k,ishift=A->ishift,isdiag,nonzeros=A->nnzeros,info; 00314 int nn1=0,nn2=0; 00315 double *ss1=0,*ss2=0; 00316 const int *ind=A->ind; 00317 00318 if (A->factored) return 0; 00319 memset((void*)iptr,0,3*n*sizeof(int)); 00320 /* Find number of nonzeros in each row */ 00321 for (isdiag=1,k=0; k<nonzeros; k++){ 00322 getij(ind[k]-ishift,&i,&j); 00323 iptr[i]++; 00324 if (i!=j) {isdiag=0;iptr[j]++;} 00325 } 00326 00327 if (isdiag){ A->factored=1; return 0;} 00328 /* Find most nonzeros per row */ 00329 for (j=0,i=0; i<n; i++){ if (iptr[i]>j) j=iptr[i]; } 00330 if (j<2){ A->factored=2; return 0; } 00331 info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info); 00332 A->factored=3; 00333 return 0; 00334 } 00335 00336 static int VechMatGetRank(void *AA,int *rank,int n){ 00337 vechmat* A=(vechmat*)AA; 00338 switch (A->factored){ 00339 case 1: 00340 *rank=A->nnzeros; 00341 break; 00342 case 2: 00343 *rank=2*A->nnzeros; 00344 break; 00345 case 3: 00346 *rank=A->Eig->neigs; 00347 break; 00348 default: 00349 DSDPSETERR(1,"Vech Matrix not factored yet\n"); 00350 } 00351 return 0; 00352 } 00353 00354 static int VechMatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indx[], int *nind){ 00355 vechmat* A=(vechmat*)AA; 00356 const double *val=A->val,tt=sqrt(0.5); 00357 int info,i,j,k,ishift=A->ishift; 00358 const int *ind=A->ind; 00359 00360 *nind=0; 00361 switch (A->factored){ 00362 case 1: 00363 memset(vv,0,n*sizeof(double)); 00364 getij(ind[rank]-ishift,&i,&j); 00365 vv[i]=1.0; 00366 *eigenvalue=val[rank]*A->alpha; 00367 *nind=1; 00368 indx[0]=i; 00369 break; 00370 case 2: 00371 memset(vv,0,n*sizeof(double)); 00372 k=rank/2; 00373 getij(ind[k]-ishift,&i,&j); 00374 if (i==j){ 00375 if (k*2==rank){ 00376 vv[i]=1.0; *eigenvalue=val[k]*A->alpha; 00377 *nind=1; 00378 indx[0]=i; 00379 } else { 00380 *eigenvalue=0; 00381 } 00382 } else { 00383 if (k*2==rank){ 00384 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha; 00385 *nind=2; 00386 indx[0]=i; indx[1]=j; 00387 } else { 00388 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha; 00389 *nind=2; 00390 indx[0]=i; indx[1]=j; 00391 } 00392 } 00393 break; 00394 case 3: 00395 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info); 00396 *eigenvalue=*eigenvalue*A->alpha; 00397 break; 00398 default: 00399 DSDPSETERR(1,"Vech Matrix not factored yet\n"); 00400 } 00401 00402 return 0; 00403 } 00404 00405 static int VechMatView(void* AA){ 00406 vechmat* A=(vechmat*)AA; 00407 int info,i=0,j,k,rank=0,ishift=A->ishift,nnz=A->nnzeros; 00408 const int *ind=A->ind; 00409 const double *val=A->val; 00410 for (k=0; k<nnz; k++){ 00411 getij(ind[k]-ishift,&i,&j); 00412 printf("Row: %d, Column: %d, Value: %10.8f \n",i,j,A->alpha*val[k]); 00413 } 00414 if (A->factored>0){ 00415 info=VechMatGetRank(AA,&rank,A->n);DSDPCHKERR(info); 00416 printf("Detected Rank: %d\n",rank); 00417 } 00418 return 0; 00419 } 00420 00421 00422 static struct DSDPDataMat_Ops vechmatops; 00423 static const char *datamatname="STANDARD VECH MATRIX"; 00424 00425 static int VechMatOpsInitialize(struct DSDPDataMat_Ops *sops){ 00426 int info; 00427 if (sops==NULL) return 0; 00428 info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info); 00429 sops->matvecvec=VechMatVecVec; 00430 sops->matdot=VechMatDot; 00431 sops->matfnorm2=VechMatFNorm2; 00432 sops->mataddrowmultiple=VechMatAddRowMultiple; 00433 sops->mataddallmultiple=VechMatAddMultiple; 00434 sops->matview=VechMatView; 00435 sops->matdestroy=VechMatDestroy; 00436 sops->matfactor2=VechMatFactor; 00437 sops->matgetrank=VechMatGetRank; 00438 sops->matgeteig=VechMatGetEig; 00439 sops->matrownz=VechMatGetRowNnz; 00440 sops->matnnz=VechMatCountNonzeros; 00441 sops->id=3; 00442 sops->matname=datamatname; 00443 return 0; 00444 } 00445 00458 #undef __FUNCT__ 00459 #define __FUNCT__ "DSDPGetVechMat" 00460 int DSDPGetVechMat(int n,int ishift,double alpha, const int ind[], const double val[],int nnz, struct DSDPDataMat_Ops**sops, void**smat){ 00461 int info,i,j,k,itmp,nn=n*(n+1)/2; 00462 double dtmp; 00463 vechmat* AA; 00464 DSDPFunctionBegin; 00465 for (k=0;k<nnz;++k){ 00466 itmp=ind[k]-ishift; 00467 if (itmp>=nn){ 00468 getij(itmp,&i,&j); 00469 /* 00470 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); 00471 */ 00472 DSDPSETERR3(2,"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn); 00473 } else if (itmp<0){ 00474 DSDPSETERR1(2,"Illegal index value: %d. Must be >= 0\n",itmp); 00475 } 00476 } 00477 for (k=0;k<nnz;++k) dtmp=val[k]; 00478 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info); 00479 AA->factored=0; 00480 AA->Eig=0; 00481 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info); 00482 if (sops){*sops=&vechmatops;} 00483 if (smat){*smat=(void*)AA;} 00484 DSDPFunctionReturn(0); 00485 } 00486 00487 00488 #undef __FUNCT__ 00489 #define __FUNCT__ "VechMatComputeEigs" 00490 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){ 00491 00492 int i,j,k,nsub,neigs,info,*iptr,*perm,*invp; 00493 long int *i2darray=(long int*)DD; 00494 int ownarray1=0,ownarray2=0,ownarray3=0; 00495 int ishift=AA->ishift,nonzeros=AA->nnzeros; 00496 const int *ind=AA->ind; 00497 const double *val=AA->val; 00498 double *dmatarray=ss1,*dworkarray=ss2,maxeig,eps=1.0e-12,eps2=1.0e-12; 00499 00500 iptr=iiptr; perm=iptr+n; invp=perm+n; 00501 /* These operations were done before calling this routine * / 00502 / * Integer arrays corresponding to rows with nonzeros and inverse map * / 00503 memset((void*)iiptr,0,3*n*sizeof(int)); 00504 00505 / * Find number of nonzeros in each row * / 00506 for (i=0,k=0; k<nonzeros; k++){ 00507 getij(ind[k],i,n,&i,&j); 00508 iptr[i]++; iptr[j]++; 00509 } 00510 */ 00511 /* Number of rows with a nonzero. Order the rows with nonzeros. */ 00512 for (nsub=0,i=0; i<n; i++){ 00513 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;} 00514 } 00515 00516 /* create a dense array in which to put numbers */ 00517 if (nsub*nsub>nn1){ 00518 DSDPCALLOC2(&dmatarray,double,(nsub*nsub),&info); DSDPCHKERR(info); 00519 ownarray1=1; 00520 } 00521 memset((void*)dmatarray,0,nsub*nsub*sizeof(double)); 00522 if (nsub*nsub>nn2){ 00523 DSDPCALLOC2(&dworkarray,double,(nsub*nsub),&info); DSDPCHKERR(info); 00524 ownarray2=1; 00525 } 00526 00527 if (nsub*nsub*sizeof(long int)>nn0*sizeof(double)){ 00528 DSDPCALLOC2(&i2darray,long int,(nsub*nsub),&info); DSDPCHKERR(info); 00529 ownarray3=1; 00530 } 00531 00532 for (i=0,k=0; k<nonzeros; k++){ 00533 getij(ind[k]-ishift,&i,&j); 00534 dmatarray[perm[i]*nsub+perm[j]] += val[k]; 00535 if (i!=j){ 00536 dmatarray[perm[j]*nsub+perm[i]] += val[k]; 00537 } 00538 } 00539 /* Call LAPACK to compute the eigenvalues */ 00540 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub, 00541 W,n,WORK,n1,iiptr+3*n,n2-3*n); 00542 if (info){ 00543 memset((void*)dmatarray,0,nsub*nsub*sizeof(double)); 00544 for (i=0,k=0; k<nonzeros; k++){ 00545 getij(ind[k]-ishift,&i,&j); 00546 dmatarray[perm[i]*nsub+perm[j]] += val[k]; 00547 if (i!=j){dmatarray[perm[j]*nsub+perm[i]] += val[k];} 00548 } 00549 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub, 00550 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info); 00551 } 00552 for (maxeig=0,i=0;i<nsub;i++){ 00553 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); } 00554 } 00555 memset((void*)iptr,0,nsub*sizeof(int)); 00556 /* Compute sparsity pattern for eigenvalue and eigenvector structures */ 00557 /* Count the nonzero eigenvalues */ 00558 for (neigs=0,k=0; k<nsub; k++){ 00559 if (fabs(W[k]) /* /maxeig */ > eps){ 00560 for (j=0;j<nsub;j++){ 00561 if (fabs(dmatarray[nsub*k+j]) >= eps2){iptr[neigs]++; 00562 } else { dmatarray[nsub*k+j]=0.0;} 00563 } 00564 neigs++; 00565 /* 00566 } else if (fabs(W[k])>1.0e-100){ 00567 printf("SKIPPING EIGENVALUE: %4.4e, max is : %4.4e\n",W[k],maxeig); 00568 */ 00569 } 00570 } 00571 00572 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info); 00573 DSDPLogInfo(0,49," Data Mat has %d eigenvalues: \n",neigs); 00574 /* Copy into structure */ 00575 for (neigs=0,i=0; i<nsub; i++){ 00576 if (fabs(W[i]) > eps){ 00577 info=EigMatSetEig(AA->Eig,neigs,W[i],invp,dmatarray+nsub*i,nsub,n);DSDPCHKERR(info); 00578 neigs++; 00579 } 00580 } 00581 00582 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);} 00583 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);} 00584 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);} 00585 return 0; 00586 } 00587