DSDP
|
00001 #include "dsdpcone_impl.h" 00002 #include "dsdpsys.h" 00003 #include "dsdpbasictypes.h" 00009 struct RDCone{ 00010 double primalr,dualr,x,logr; 00011 DSDPPenalty UsePenalty; 00012 DSDP dsdp; /* Really only need the Penalty flag, which is here. */ 00013 }; 00014 00015 typedef struct RDCone RCone; 00016 00017 #undef __FUNCT__ 00018 #define __FUNCT__ "DSDPRHessian" 00019 static int DSDPRHessian( void *dspcone, double mu, DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2){ 00020 RCone *K=(RCone*)dspcone; 00021 int info,m; 00022 double rr,sr,rssr; 00023 DSDPFunctionBegin; 00024 if (K->dualr ){ 00025 info=DSDPVecGetSize(vrhs2,&m);DSDPCHKERR(info); 00026 info=DSDPSchurMatVariableCompute(M,m-1,&rr);DSDPCHKERR(info); 00027 if (rr){ 00028 sr=-mu*rr/(K->dualr); 00029 rssr=mu*rr/(K->dualr*K->dualr); 00030 info=DSDPVecAddR(vrhs2,sr);DSDPCHKERR(info); 00031 info=DSDPSchurMatAddDiagonalElement(M,m-1,rssr);DSDPCHKERR(info); 00032 } 00033 } 00034 DSDPFunctionReturn(0); 00035 } 00036 00037 00038 #undef __FUNCT__ 00039 #define __FUNCT__ "DSDPRHS" 00040 static int DSDPRHS( void *dspcone, double mu, DSDPVec vrow,DSDPVec vrhs1,DSDPVec vrhs2){ 00041 RCone *K=(RCone*)dspcone; 00042 int info; 00043 double rr,sr; 00044 DSDPFunctionBegin; 00045 if (K->dualr ){ 00046 sr=-mu/(K->dualr); 00047 info=DSDPVecGetR(vrow,&rr);DSDPCHKERR(info); 00048 info=DSDPVecAddR(vrhs2,rr*sr);DSDPCHKERR(info); 00049 } 00050 DSDPFunctionReturn(0); 00051 } 00052 00053 00054 #undef __FUNCT__ 00055 #define __FUNCT__ "DSDPSetupRCone" 00056 static int DSDPSetupRCone(void* dspcone,DSDPVec y){ 00057 DSDPFunctionBegin; 00058 DSDPFunctionReturn(0); 00059 } 00060 00061 #undef __FUNCT__ 00062 #define __FUNCT__ "DSDPSetupRCone2" 00063 static int DSDPSetupRCone2(void* dspcone, DSDPVec y, DSDPSchurMat M){ 00064 DSDPFunctionBegin; 00065 DSDPFunctionReturn(0); 00066 } 00067 00068 00069 #undef __FUNCT__ 00070 #define __FUNCT__ "DSDPDestroyRCone" 00071 static int DSDPDestroyRCone(void* dspcone){ 00072 int info; 00073 DSDPFunctionBegin; 00074 DSDPFREE(&dspcone,&info);DSDPCHKERR(info); 00075 DSDPFunctionReturn(0); 00076 } 00077 00078 00079 #undef __FUNCT__ 00080 #define __FUNCT__ "DSDPRSize" 00081 static int DSDPRSize(void*dspcone,double *n){ 00082 RCone *K=(RCone*)dspcone; 00083 DSDPFunctionBegin; 00084 if (K->dualr){*n=1;} 00085 else {*n=0;} 00086 DSDPFunctionReturn(0); 00087 } 00088 00089 #undef __FUNCT__ 00090 #define __FUNCT__ "DSDPRSparsity" 00091 static int DSDPRSparsity(void*dspcone,int row, int *tnnz, int rnnz[], int m){ 00092 DSDPFunctionBegin; 00093 *tnnz=0; 00094 DSDPFunctionReturn(0); 00095 } 00096 00097 00098 #undef __FUNCT__ 00099 #define __FUNCT__ "DSDPComputeRS" 00100 static int DSDPComputeRS(void *dspcone, DSDPVec Y, DSDPDualFactorMatrix flag, DSDPTruth *ispsdefinite){ 00101 RCone *K=(RCone*)dspcone; 00102 int info; 00103 double rbeta; 00104 DSDPFunctionBegin; 00105 info=DSDPVecGetR(Y,&rbeta); DSDPCHKERR(info); 00106 if (K->UsePenalty==DSDPAlways){ 00107 if (rbeta<0){ *ispsdefinite=DSDP_TRUE; } else { *ispsdefinite=DSDP_FALSE;} 00108 } else { 00109 if (rbeta>0) rbeta=0; 00110 *ispsdefinite=DSDP_TRUE; 00111 } 00112 if (flag==DUAL_FACTOR){ K->dualr=rbeta; } else { K->primalr=rbeta;} 00113 DSDPFunctionReturn(0); 00114 } 00115 #undef __FUNCT__ 00116 #define __FUNCT__ "DSDPInvertRS" 00117 static int DSDPInvertRS(void *dspcone){ 00118 DSDPFunctionBegin; 00119 DSDPFunctionReturn(0); 00120 } 00121 00122 00123 00124 #undef __FUNCT__ 00125 #define __FUNCT__ "DSDPComputeRStepLength" 00126 static int DSDPComputeRStepLength(void *dspcone, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){ 00127 RCone *K=(RCone*)dspcone; 00128 double r,rbeta,msteplength=1.0e100,rt=1.0e30; 00129 int info; 00130 DSDPFunctionBegin; 00131 00132 info=DSDPVecGetR(DY,&rbeta); DSDPCHKERR(info); 00133 if (flag==DUAL_FACTOR){ r=K->dualr; } else { r=K->primalr;} 00134 if (r * rbeta<0) rt=-r/rbeta; 00135 00136 if (K->UsePenalty==DSDPAlways){msteplength=rt;} 00137 else if (flag==PRIMAL_FACTOR){ msteplength=rt;} 00138 else if (flag==DUAL_FACTOR){msteplength=rt/0.94;} 00139 00140 *maxsteplength=msteplength; 00141 DSDPFunctionReturn(0); 00142 } 00143 00144 #undef __FUNCT__ 00145 #define __FUNCT__ "DSDPRX" 00146 static int DSDPRX( void *dspcone, double mu, DSDPVec y, DSDPVec dy, DSDPVec AX,double *tracexs){ 00147 RCone *K=(RCone*)dspcone; 00148 int info; 00149 double rr,dr,trxs,r=K->dualr; 00150 DSDPFunctionBegin; 00151 00152 info=DSDPVecGetR(y,&rr); DSDPCHKERR(info); 00153 info=DSDPVecGetR(dy,&dr); DSDPCHKERR(info); 00154 if (K->dualr){ 00155 r=-1.0/rr; 00156 K->x=mu*(r-r*dr*r); 00157 trxs=K->x/r; 00158 DSDPLogInfo(0,2,"RESIDUAL X (Minimum Penalty Parameter): %4.4e, Trace(XS): %4.4e\n",K->x,trxs); 00159 /* *tracexs=trxs */ 00160 } else { 00161 K->x=0.0; 00162 } 00163 DSDPFunctionReturn(0); 00164 } 00165 00166 #undef __FUNCT__ 00167 #define __FUNCT__ "DSDPSetX" 00168 static int DSDPSetX( void *dspcone, double mu, DSDPVec y, DSDPVec dy){ 00169 RCone *K=(RCone*)dspcone; 00170 int info; 00171 double rr,dr,trxs,r; 00172 DSDPFunctionBegin; 00173 00174 info=DSDPVecGetR(y,&rr); DSDPCHKERR(info); 00175 info=DSDPVecGetR(dy,&dr); DSDPCHKERR(info); 00176 if (rr){ 00177 r=-1.0/rr; 00178 K->x=mu*(r-r*dr*r); 00179 trxs=K->x/r; 00180 DSDPLogInfo(0,2,"RESIDUAL X (Minimum Penalty Parameter): %4.4e, Trace(XS): %4.4e\n",K->x,trxs); 00181 /* *tracexs=trxs */ 00182 } else { 00183 K->x=0.0; 00184 } 00185 DSDPFunctionReturn(0); 00186 } 00187 00188 #undef __FUNCT__ 00189 #define __FUNCT__ "DSDPComputeRLog" 00190 static int DSDPComputeRLog(void *dspcone, double *logobj, double *logdet){ 00191 RCone *K=(RCone*)dspcone; 00192 DSDPFunctionBegin; 00193 *logdet=K->logr; 00194 *logobj=0; 00195 if (K->dualr<0){ 00196 *logdet=log(-K->dualr); 00197 K->logr=log(-K->dualr); 00198 } 00199 DSDPFunctionReturn(0); 00200 } 00201 00202 #undef __FUNCT__ 00203 #define __FUNCT__ "DSDPRANorm2" 00204 static int DSDPRANorm2(void *dspcone,DSDPVec Anorm2){ 00205 DSDPFunctionBegin; 00206 DSDPFunctionReturn(0); 00207 } 00208 00209 #undef __FUNCT__ 00210 #define __FUNCT__ "DSDPRMultiplyAdd" 00211 static int DSDPRMultiplyAdd(void *dspcone,double mu,DSDPVec vrow,DSDPVec vin,DSDPVec vout){ 00212 RCone *K=(RCone*)dspcone; 00213 int info; 00214 double v1,v2,rssr; 00215 DSDPFunctionBegin; 00216 if (K->dualr){ 00217 info=DSDPVecGetR(vrow,&v1);DSDPCHKERR(info); 00218 info=DSDPVecGetR(vin,&v2);DSDPCHKERR(info); 00219 rssr=v1*v2*mu/(K->dualr*K->dualr); 00220 info=DSDPVecAddR(vout,rssr);DSDPCHKERR(info); 00221 } 00222 DSDPFunctionReturn(0); 00223 } 00224 00225 #undef __FUNCT__ 00226 #define __FUNCT__ "DSDPRMonitor" 00227 static int DSDPRMonitor( void *dspcone, int tag){ 00228 DSDPFunctionBegin; 00229 DSDPFunctionReturn(0); 00230 } 00231 00232 static struct DSDPCone_Ops kops; 00233 static const char* matname="R Cone"; 00234 00235 #undef __FUNCT__ 00236 #define __FUNCT__ "RConeOperationsInitialize" 00237 static int RConeOperationsInitialize(struct DSDPCone_Ops* coneops){ 00238 int info; 00239 if (coneops==NULL) return 0; 00240 info=DSDPConeOpsInitialize(coneops); DSDPCHKERR(info); 00241 coneops->conehessian=DSDPRHessian; 00242 coneops->conesetup=DSDPSetupRCone; 00243 coneops->conesetup2=DSDPSetupRCone2; 00244 coneops->conedestroy=DSDPDestroyRCone; 00245 coneops->conecomputes=DSDPComputeRS; 00246 coneops->coneinverts=DSDPInvertRS; 00247 coneops->conesetxmaker=DSDPSetX; 00248 coneops->conecomputex=DSDPRX; 00249 coneops->conerhs=DSDPRHS; 00250 coneops->conemaxsteplength=DSDPComputeRStepLength; 00251 coneops->conelogpotential=DSDPComputeRLog; 00252 coneops->conesize=DSDPRSize; 00253 coneops->conesparsity=DSDPRSparsity; 00254 coneops->coneanorm2=DSDPRANorm2; 00255 coneops->conemonitor=DSDPRMonitor; 00256 coneops->conehmultiplyadd=DSDPRMultiplyAdd; 00257 coneops->id=19; 00258 coneops->name=matname; 00259 return 0; 00260 } 00261 00262 /* 00263 \fn int RConeSetType(RCone *rcone, DSDPPenalty UsePenalty); 00264 \brief Set penalty type. 00265 00266 \param dsdp the solver 00267 \param UsePenalty true or false 00268 */ 00269 #undef __FUNCT__ 00270 #define __FUNCT__ "RConeSetType" 00271 int RConeSetType(RCone *rcone, DSDPPenalty UsePenalty){ 00272 DSDPFunctionBegin; 00273 rcone->UsePenalty=UsePenalty; 00274 DSDPFunctionReturn(0); 00275 } 00276 /* 00277 \fn int RConeGetRX(RCone *rcone, double *rx); 00278 \brief Get slack of trace of matrix 00279 00280 \param rcone cone 00281 \param rx dual of r. 00282 Accurate only when r > 0. 00283 00284 */ 00285 #undef __FUNCT__ 00286 #define __FUNCT__ "RConeGetRX" 00287 int RConeGetRX(RCone *rcone, double *xtrace){ 00288 DSDPFunctionBegin; 00289 *xtrace=rcone->x; 00290 DSDPFunctionReturn(0); 00291 } 00292 00300 #undef __FUNCT__ 00301 #define __FUNCT__ "DSDPAddRCone" 00302 int DSDPAddRCone(DSDP dsdp, RCone **rrcone){ 00303 RCone *rcone; 00304 DSDPPenalty UsePenalty=DSDPInfeasible; 00305 int info; 00306 DSDPFunctionBegin; 00307 info=RConeOperationsInitialize(&kops); DSDPCHKERR(info); 00308 DSDPCALLOC1(&rcone,RCone,&info); DSDPCHKERR(info); 00309 info=RConeSetType(rcone,UsePenalty); DSDPCHKERR(info); 00310 rcone->dsdp=dsdp; 00311 rcone->logr=0; 00312 *rrcone=rcone; 00313 info=DSDPAddCone(dsdp,&kops,(void*)rcone); DSDPCHKERR(info); 00314 DSDPFunctionReturn(0); 00315 } 00316