DSDP
|
00001 #include "dsdpsdp.h" 00002 #include "dsdpcone_impl.h" 00003 #include "dsdpsys.h" 00007 static int SDPConeOperationsInitialize(struct DSDPCone_Ops*); 00008 00009 static int KSDPConeSetup(void*,DSDPVec); 00010 static int KSDPConeSetup2(void*, DSDPVec, DSDPSchurMat); 00011 static int KSDPConeSize(void*,double*); 00012 static int KSDPConeSparsity(void*,int,int[],int[],int); 00013 static int KSDPConeComputeHessian(void*,double,DSDPSchurMat,DSDPVec,DSDPVec); 00014 static int KSDPConeComputeMaxStepLength(void*, DSDPVec, DSDPDualFactorMatrix, double *); 00015 static int KSDPConeComputeSS(void*, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *); 00016 static int KSDPConeComputeLogSDeterminant(void *, double *, double*); 00017 static int KSDPConeComputeXX(void*, double, DSDPVec,DSDPVec,DSDPVec,double*); 00018 static int KSDPConeDestroy(void*); 00019 00020 #undef __FUNCT__ 00021 #define __FUNCT__ "KSDPConeComputeHessian" 00022 static int KSDPConeComputeHessian( void *K, double mu, DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2){ 00023 int info; 00024 SDPCone sdpcone=(SDPCone)K; 00025 DSDPFunctionBegin; 00026 info=SDPConeComputeHessian(sdpcone,mu,M,vrhs1,vrhs2);DSDPCHKERR(info); 00027 DSDPFunctionReturn(0); 00028 } 00029 00030 #undef __FUNCT__ 00031 #define __FUNCT__ "KDPConeMultiply" 00032 static int KSDPConeMultiply( void *K, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout){ 00033 int kk,info; 00034 SDPCone sdpcone=(SDPCone)K; 00035 SDPConeValid(sdpcone); 00036 DSDPFunctionBegin; 00037 for (kk=0; kk<sdpcone->nblocks; kk++){ 00038 info=SDPConeMultiply(sdpcone,kk,mu,vrow,vin,vout);DSDPCHKBLOCKERR(kk,info); 00039 } 00040 DSDPFunctionReturn(0); 00041 } 00042 00043 #undef __FUNCT__ 00044 #define __FUNCT__ "KDPConeRHS " 00045 static int KSDPConeRHS( void *K, double mu, DSDPVec vrow, DSDPVec vrhs1, DSDPVec vrhs2){ 00046 int kk,info; 00047 SDPCone sdpcone=(SDPCone)K; 00048 DSDPFunctionBegin; 00049 SDPConeValid(sdpcone); 00050 for (kk=0; kk<sdpcone->nblocks; kk++){ 00051 if (sdpcone->blk[kk].n<1) continue; 00052 info=SDPConeComputeRHS(sdpcone,kk,mu,vrow,vrhs1,vrhs2); DSDPCHKBLOCKERR(kk,info); 00053 } 00054 DSDPFunctionReturn(0); 00055 } 00056 00057 00058 #undef __FUNCT__ 00059 #define __FUNCT__ "KSDPConeSetup" 00060 static int KSDPConeSetup(void* K, DSDPVec y){ 00061 int info; 00062 SDPCone sdpcone=(SDPCone)K; 00063 DSDPFunctionBegin; 00064 SDPConeValid(sdpcone); 00065 info=SDPConeSetup(sdpcone,y);DSDPCHKERR(info); 00066 DSDPFunctionReturn(0); 00067 } 00068 00069 #undef __FUNCT__ 00070 #define __FUNCT__ "KSDPConeSetup2" 00071 static int KSDPConeSetup2(void* K, DSDPVec y, DSDPSchurMat M){ 00072 int info; 00073 SDPCone sdpcone=(SDPCone)K; 00074 DSDPFunctionBegin; 00075 info=SDPConeSetup2(sdpcone,y,M); DSDPCHKERR(info); 00076 DSDPFunctionReturn(0); 00077 } 00078 00079 00080 #undef __FUNCT__ 00081 #define __FUNCT__ "KSDPConeDestroy" 00082 static int KSDPConeDestroy(void* K){ 00083 int info; 00084 SDPCone sdpcone=(SDPCone)K; 00085 DSDPFunctionBegin; 00086 SDPConeValid(sdpcone); 00087 info=SDPConeDestroy(sdpcone);DSDPCHKERR(info); 00088 DSDPFunctionReturn(0); 00089 } 00090 00091 00092 #undef __FUNCT__ 00093 #define __FUNCT__ "KSDPConeSize" 00094 static int KSDPConeSize(void* K,double *n){ 00095 SDPCone sdpcone=(SDPCone)K; 00096 DSDPFunctionBegin; 00097 SDPConeValid(sdpcone); 00098 *n=sdpcone->nn; 00099 DSDPFunctionReturn(0); 00100 } 00101 00102 #undef __FUNCT__ 00103 #define __FUNCT__ "KSDPConeSparsity" 00104 static int KSDPConeSparsity(void *K,int row, int *tnnz, int rnnz[], int m){ 00105 SDPCone sdpcone=(SDPCone)K; 00106 SDPblk *blk=sdpcone->blk; 00107 int info,j,kk; 00108 int nnzblocks=sdpcone->ATR.nnzblocks[row],*nzblocks=sdpcone->ATR.nzblocks[row]; 00109 DSDPFunctionBegin; 00110 SDPConeValid(sdpcone); 00111 for (j=0; j<nnzblocks; j++){ 00112 kk=nzblocks[j]; 00113 if (blk[kk].n<1) continue; 00114 info=DSDPBlockDataMarkNonzeroMatrices(&blk[kk].ADATA,rnnz);DSDPCHKBLOCKERR(kk,info); 00115 } 00116 DSDPFunctionReturn(0); 00117 } 00118 00119 #undef __FUNCT__ 00120 #define __FUNCT__ "KSDPConeComputeSS" 00121 static int KSDPConeComputeSS(void *K, DSDPVec Y, DSDPDualFactorMatrix flag, DSDPTruth *ispsdefinite){ 00122 int kk,info; 00123 SDPCone sdpcone=(SDPCone)K; 00124 SDPblk *blk=sdpcone->blk; 00125 DSDPTruth psdefinite; 00126 DSDPDualMat SS; 00127 00128 DSDPFunctionBegin; 00129 psdefinite = DSDP_TRUE; 00130 for (kk=sdpcone->nblocks-1; kk>=0 && psdefinite == DSDP_TRUE; kk--){ 00131 if (blk[kk].n<1) continue; 00132 if (flag==DUAL_FACTOR) SS=blk[kk].S; 00133 else SS=blk[kk].SS; 00134 00135 switch (sdpcone->optype){ 00136 default: 00137 info=SDPConeComputeSS(sdpcone,kk,Y,blk[kk].T);DSDPCHKBLOCKERR(kk,info); 00138 info=DSDPDualMatSetArray(SS,blk[kk].T); DSDPCHKBLOCKERR(kk,info); 00139 info=DSDPDualMatCholeskyFactor(SS,&psdefinite); DSDPCHKBLOCKERR(kk,info); 00140 if (psdefinite == DSDP_FALSE){ 00141 if (flag==DUAL_FACTOR){ 00142 DSDPLogInfo(0,2,"Dual SDP Block %2.0f not PSD\n",kk); 00143 } else { 00144 DSDPLogInfo(0,2,"Primal SDP Block %2.0f not PSD\n",kk); 00145 } 00146 } 00147 break; 00148 } 00149 } 00150 *ispsdefinite=psdefinite; 00151 if (flag==DUAL_FACTOR && psdefinite==DSDP_TRUE){ 00152 info=DSDPVecCopy(Y,sdpcone->YY);DSDPCHKERR(info); 00153 } 00154 DSDPFunctionReturn(0); 00155 } 00156 00157 #undef __FUNCT__ 00158 #define __FUNCT__ "KSDPConeInvertSS" 00159 static int KSDPConeInvertSS(void *K){ 00160 int kk,info; 00161 SDPCone sdpcone=(SDPCone)K; 00162 DSDPDualMat SS; 00163 00164 DSDPFunctionBegin; 00165 SDPConeValid(sdpcone); 00166 for (kk=0;kk<sdpcone->nblocks;kk++){ 00167 if (sdpcone->blk[kk].n<1) continue; 00168 SS=sdpcone->blk[kk].S; 00169 info=DSDPDualMatInvert(SS);DSDPCHKBLOCKERR(kk,info); 00170 } 00171 DSDPFunctionReturn(0); 00172 } 00173 00174 00175 00176 #undef __FUNCT__ 00177 #define __FUNCT__ "KSDPConeComputeMaxStepLength" 00178 static int KSDPConeComputeMaxStepLength(void *K, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){ 00179 int kk,info; 00180 double smaxstep,maxmaxstep=1.0e20; 00181 SDPCone sdpcone=(SDPCone)K; 00182 DSDPDualMat SS; 00183 SDPblk *blk=sdpcone->blk; 00184 DSDPDSMat DS; 00185 DSDPVMat T; 00186 00187 DSDPFunctionBegin; 00188 SDPConeValid(sdpcone); 00189 for (kk=0; kk<sdpcone->nblocks; kk++){ 00190 if (blk[kk].n<1) continue; 00191 if (flag==DUAL_FACTOR) SS=blk[kk].S; 00192 else SS=blk[kk].SS; 00193 DS=blk[kk].DS; T=blk[kk].T; 00194 00195 info=SDPConeComputeSS(sdpcone,kk,DY,T);DSDPCHKBLOCKERR(kk,info); 00196 info=DSDPDSMatSetArray(DS,T); DSDPCHKBLOCKERR(kk,info); 00197 00198 info=DSDPLanczosStepSize( &blk[kk].Lanczos,blk[kk].W,blk[kk].W2,SS,DS,&smaxstep );DSDPCHKBLOCKERR(kk,info); 00199 DSDPLogInfo(0,12,"Block %d, PD %d, maxstepsize: %4.4e\n",kk,flag,smaxstep); 00200 maxmaxstep=DSDPMin(smaxstep,maxmaxstep); 00201 } 00202 *maxsteplength=maxmaxstep; 00203 DSDPFunctionReturn(0); 00204 } 00205 00206 00207 00208 #undef __FUNCT__ 00209 #define __FUNCT__ "KSDPConeAddANorm2" 00210 static int KSDPConeAddANorm2(void *K, DSDPVec ANorm2){ 00211 int kk,info; 00212 SDPCone sdpcone=(SDPCone)K; 00213 SDPblk *blk=sdpcone->blk; 00214 00215 DSDPFunctionBegin; 00216 SDPConeValid(sdpcone); 00217 for (kk=0; kk<sdpcone->nblocks; kk++){ 00218 if (blk[kk].n<1) continue; 00219 info=DSDPBlockANorm2( &blk[kk].ADATA,ANorm2,blk[kk].n); DSDPCHKBLOCKERR(kk,info); 00220 } 00221 DSDPFunctionReturn(0); 00222 } 00223 00224 00225 00226 #undef __FUNCT__ 00227 #define __FUNCT__ "KSDPConeSetX" 00228 static int KSDPConeSetX(void *K, double mu, DSDPVec Y,DSDPVec DY){ 00229 SDPCone sdpcone=(SDPCone)K; 00230 int info; 00231 DSDPFunctionBegin; 00232 SDPConeValid(sdpcone); 00233 info=DSDPVecCopy(Y,sdpcone->YX);DSDPCHKERR(info); 00234 info=DSDPVecCopy(DY,sdpcone->DYX);DSDPCHKERR(info); 00235 sdpcone->xmakermu=mu; 00236 DSDPFunctionReturn(0); 00237 } 00238 00239 00240 #undef __FUNCT__ 00241 #define __FUNCT__ "KSDPConeComputeXX" 00242 static int KSDPConeComputeXX(void *K, double mu, DSDPVec Y,DSDPVec DY,DSDPVec AX,double* tracexs){ 00243 00244 SDPCone sdpcone=(SDPCone)K; 00245 int info,kk; 00246 double xnorm,trxs,xtrace; 00247 DSDPVMat X; 00248 00249 DSDPFunctionBegin; 00250 SDPConeValid(sdpcone); 00251 info=KSDPConeSetX(K,mu,Y,DY);DSDPCHKERR(info); 00252 for (kk=0; kk<sdpcone->nblocks; kk++){ 00253 if (sdpcone->blk[kk].n<1) continue; 00254 X=sdpcone->blk[kk].T; 00255 info=SDPConeComputeX3(sdpcone,kk,mu,Y,DY,X);DSDPCHKBLOCKERR(kk,info); 00256 info=SDPConeComputeXDot(sdpcone,kk,Y,X,AX,&xtrace,&xnorm,&trxs);DSDPCHKBLOCKERR(kk,info); 00257 *tracexs+=trxs; 00258 DSDPLogInfo(0,10,"SDP Block %d: norm(X): %4.2e, trace(X): %4.2e, trace(XS): %4.2e.\n",kk,xnorm,xtrace,trxs); 00259 } 00260 DSDPFunctionReturn(0); 00261 } 00262 00263 00264 #undef __FUNCT__ 00265 #define __FUNCT__ "KSDPConeComputeLogSDeterminant" 00266 static int KSDPConeComputeLogSDeterminant(void *K, double *logdetobj, double *logdet){ 00267 int kk,info; 00268 double dlogdet=0,dlogdet2=0,dd; 00269 SDPCone sdpcone=(SDPCone)K; 00270 SDPblk *blk=sdpcone->blk; 00271 00272 DSDPFunctionBegin; 00273 SDPConeValid(sdpcone); 00274 for (kk=0; kk<sdpcone->nblocks; kk++){ 00275 if (blk[kk].n<1) continue; 00276 info=DSDPDualMatLogDeterminant(blk[kk].S,&dd);DSDPCHKBLOCKERR(kk,info); 00277 dlogdet+=dd*blk[kk].gammamu; 00278 dlogdet2+=dd*blk[kk].bmu; 00279 } 00280 *logdet=dlogdet; 00281 *logdetobj=dlogdet2; 00282 DSDPFunctionReturn(0); 00283 } 00284 00285 00286 #undef __FUNCT__ 00287 #define __FUNCT__ "KSDPConeMonitor" 00288 int KSDPConeMonitor(void *K, int tag){ 00289 DSDPFunctionBegin; 00290 DSDPFunctionReturn(0); 00291 } 00292 00293 static struct DSDPCone_Ops kops; 00294 static const char *sdpconename ="SDP Cone"; 00295 00296 #undef __FUNCT__ 00297 #define __FUNCT__ "SDPConeOperationsInitialize" 00298 static int SDPConeOperationsInitialize(struct DSDPCone_Ops* coneops){ 00299 int info; 00300 if (coneops==NULL) return 0; 00301 info=DSDPConeOpsInitialize(coneops); DSDPCHKERR(info); 00302 coneops->conehessian=KSDPConeComputeHessian; 00303 coneops->conerhs=KSDPConeRHS; 00304 coneops->conesetup=KSDPConeSetup; 00305 coneops->conesetup2=KSDPConeSetup2; 00306 coneops->conedestroy=KSDPConeDestroy; 00307 coneops->conecomputes=KSDPConeComputeSS; 00308 coneops->coneinverts=KSDPConeInvertSS; 00309 coneops->conesetxmaker=KSDPConeSetX; 00310 coneops->conecomputex=KSDPConeComputeXX; 00311 coneops->conemaxsteplength=KSDPConeComputeMaxStepLength; 00312 coneops->conelogpotential=KSDPConeComputeLogSDeterminant; 00313 coneops->conesize=KSDPConeSize; 00314 coneops->conesparsity=KSDPConeSparsity; 00315 coneops->conehmultiplyadd=KSDPConeMultiply; 00316 coneops->coneanorm2=KSDPConeAddANorm2; 00317 coneops->conemonitor=KSDPConeMonitor; 00318 coneops->id=1; 00319 coneops->name=sdpconename; 00320 return 0; 00321 } 00322 00323 #undef __FUNCT__ 00324 #define __FUNCT__ "DSDPAddSDP" 00325 00331 int DSDPAddSDP(DSDP dsdp,SDPCone sdpcone){ 00332 int info; 00333 DSDPFunctionBegin; 00334 SDPConeValid(sdpcone); 00335 info=SDPConeOperationsInitialize(&kops); DSDPCHKERR(info); 00336 info=DSDPAddCone(dsdp,&kops,(void*)sdpcone); DSDPCHKERR(info); 00337 DSDPFunctionReturn(0); 00338 } 00339