sdpconesetup.c

Go to the documentation of this file.
00001 #include "dsdpsdp.h"
00002 #include "dsdpsys.h"
00008 #undef __FUNCT__  
00009 #define __FUNCT__ "DSDPDataTransposeInitialize"
00010 
00015 int DSDPDataTransposeInitialize(DSDPDataTranspose *ATranspose){
00016   DSDPFunctionBegin;
00017   ATranspose->nnzblocks=0;
00018   ATranspose->nzblocks=0;
00019   ATranspose->idA=0;
00020   ATranspose->idAP=0;
00021   ATranspose->ttnzmat=0;
00022   ATranspose->nnzblocks=0;
00023   DSDPFunctionReturn(0);
00024 }
00025 
00026 #undef __FUNCT__
00027 #define __FUNCT__ "DSDPDataTransposeSetup"
00028 
00036 int DSDPDataTransposeSetup(DSDPDataTranspose *ATranspose, SDPblk *blk, int nblocks, int m){
00037   
00038   int i,ii,kk,vvar,info;
00039   int nnzmats,tnzmats=0;
00040   DSDPFunctionBegin;
00041 
00042   info=DSDPDataTransposeTakeDown(ATranspose);DSDPCHKERR(info);
00043   /* Determine sparsity pattern of SDP Data Matrices */
00044 
00045   DSDPCALLOC2(&ATranspose->nnzblocks,int,(m),&info);DSDPCHKERR(info);
00046   DSDPCALLOC2(&ATranspose->nzblocks,int*,(m),&info);DSDPCHKERR(info);
00047   DSDPCALLOC2(&ATranspose->idA,int*,(m),&info);DSDPCHKERR(info);
00048   ATranspose->m=m;
00049   for (i=0;i<m;i++){ ATranspose->nnzblocks[i]=0; }
00050   for (kk=0; kk<nblocks; kk++){
00051     info=DSDPBlockDataMarkNonzeroMatrices(&blk[kk].ADATA,ATranspose->nnzblocks);DSDPCHKERR(info);
00052   }
00053   for (tnzmats=0,i=0;i<m;i++){ tnzmats += ATranspose->nnzblocks[i];}
00054 
00055   DSDPCALLOC2(&ATranspose->ttnzmat,int,tnzmats,&info);DSDPCHKERR(info);
00056   ATranspose->nzblocks[0]=ATranspose->ttnzmat;
00057   for (i=1;i<m;i++){
00058     ATranspose->nzblocks[i]=ATranspose->nzblocks[i-1]+ATranspose->nnzblocks[i-1];
00059   }
00060 
00061   DSDPCALLOC2(&ATranspose->idAP,int,tnzmats,&info);DSDPCHKERR(info);
00062   ATranspose->idA[0]=ATranspose->idAP;
00063   for (i=1;i<m;i++){
00064     ATranspose->idA[i]=ATranspose->idA[i-1]+ATranspose->nnzblocks[i-1];
00065   }
00066 
00067   for (i=0;i<m;i++){ATranspose->nnzblocks[i]=0;}
00068   for (kk=0; kk<nblocks; kk++){
00069     info=DSDPBlockCountNonzeroMatrices(&blk[kk].ADATA,&nnzmats);DSDPCHKERR(info);
00070     for (i=0;i<nnzmats;i++){
00071       info=DSDPBlockGetMatrix(&blk[kk].ADATA,i,&ii,0,0);DSDPCHKERR(info);
00072       vvar=ATranspose->nnzblocks[ii];
00073       ATranspose->nzblocks[ii][vvar]=kk;
00074       ATranspose->idA[ii][vvar]=i;
00075       ATranspose->nnzblocks[ii]++;
00076     }
00077   }
00078 
00079   DSDPFunctionReturn(0);
00080 }
00081 
00082 #undef __FUNCT__  
00083 #define __FUNCT__ "DSDPDataTransposeTakeDown"
00084 
00089 int DSDPDataTransposeTakeDown(DSDPDataTranspose *ATranspose){
00090   int info;
00091   DSDPFunctionBegin;
00092   DSDPFREE(&ATranspose->ttnzmat,&info);DSDPCHKERR(info);
00093   DSDPFREE(&ATranspose->idAP,&info);DSDPCHKERR(info);
00094   DSDPFREE(&ATranspose->nzblocks,&info);DSDPCHKERR(info);
00095   DSDPFREE(&ATranspose->nnzblocks,&info);DSDPCHKERR(info);
00096   DSDPFREE(&ATranspose->idA,&info);DSDPCHKERR(info);
00097   info=DSDPDataTransposeInitialize(ATranspose);DSDPCHKERR(info);
00098   DSDPFunctionReturn(0);
00099 }
00100 
00101 #undef __FUNCT__  
00102 #define __FUNCT__ "DSDPCreateSDPCone"
00103 
00113 int DSDPCreateSDPCone(DSDP dsdp, int blocks, SDPCone* dspcone){
00114   int i,info;
00115   SDPCone sdpcone;
00116 
00117   DSDPFunctionBegin;
00118   DSDPCALLOC1(&sdpcone,struct SDPCone_C,&info);DSDPCHKERR(info);
00119   *dspcone=sdpcone;
00120   sdpcone->keyid=SDPCONEKEY;
00121   info=DSDPAddSDP(dsdp,sdpcone);DSDPCHKERR(info);
00122 
00123   info=DSDPGetNumberOfVariables(dsdp,&sdpcone->m);DSDPCHKERR(info);
00124   DSDPCALLOC2(&sdpcone->blk,SDPblk,blocks,&info); DSDPCHKERR(info);
00125   for (i=0;i<blocks; i++){
00126     info=DSDPBlockInitialize(&sdpcone->blk[i]); DSDPCHKBLOCKERR(i,info);
00127   }
00128 
00129   sdpcone->nblocks=blocks;
00130   sdpcone->optype=3;
00131   info=DSDPUseDefaultDualMatrix(sdpcone); DSDPCHKERR(info);
00132 
00133   sdpcone->nn=0;
00134   sdpcone->dsdp=dsdp;
00135   info=DSDPDataTransposeInitialize(&sdpcone->ATR); DSDPCHKERR(info);
00136   info=DSDPBlockEventZero();DSDPCHKERR(info);
00137   info=DSDPDualMatEventZero();DSDPCHKERR(info);
00138   info=DSDPVMatEventZero();DSDPCHKERR(info);
00139   DSDPFunctionReturn(0);
00140 }
00141 
00142 
00143 int DSDPCreateS(DSDPBlockData*,char,int,DSDPVec,DSDPVMat,SDPConeVec,SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
00144 
00145 #undef __FUNCT__ 
00146 #define __FUNCT__ "DSDPBlockSetup"
00147 
00154 int DSDPBlockSetup(SDPblk *blk, int blockj, DSDPVec WY){ 
00155   int n,info,trank,flag;
00156   DSDPFunctionBegin;
00157   /*  
00158   info=DSDPBlockTakeDown(blk);DSDPCHKERR(info);
00159   */
00160   n=blk->n;
00161   info=DSDPVMatExist(blk->T,&flag);DSDPCHKERR(info);
00162   if (flag==0){
00163     info=DSDPMakeVMat(blk->format,n,&blk->T);DSDPCHKERR(info);
00164   }
00165 
00166   info = DSDPIndexCreate(blk->n,&blk->IS);DSDPCHKERR(info);
00167   info = SDPConeVecCreate(blk->n,&blk->W);DSDPCHKERR(info);
00168   info = SDPConeVecDuplicate(blk->W,&blk->W2);DSDPCHKERR(info);
00169 
00170   /* Build Lanczos structures */
00171   info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,20); DSDPCHKERR(info);
00172   if (n>30){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,20); DSDPCHKERR(info);}
00173   if (n>300){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,40); DSDPCHKERR(info);}
00174   if (n>1000){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,50); DSDPCHKERR(info);}
00175   if (1){
00176     info=DSDPFastLanczosSetup(&blk->Lanczos,blk->W);DSDPCHKERR(info);
00177     DSDPLogInfo(0,19,"SDP Block %d using Fast Lanczos\n",blockj);
00178   } else {
00179     info=DSDPRobustLanczosSetup(&blk->Lanczos,blk->W);DSDPCHKERR(info);
00180     DSDPLogInfo(0,19,"SDP Block %d using Full Lanczos\n",blockj);
00181   }
00182 
00183   /* Eigenvalues and Eigenvectors */
00184   info=DSDPBlockFactorData(&blk->ADATA,blk->T,blk->W);DSDPCHKERR(info);
00185   info=DSDPBlockDataRank(&blk->ADATA,&trank,n);DSDPCHKERR(info);
00186 
00187   info=DSDPCreateS(&blk->ADATA,blk->format,trank,WY,blk->T,blk->W,blk->W2,&blk->S,&blk->SS,&blk->DS,0);DSDPCHKERR(info);
00188   
00189   DSDPFunctionReturn(0);
00190 }
00191 
00192 #undef __FUNCT__  
00193 #define __FUNCT__ "SDPConeBlockNNZ"
00194 int SDPConeBlockNNZ(SDPblk *blk,int m){
00195   int i,ii,n,info,nnz,nnzmats,tnnzmats,tnnz=0;
00196   double scl;
00197   DSDPDataMat AA;
00198   DSDPFunctionBegin;
00199   nnzmats=blk->ADATA.nnzmats;tnnzmats=nnzmats;
00200   n=blk->n;
00201 
00202   for (i=0;i<nnzmats;i++){
00203     info=DSDPBlockGetMatrix(&blk->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
00204     if (ii==0){tnnzmats--; continue;}
00205     if (ii==m-1){continue;}
00206     info = DSDPDataMatCountNonzeros(AA,&nnz,n); DSDPCHKERR(info);
00207     tnnz+= (nnz*(tnnzmats-i));
00208   }
00209   if (tnnzmats>1){ tnnz=tnnz/((tnnzmats)*(tnnzmats+1)/2); }
00210   if (tnnz<1)  tnnz = 1; 
00211   blk->nnz=tnnz;
00212   DSDPFunctionReturn(0);
00213 }
00214 
00215 #undef __FUNCT__  
00216 #define __FUNCT__ "SDPConeSetup2"
00217 
00224 int SDPConeSetup2(SDPCone sdpcone, DSDPVec yy0, DSDPSchurMat M){  
00225   int kk,n,m,info;
00226   double nn=0;
00227   SDPblk *blk;
00228   DSDPFunctionBegin;
00229   info=DSDPVecGetSize(yy0,&m);DSDPCHKERR(info);
00230   for (kk=0; kk<sdpcone->nblocks; kk++){
00231     blk=&sdpcone->blk[kk];
00232     n=blk->n;
00233     info=SDPConeBlockNNZ(blk,m);DSDPCHKERR(info);
00234     info=DSDPBlockSetup(blk,kk,sdpcone->Work);DSDPCHKERR(info);
00235     nn+=n*blk->gammamu;
00236   }
00237   sdpcone->nn=(int)nn;
00238   DSDPFunctionReturn(0);
00239 }
00240 
00241 #undef __FUNCT__  
00242 #define __FUNCT__ "SDPConeSetup"
00243 
00249 int SDPConeSetup(SDPCone sdpcone, DSDPVec yy0){
00250   int kk,n,m,info;
00251   DSDPFunctionBegin;
00252 
00253   info = DSDPVecGetSize(yy0,&m);DSDPCHKERR(info);
00254   if (m!=sdpcone->m+2){DSDPSETERR(8,"CHECK DIMENSION\n");}
00255   info = DSDPVecDuplicate(yy0,&sdpcone->Work);DSDPCHKERR(info);
00256   info = DSDPVecDuplicate(yy0,&sdpcone->Work2);DSDPCHKERR(info);
00257   info = DSDPVecDuplicate(yy0,&sdpcone->YY);DSDPCHKERR(info);
00258   info = DSDPVecDuplicate(yy0,&sdpcone->YX);DSDPCHKERR(info);
00259   info = DSDPVecDuplicate(yy0,&sdpcone->DYX);DSDPCHKERR(info);
00260   for (kk=0; kk<sdpcone->nblocks; kk++){
00261     n=sdpcone->blk[kk].n;
00262     info=SDPConeSetRIdentity(sdpcone,kk,n,1.0);DSDPCHKERR(info);
00263   }
00264   
00265   info=DSDPDataTransposeSetup(&sdpcone->ATR,sdpcone->blk,sdpcone->nblocks,m); DSDPCHKERR(info);
00266   info=DSDPBlockEventInitialize();DSDPCHKERR(info);
00267   info=DSDPDualMatEventInitialize();DSDPCHKERR(info);
00268   info=DSDPVMatEventInitialize();DSDPCHKERR(info);
00269   DSDPFunctionReturn(0);
00270 }
00271 
00272 #undef __FUNCT__ 
00273 #define __FUNCT__ "DSDPBlockInitialize"
00274 
00279 int DSDPBlockInitialize(SDPblk *blk){
00280   int info;
00281   DSDPFunctionBegin;
00282   blk->n=0;
00283   blk->format='N';
00284   blk->gammamu=1.0;
00285   blk->bmu=0.0;
00286   blk->nnz=10000000;
00287 
00288   info = DSDPDualMatInitialize(&blk->S); DSDPCHKERR(info);
00289   info = DSDPDualMatInitialize(&blk->SS); DSDPCHKERR(info);
00290   info = DSDPDSMatInitialize(&blk->DS); DSDPCHKERR(info);
00291   info = DSDPVMatInitialize(&blk->T); DSDPCHKERR(info);
00292   info = DSDPLanczosInitialize(&blk->Lanczos); DSDPCHKERR(info);
00293   info = DSDPBlockDataInitialize(&blk->ADATA); DSDPCHKERR(info);
00294   info = DSDPIndexInitialize(&blk->IS); DSDPCHKERR(info);
00295   DSDPFunctionReturn(0);
00296 }
00297   
00298 #undef __FUNCT__  
00299 #define __FUNCT__ "DSDPBlockTakeDown"
00300 
00305 int DSDPBlockTakeDown(SDPblk *blk){
00306   int info;
00307   DSDPFunctionBegin;
00308   if (!blk){DSDPFunctionReturn(0);}
00309   info=DSDPBlockTakeDownData(&blk->ADATA);DSDPCHKERR(info);
00310   info=SDPConeVecDestroy(&blk->W);DSDPCHKERR(info);
00311   info=SDPConeVecDestroy(&blk->W2);DSDPCHKERR(info);
00312   info=DSDPIndexDestroy(&blk->IS);DSDPCHKERR(info);
00313   info=DSDPLanczosDestroy(&blk->Lanczos);DSDPCHKERR(info);
00314   info=DSDPDualMatDestroy(&blk->SS);DSDPCHKERR(info);
00315   info=DSDPDualMatDestroy(&blk->S);DSDPCHKERR(info);
00316   info=DSDPDSMatDestroy(&blk->DS);DSDPCHKERR(info);
00317   info=DSDPVMatDestroy(&blk->T);DSDPCHKERR(info);
00318   DSDPFunctionReturn(0);
00319 }
00320 
00321 #undef __FUNCT__  
00322 #define __FUNCT__ "DSDPConeTakeDown"
00323 
00328 int DSDPConeTakeDown(SDPCone sdpcone){
00329   int blockj,info;
00330   DSDPFunctionBegin;
00331   for (blockj=0; blockj<sdpcone->nblocks; blockj++){
00332     info=DSDPBlockTakeDown(&sdpcone->blk[blockj]);DSDPCHKERR(info);
00333   }
00334   info=DSDPVecDestroy(&sdpcone->Work);DSDPCHKERR(info);
00335   info=DSDPVecDestroy(&sdpcone->Work2);DSDPCHKERR(info);
00336   info=DSDPVecDestroy(&sdpcone->YY);DSDPCHKERR(info);
00337   info=DSDPVecDestroy(&sdpcone->YX);DSDPCHKERR(info);
00338   info=DSDPVecDestroy(&sdpcone->DYX);DSDPCHKERR(info);
00339   info=DSDPDataTransposeTakeDown(&sdpcone->ATR);DSDPCHKERR(info);
00340   DSDPFunctionReturn(0);
00341 }
00342 
00343 #undef __FUNCT__  
00344 #define __FUNCT__ "SDPConeDestroy"
00345 
00350 int SDPConeDestroy(SDPCone sdpcone){
00351   int blockj,info;
00352   DSDPFunctionBegin;
00353   info=DSDPConeTakeDown(sdpcone);DSDPCHKERR(info);
00354   for (blockj=0; blockj<sdpcone->nblocks; blockj++){
00355     info=DSDPBlockDataDestroy(&sdpcone->blk[blockj].ADATA);DSDPCHKERR(info);
00356   }
00357   DSDPFREE(&sdpcone->blk,&info);DSDPCHKERR(info);
00358   DSDPFREE(&sdpcone,&info);DSDPCHKERR(info);
00359   info=DSDPBlockEventZero();DSDPCHKERR(info);
00360   info=DSDPDualMatEventZero();DSDPCHKERR(info);
00361   info=DSDPVMatEventZero();DSDPCHKERR(info);
00362   DSDPFunctionReturn(0);
00363 }
00364 

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