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