dsdpblock.c

Go to the documentation of this file.
00001 #include "dsdpsdp.h"
00002 #include "dsdpsys.h"
00003 
00004 static int sdpvecvecevent=0,sdpdotevent=0;
00010 #undef __FUNCT__
00011 #define __FUNCT__ "DSDPBlockASum"
00012 
00020 int DSDPBlockASum(DSDPBlockData *ADATA, double aa, DSDPVec Yk, DSDPVMat XX){
00021 
00022   double *xx,ytmp,scl=ADATA->scl;
00023   int    ii,vari,n,nn,info;
00024 
00025   DSDPFunctionBegin;
00026   info=DSDPVMatGetSize(XX, &n); DSDPCHKERR(info);
00027   info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info);
00028   for (ii=0;ii<ADATA->nnzmats;ii++){
00029     vari=ADATA->nzmat[ii];
00030     info=DSDPVecGetElement(Yk,vari,&ytmp);DSDPCHKVARERR(vari,info);
00031     if (ytmp==0) continue;
00032     info = DSDPDataMatAddMultiple(ADATA->A[ii], -aa*scl*ytmp, xx,nn,n); DSDPCHKVARERR(vari,info);
00033   }
00034   info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info);
00035   DSDPFunctionReturn(0);
00036 }
00037 
00038 #undef __FUNCT__
00039 #define __FUNCT__ "DSDPBlockADot"
00040 
00049 int DSDPBlockADot(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, DSDPVMat X, DSDPVec AX){
00050 
00051   int    ii,vari,n,nn,info;
00052   double *x,sum=0,aalpha=0,scl=ADATA->scl;
00053 
00054   DSDPFunctionBegin;
00055   DSDPEventLogBegin(sdpdotevent);
00056   info=DSDPVMatScaleDiagonal(X,0.5); DSDPCHKERR(info);
00057   info=DSDPVMatGetSize(X, &n); DSDPCHKERR(info);
00058   info=DSDPVMatGetArray(X, &x, &nn); DSDPCHKERR(info);
00059   for (ii=0;ii<ADATA->nnzmats; ii++){  /* Matrix Entries */
00060     vari=ADATA->nzmat[ii];
00061     info=DSDPVecGetElement(Alpha,vari,&aalpha);DSDPCHKVARERR(vari,info);
00062     if (aalpha==0.0) continue;
00063     info=DSDPDataMatDot(ADATA->A[ii],x,nn,n,&sum);DSDPCHKVARERR(vari,info);
00064     info=DSDPVecAddElement(AX,vari,aa*aalpha*sum*scl);DSDPCHKVARERR(vari,info);
00065   }
00066   info=DSDPVMatRestoreArray(X, &x, &nn); DSDPCHKERR(info);
00067   info=DSDPVMatScaleDiagonal(X,2.0); DSDPCHKERR(info);
00068   DSDPEventLogEnd(sdpdotevent);
00069   DSDPFunctionReturn(0);
00070 }
00071 
00072 #undef __FUNCT__
00073 #define __FUNCT__ "DSDPBlockvAv"
00074 
00084 int DSDPBlockvAv(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, SDPConeVec V, DSDPVec VAV){
00085 
00086   int    ii,vari,info;
00087   double sum=0,aalpha=0,scl=ADATA->scl;
00088 
00089   DSDPFunctionBegin;
00090   DSDPEventLogBegin(sdpvecvecevent);
00091   if (aa==0){DSDPFunctionReturn(0);}
00092   for (ii=0;ii<ADATA->nnzmats; ii++){  /* Matrix Entries */
00093     vari=ADATA->nzmat[ii];
00094     info=DSDPVecGetElement(Alpha,vari,&aalpha);DSDPCHKVARERR(vari,info);
00095     if (aalpha==0.0) continue;
00096     info=DSDPDataMatVecVec(ADATA->A[ii],V,&sum);DSDPCHKVARERR(vari,info);
00097     info=DSDPVecAddElement(VAV,vari,aa*aalpha*sum*scl);DSDPCHKVARERR(vari,info);
00098   }
00099   DSDPEventLogEnd(sdpvecvecevent);
00100   DSDPFunctionReturn(0);
00101 }
00102 
00103 #undef __FUNCT__
00104 #define __FUNCT__ "DSDPBlockFactorData"
00105 
00113 int DSDPBlockFactorData(DSDPBlockData *ADATA, 
00114                         DSDPVMat X, SDPConeVec W){
00115   
00116   int    ii,vari,n,nn,info,*iwork3n,i13,n26;
00117   double *x,*dwork3n;
00118   DSDPDataMat AA;
00119 
00120   DSDPFunctionBegin;
00121   info=DSDPVMatGetSize(X, &n); DSDPCHKERR(info);
00122   i13=13*n+1;n26=26*n+1;
00123   DSDPCALLOC2(&dwork3n,double,n26,&info);DSDPCHKERR(info);
00124   DSDPCALLOC2(&iwork3n,int,i13,&info);DSDPCHKERR(info);
00125   info=DSDPVMatGetArray(X, &x, &nn); DSDPCHKERR(info);
00126   for (ii=0;ii<ADATA->nnzmats; ii++){  /* Matrix Entries */
00127     info=DSDPBlockGetMatrix(ADATA,ii,&vari,0,&AA);DSDPCHKVARERR(vari,info);
00128     DSDPLogInfo(0,39,"SDP Data Mat Setup: %d\n",vari);
00129     if (vari==0) continue;
00130     info=DSDPDataMatFactor(AA,W,x,nn,dwork3n,n26,iwork3n,i13); DSDPCHKVARERR(vari,info);
00131   }
00132   info=DSDPVMatRestoreArray(X, &x, &nn); DSDPCHKERR(info);
00133   DSDPFREE(&dwork3n,&info);DSDPCHKERR(info);
00134   DSDPFREE(&iwork3n,&info);DSDPCHKERR(info);
00135   DSDPFunctionReturn(0);
00136 }
00137 
00138 
00139 #undef __FUNCT__
00140 #define __FUNCT__ "DSDPBlockEventZero"
00141 int DSDPBlockEventZero(void){
00142   DSDPFunctionBegin;
00143   sdpvecvecevent=0;sdpdotevent=0;
00144   DSDPFunctionReturn(0);
00145 }
00146 
00147 #undef __FUNCT__
00148 #define __FUNCT__ "DSDPBlockEventInitialize"
00149 int DSDPBlockEventInitialize(void){
00150   DSDPFunctionBegin;
00151   if (sdpvecvecevent==0){DSDPEventLogRegister("SDP VecMatVec",&sdpvecvecevent);}
00152   if (sdpdotevent==0){DSDPEventLogRegister("SDP Dot",&sdpdotevent);}
00153   DSDPFunctionReturn(0);
00154 }
00155 
00156 #undef __FUNCT__
00157 #define __FUNCT__ "DSDPBlockDataInitialize"
00158 
00163 int DSDPBlockDataInitialize(DSDPBlockData *ADATA){
00164   DSDPFunctionBegin;
00165   ADATA->nnzmats=0;
00166   ADATA->maxnnzmats=0;
00167   ADATA->nzmat=0;
00168   ADATA->A=0;
00169   ADATA->r=1.0;
00170   ADATA->scl=1.0;
00171   /*  ADATA->n=0; */
00172   DSDPFunctionReturn(0);
00173 }
00174 
00175 #undef __FUNCT__
00176 #define __FUNCT__ "DSDPBlockTakeDownData"
00177 
00182 int DSDPBlockTakeDownData(DSDPBlockData *ADATA){
00183   DSDPFunctionBegin;
00184   DSDPFunctionReturn(0);
00185 }
00186 
00187 
00188 #undef __FUNCT__
00189 #define __FUNCT__ "DSDPBlockDataDestroy"
00190 
00195 int DSDPBlockDataDestroy(DSDPBlockData *ADATA){
00196   int ii,vari,info;
00197   DSDPFunctionBegin;
00198   if (!ADATA){DSDPFunctionReturn(0);}
00199   DSDPLogInfo(0,18,"Destroying All Existing Data Matrices \n");
00200   for (ii=0; ii<ADATA->nnzmats; ii++){
00201     vari=ADATA->nzmat[ii];
00202     info = DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKVARERR(vari,info);
00203     ADATA->nzmat[ii]=0;
00204   }
00205   ADATA->nnzmats=0;
00206   info=DSDPBlockTakeDownData(ADATA);DSDPCHKERR(info);
00207   DSDPFREE(&ADATA->nzmat,&info);DSDPCHKERR(info);
00208   DSDPFREE(&ADATA->A,&info);DSDPCHKERR(info);
00209   info=DSDPBlockDataInitialize(ADATA);DSDPCHKERR(info);
00210   DSDPFunctionReturn(0);
00211 }
00212 
00213 #undef __FUNCT__  
00214 #define __FUNCT__ "DSDPBlockDataAllocate"
00215 
00221 int DSDPBlockDataAllocate(DSDPBlockData *ADATA, int nnz){
00222   int j,info,*nzmat;
00223   DSDPDataMat *A;
00224   DSDPFunctionBegin;
00225   if (!ADATA){DSDPFunctionReturn(0);}
00226   if (nnz<=ADATA->maxnnzmats){DSDPFunctionReturn(0);}
00227   DSDPLogInfo(0,18,"REALLOCATING SPACE FOR %d SDP BLOCK MATRICES! Previously allocated: %d \n",nnz,ADATA->maxnnzmats);
00228   DSDPCALLOC2(&A,struct DSDPDataMat_C,nnz,&info);DSDPCHKERR(info);
00229   DSDPCALLOC2(&nzmat,int,nnz,&info);DSDPCHKERR(info);
00230   for (j=0;j<nnz;j++){nzmat[j]=0;}
00231   for (j=0;j<nnz;j++){info = DSDPDataMatInitialize(&A[j]);DSDPCHKERR(info);}
00232   if (ADATA->maxnnzmats>0){
00233     for (j=0;j<ADATA->nnzmats;j++){nzmat[j]=ADATA->nzmat[j];}
00234     for (j=0;j<ADATA->nnzmats;j++){A[j]=ADATA->A[j];}
00235     DSDPFREE(&ADATA->A,&info);DSDPCHKERR(info);
00236     DSDPFREE(&ADATA->nzmat,&info);DSDPCHKERR(info);
00237   } else {
00238     ADATA->nnzmats=0;
00239   }
00240   ADATA->maxnnzmats=nnz;
00241   ADATA->nzmat=nzmat;
00242   ADATA->A=A;
00243   DSDPFunctionReturn(0);
00244 }
00245 
00246 #undef __FUNCT__
00247 #define __FUNCT__ "DSDPBlockDataMarkNonzeroMatrices"
00248 
00254 int DSDPBlockDataMarkNonzeroMatrices(DSDPBlockData *ADATA,int*annz){
00255   int i;
00256   DSDPFunctionBegin;
00257   for (i=0; i<ADATA->nnzmats; i++){
00258     annz[ADATA->nzmat[i]]++;
00259   }
00260   DSDPFunctionReturn(0);
00261 }
00262 
00263 #undef __FUNCT__
00264 #define __FUNCT__ "DSDPBlockCountNonzerosMatrices"
00265 
00272 int DSDPBlockCountNonzeroMatrices(DSDPBlockData *ADATA,int*nzmats){
00273   DSDPFunctionBegin;
00274   *nzmats=ADATA->nnzmats;
00275   DSDPFunctionReturn(0);
00276 }
00277 
00278 #undef __FUNCT__
00279 #define __FUNCT__ "DSDPBlockDataRank"
00280 int DSDPBlockDataRank(DSDPBlockData *ADATA, int *trank, int n){
00281   int    ii,vari,info,ri,r2=0;
00282   DSDPDataMat AA;
00283 
00284   DSDPFunctionBegin;
00285   for (ii=0;ii<ADATA->nnzmats;ii++){
00286     info=DSDPBlockGetMatrix(ADATA,ii,&vari,0,&AA);DSDPCHKVARERR(vari,info);
00287     if (vari==0) continue;
00288     info=DSDPDataMatGetRank(AA,&ri,n); DSDPCHKVARERR(vari,info);
00289     r2+=ri;
00290   }
00291   *trank=r2;
00292   DSDPFunctionReturn(0);
00293 }
00294 
00295 #undef __FUNCT__
00296 #define __FUNCT__ "DSDPBlockGetMatrix"
00297 
00307 int DSDPBlockGetMatrix(DSDPBlockData *ADATA,int id, int *vari, double *scl, DSDPDataMat *A){
00308   DSDPFunctionBegin;
00309   if (id>=0 && id < ADATA->nnzmats){
00310     if (vari) *vari=ADATA->nzmat[id];
00311     if (scl) *scl=ADATA->scl;
00312     if (A) *A=ADATA->A[id];
00313   } else {
00314     DSDPSETERR2(2,"Invalid Matrix request.  0 <= %d < %d\n",id,ADATA->nnzmats);
00315   }
00316   DSDPFunctionReturn(0);
00317 }
00318 
00319 #undef __FUNCT__
00320 #define __FUNCT__ "DSDPBlockDataRowSparsity"
00321 
00330 int DSDPBlockDataRowSparsity(DSDPBlockData *ADATA,int row, int ai[], int rnnz[],int n){
00331   int info,i,vari,rn;
00332   DSDPFunctionBegin;
00333   if (ai){
00334     for (i=0; i<ADATA->nnzmats; i++){
00335       vari=ADATA->nzmat[i];
00336       if (ai[vari]==0){continue;}
00337       info=DSDPDataMatGetRowNonzeros(ADATA->A[i],row, n, rnnz, &rn); DSDPCHKVARERR(vari,info);
00338     }
00339   }
00340   DSDPFunctionReturn(0);
00341 }
00342 
00343 #undef __FUNCT__  
00344 #define __FUNCT__ "DSDPBlockRemoveDataMatrix"
00345 
00351 int DSDPBlockRemoveDataMatrix(DSDPBlockData *ADATA,int vari){ 
00352   int info,ii,k;
00353   DSDPFunctionBegin;
00354   for (ii=0;ii<ADATA->nnzmats;ii++){
00355     if (ADATA->nzmat[ii]==vari){
00356       info=DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKVARERR(vari,info);
00357       info=DSDPSetDataMatZero(&ADATA->A[ii]);DSDPCHKVARERR(vari,info);
00358       for (k=ii;k<ADATA->nnzmats;k++){
00359         ADATA->A[k]=ADATA->A[k+1];
00360         ADATA->nzmat[k]=ADATA->nzmat[k+1];
00361       }
00362       ADATA->nnzmats--;
00363       info=DSDPSetDataMatZero(&ADATA->A[ADATA->nnzmats]);DSDPCHKERR(info);
00364       DSDPFunctionReturn(0);
00365     }
00366   }
00367   DSDPFunctionReturn(0);
00368 }
00369 
00370 #undef __FUNCT__  
00371 #define __FUNCT__ "DSDPBlockAddDataMatrix"
00372 
00381 int DSDPBlockAddDataMatrix(DSDPBlockData *ADATA,int vari, struct DSDPDataMat_Ops* dsdpdataops, void* data){ 
00382   int info,ii;
00383   DSDPFunctionBegin;
00384   if (ADATA->nnzmats>=ADATA->maxnnzmats){
00385     info=DSDPBlockDataAllocate(ADATA,2*ADATA->maxnnzmats+7);DSDPCHKERR(info);
00386   }
00387   ii=ADATA->nnzmats;
00388   info=DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKERR(info);
00389   info=DSDPDataMatSetData(&ADATA->A[ii], dsdpdataops, data);DSDPCHKVARERR(vari,info);
00390   ADATA->nzmat[ii]=vari;
00391   ADATA->nnzmats++;
00392   DSDPFunctionReturn(0);
00393 }
00394 
00395 #undef __FUNCT__  
00396 #define __FUNCT__ "DSDPBlockSetDataMatrix"
00397 
00406 int DSDPBlockSetDataMatrix(DSDPBlockData *ADATA,int vari, struct DSDPDataMat_Ops* dsdpdataops, void* data){ 
00407   int info;
00408   DSDPFunctionBegin;
00409   info=DSDPBlockRemoveDataMatrix(ADATA,vari);DSDPCHKERR(info);
00410   info=DSDPBlockAddDataMatrix(ADATA,vari,dsdpdataops,data);DSDPCHKERR(info);
00411   DSDPFunctionReturn(0);
00412 }
00413 
00414 #undef __FUNCT__
00415 #define __FUNCT__ "DSDPBlockNorm2"
00416 int DSDPBlockNorm2(DSDPBlockData *ADATA, int n){
00417   double fn2,tt=0;
00418   int    ii,info;
00419   DSDPFunctionBegin;
00420   for (ii=0;ii<ADATA->nnzmats;ii++){
00421     info=DSDPDataMatFNorm2(ADATA->A[ii],n,&fn2); DSDPCHKERR(info);
00422     tt+=fn2;
00423   }
00424   DSDPFunctionReturn(0);
00425 }
00426 
00427 #undef __FUNCT__
00428 #define __FUNCT__ "DSDPBlockANorm2"
00429 int DSDPBlockANorm2(DSDPBlockData *ADATA, DSDPVec ANorm2, int n){
00430 
00431   double fn2,scl=ADATA->scl;
00432   int    ii,vari,info;
00433 
00434   DSDPFunctionBegin;
00435   info=DSDPBlockNorm2(ADATA,n);DSDPCHKERR(info);
00436   scl=ADATA->scl;
00437   for (ii=0;ii<ADATA->nnzmats;ii++){
00438     vari=ADATA->nzmat[ii];
00439     info=DSDPDataMatFNorm2(ADATA->A[ii],n,&fn2); DSDPCHKVARERR(vari,info);
00440     info=DSDPVecAddElement(ANorm2,vari,fn2*scl);DSDPCHKVARERR(vari,info);
00441   }
00442   DSDPFunctionReturn(0);
00443 }
00444 
00445 
00446 #undef __FUNCT__
00447 #define __FUNCT__ "DSDPBlockView"
00448 
00454 int DSDPBlockView(DSDPBlockData *ADATA){
00455   int    ii,kk;
00456 
00457   DSDPFunctionBegin;
00458   for (ii=0;ii<ADATA->nnzmats;ii++){
00459     kk=ADATA->nzmat[ii];
00460     if (kk==0){ printf("+ C\n");}
00461     else { printf(" - A[%d] y%d\n",kk,kk);}
00462   }
00463   printf(" = S >= 0\n");
00464   DSDPFunctionReturn(0);
00465 }
00466 #undef __FUNCT__
00467 #define __FUNCT__ "DSDPBlockView2"
00468 
00474 int DSDPBlockView2(DSDPBlockData *ADATA){
00475   int    ii,vari,info;
00476 
00477   DSDPFunctionBegin;
00478   for (ii=0;ii<ADATA->nnzmats;ii++){
00479     vari=ADATA->nzmat[ii];
00480     printf("A[%d] y%d \n",vari,vari);
00481     info=DSDPDataMatView(ADATA->A[ii]); DSDPCHKERR(info);
00482   }
00483   DSDPFunctionReturn(0);
00484 }
00485 
00486 
00487 #undef __FUNCT__
00488 #define __FUNCT__ "DSDPDataMatCheck"
00489 
00498 int DSDPDataMatCheck(DSDPDataMat AA, SDPConeVec W, DSDPIndex IS, DSDPVMat XX){
00499 
00500   double *xx,ack,vAv=0,esum=0,sum,eignorm,fnorm22,dnorm,scl=1;
00501   int    k,n,nn,rank,info;
00502 
00503   DSDPFunctionBegin;
00504   info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info);
00505 
00506   info=DSDPVMatZeroEntries(XX);DSDPCHKERR(info);
00507   info=DSDPDataMatGetRank(AA,&rank,n);DSDPCHKERR(info);
00508   for (k=0; k<rank; k++){
00509     info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKERR(info);
00510     info=SDPConeVecDot(W,W,&eignorm);DSDPCHKERR(info);
00511     info=DSDPVMatAddOuterProduct(XX,scl*ack,W);DSDPCHKERR(info);
00512     info=DSDPDataMatVecVec(AA,W,&sum);DSDPCHKERR(info);
00513     vAv+=ack*ack*eignorm*eignorm*scl;
00514   }
00515   info=DSDPDataMatFNorm2(AA,n,&fnorm22); DSDPCHKERR(info);
00516   
00517   info=DSDPVMatScaleDiagonal(XX,0.5); DSDPCHKERR(info);
00518   info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info);
00519   info=DSDPDataMatDot(AA,xx,nn,n,&esum);DSDPCHKERR(info);
00520   info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info);
00521   info=DSDPVMatScaleDiagonal(XX,2.0); DSDPCHKERR(info);
00522   
00523   info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info);
00524   info=DSDPDataMatAddMultiple(AA,-1.0,xx,nn,n); DSDPCHKERR(info);
00525   info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info);
00526   if (0==1){info=DSDPVMatView(XX);DSDPCHKERR(info);}
00527   info=DSDPVMatNormF2(XX,&dnorm); DSDPCHKERR(info);
00528   printf("  %4.4e, %4.4e  %4.4e\n",esum,vAv,fnorm22);
00529   printf("  error1: %4.4e, error2: %4.4e,  error3: %4.4e\n",sqrt(dnorm),fabs(esum-vAv),fabs(fnorm22-vAv));
00530   if (dnorm>1)  printf("Check Add or eigs\n");
00531   if (fabs(esum-vAv) > 1.0) printf("Check vAv \n");
00532   if (fabs(fnorm22-vAv) > 1.0) printf("Check fnorm22\n");
00533   
00534   DSDPFunctionReturn(0);
00535 }
00536 

Generated on Wed Nov 11 20:41:02 2009 for DSDP by  doxygen 1.6.1