DSDP
|
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