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
00042
00043
00044
00045
00046
00047
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
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
00189
00190
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
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
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
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
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512 for (nsub=0,i=0; i<n; i++){
00513 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
00514 }
00515
00516
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
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
00557
00558 for (neigs=0,k=0; k<nsub; k++){
00559 if (fabs(W[k]) > 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
00567
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
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