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