DSDP
|
00001 #include "dsdpcone_impl.h" 00002 #include "dsdpsys.h" 00008 typedef struct { 00009 DSDPVec b,bb,T; 00010 double dmin; 00011 double pss,dss; 00012 DSDP dsdp; 00013 DSDPTruth useit; 00014 } BDCone; 00015 00016 00017 static int BComputeS(BDCone *K, DSDPVec v, double *ss){ 00018 int info; 00019 DSDPFunctionBegin; 00020 info=DSDPVecDot(K->bb,v,ss);DSDPCHKERR(info); 00021 DSDPFunctionReturn(0); 00022 } 00023 00024 #undef __FUNCT__ 00025 #define __FUNCT__ "DSDPRHessian" 00026 static int DSDPRHessian( void *dspcone, double mu, DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2){ 00027 BDCone *K=(BDCone*)dspcone; 00028 double bb,dd,ss=K->dss; 00029 int info,i,m,ncols; 00030 DSDPVec T=K->T,b=K->bb; 00031 DSDPFunctionBegin; 00032 if (K->useit){ 00033 info=DSDPVecGetSize(b,&m);DSDPCHKERR(info); 00034 for (i=0;i<m;i++){ 00035 info=DSDPVecGetElement(b,i,&bb);DSDPCHKERR(info); 00036 if (bb==0) continue; 00037 00038 info=DSDPSchurMatRowColumnScaling(M,i,T,&ncols); DSDPCHKERR(info); 00039 if (ncols==0) continue; 00040 00041 info=DSDPVecGetElement(T,i,&dd);DSDPCHKERR(info); 00042 info=DSDPVecAddElement(vrhs2,i,-bb*dd*mu/ss);DSDPCHKERR(info); 00043 00044 info=DSDPVecPointwiseMult(T,b,T);DSDPCHKERR(info); 00045 00046 info=DSDPVecScale(mu*bb/(ss*ss),T);DSDPCHKERR(info); 00047 info=DSDPSchurMatAddRow(M,i,1.0,T);DSDPCHKERR(info); 00048 /* 00049 info=DSDPSchurMatAddRow(M,i,mu*bb/(ss*ss),T);DSDPCHKERR(info); 00050 */ 00051 if (i==-m-1){info=DSDPVecView(T);} 00052 } 00053 } 00054 DSDPFunctionReturn(0); 00055 } 00056 00057 static int DSDPRRHS( void *dcone, double mu, DSDPVec vrow, DSDPVec vrhs1, DSDPVec vrhs2){ 00058 BDCone *K=(BDCone*)dcone; 00059 double bb,dd,ss=K->dss; 00060 int info,i,m; 00061 DSDPVec b=K->bb; 00062 DSDPFunctionBegin; 00063 if (K->useit){ 00064 info=DSDPVecGetSize(b,&m);DSDPCHKERR(info); 00065 for (i=0;i<m;i++){ 00066 info=DSDPVecGetElement(b,i,&bb);DSDPCHKERR(info); 00067 info=DSDPVecGetElement(vrow,i,&dd);DSDPCHKERR(info); 00068 info=DSDPVecAddElement(vrhs2,i,-bb*dd*mu/ss);DSDPCHKERR(info); 00069 } 00070 /* 00071 info=DSDPVecGetR(b,&bb);DSDPCHKERR(info); 00072 info=DSDPVecGetR(vrhs3,&dd);DSDPCHKERR(info); 00073 info=DSDPVecPointwiseMult(vrow,b,T);DSDPCHKERR(info); 00074 info=DSDPVecScale(mu*bb/(ss*ss),T);DSDPCHKERR(info); 00075 info=DSDPVecAXPY(dd,T,vrhs3);DSDPCHKERR(info); 00076 */ 00077 } 00078 DSDPFunctionReturn(0); 00079 } 00080 00081 #undef __FUNCT__ 00082 #define __FUNCT__ "DSDPSetupBCone" 00083 static int DSDPSetupBCone(void* dspcone,DSDPVec y){ 00084 DSDPFunctionBegin; 00085 DSDPFunctionReturn(0); 00086 } 00087 00088 00089 #undef __FUNCT__ 00090 #define __FUNCT__ "DSDPSetDRData" 00091 static int DSDPSetDRData(BDCone *K){ 00092 int info; 00093 DSDPFunctionBegin; 00094 info=DSDPVecCopy(K->b,K->bb);DSDPCHKERR(info); 00095 info=DSDPVecSetC(K->bb,K->dmin);DSDPCHKERR(info); 00096 info=DSDPVecSetR(K->bb,-1.0);DSDPCHKERR(info); 00097 DSDPFunctionReturn(0); 00098 } 00099 00100 #undef __FUNCT__ 00101 #define __FUNCT__ "DSDPSetupBCone2" 00102 static int DSDPSetupBCone2(void* dspcone, DSDPVec y, DSDPSchurMat M){ 00103 BDCone *K=(BDCone*)dspcone; 00104 int info; 00105 DSDPFunctionBegin; 00106 info=DSDPVecDuplicate(K->b,&K->T);DSDPCHKERR(info); 00107 info=DSDPVecDuplicate(K->b,&K->bb);DSDPCHKERR(info); 00108 info=DSDPSetDRData(K);DSDPCHKERR(info); 00109 DSDPFunctionReturn(0); 00110 } 00111 00112 00113 #undef __FUNCT__ 00114 #define __FUNCT__ "DSDPDestroyBCone" 00115 static int DSDPDestroyBCone(void* dspcone){ 00116 BDCone *K=(BDCone*)dspcone; 00117 int info; 00118 DSDPFunctionBegin; 00119 info=DSDPVecDestroy(&K->T);DSDPCHKERR(info); 00120 info=DSDPVecDestroy(&K->bb);DSDPCHKERR(info); 00121 DSDPFREE(&dspcone,&info);DSDPCHKERR(info); 00122 DSDPFunctionReturn(0); 00123 } 00124 00125 00126 #undef __FUNCT__ 00127 #define __FUNCT__ "DSDPRSize" 00128 static int DSDPRSize(void*dspcone,double*n){ 00129 DSDPFunctionBegin; 00130 *n=1.0; 00131 DSDPFunctionReturn(0); 00132 } 00133 00134 #undef __FUNCT__ 00135 #define __FUNCT__ "DSDPRSparsity" 00136 static int DSDPRSparsity(void*dspcone,int row, int *tnnz, int rnnz[], int m){ 00137 BDCone *K=(BDCone*)dspcone; 00138 int i,info;double dd; 00139 DSDPFunctionBegin; 00140 *tnnz=0; 00141 00142 info=DSDPVecGetElement(K->b,row,&dd);DSDPCHKERR(info); 00143 if (dd){ 00144 for (i=0;i<m;i++){ 00145 info=DSDPVecGetElement(K->b,i,&dd);DSDPCHKERR(info); 00146 if (dd!=0){rnnz[i]++; (*tnnz)++;} 00147 } 00148 } 00149 DSDPFunctionReturn(0); 00150 } 00151 00152 #undef __FUNCT__ 00153 #define __FUNCT__ "DSDPComputeRS" 00154 static int DSDPComputeRS(void *dspcone, DSDPVec Y, DSDPDualFactorMatrix flag, DSDPTruth *ispsdefinite){ 00155 BDCone *K=(BDCone*)dspcone; 00156 int info; 00157 double ss; 00158 DSDPFunctionBegin; 00159 info=BComputeS(K,Y,&ss);DSDPCHKERR(info); 00160 if (ss>0){ *ispsdefinite=DSDP_TRUE; } else { *ispsdefinite=DSDP_FALSE;} 00161 if (flag==DUAL_FACTOR){ K->dss=ss; } else { K->pss=ss;} 00162 DSDPLogInfo(0,2,"DOBJCone SS: %4.4e \n",ss); 00163 DSDPFunctionReturn(0); 00164 } 00165 00166 #undef __FUNCT__ 00167 #define __FUNCT__ "DSDPInvertRS" 00168 static int DSDPInvertRS(void *dspcone){ 00169 DSDPFunctionBegin; 00170 DSDPFunctionReturn(0); 00171 } 00172 00173 00174 #undef __FUNCT__ 00175 #define __FUNCT__ "DSDPComputeRStepLength" 00176 static int DSDPComputeRStepLength(void *dspcone, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){ 00177 BDCone *K=(BDCone*)dspcone; 00178 double ds,ss,rt=1.0e30; 00179 int info; 00180 00181 DSDPFunctionBegin; 00182 info=BComputeS(K,DY,&ds);DSDPCHKERR(info); 00183 if (flag==DUAL_FACTOR){ ss=K->dss; } else { ss=K->pss;} 00184 if (ds<0) rt=-ss/ds; 00185 if (K->useit){ 00186 *maxsteplength=rt; 00187 } 00188 DSDPFunctionReturn(0); 00189 } 00190 00191 #undef __FUNCT__ 00192 #define __FUNCT__ "DSDPSetX" 00193 static int DSDPSetX( void *dspcone, double mu, DSDPVec y, DSDPVec dy){ 00194 DSDPFunctionBegin; 00195 DSDPFunctionReturn(0); 00196 } 00197 #undef __FUNCT__ 00198 #define __FUNCT__ "DSDPRX" 00199 static int DSDPRX( void *dspcone, double mu, DSDPVec y, DSDPVec dy, DSDPVec AX,double*tracxs){ 00200 BDCone *K=(BDCone*)dspcone; 00201 double x,dss,ss=K->dss; 00202 int info; 00203 DSDPFunctionBegin; 00204 00205 info=BComputeS(K,y,&ss);DSDPCHKERR(info); 00206 ss=1.0/ss; 00207 info=BComputeS(K,dy,&dss);DSDPCHKERR(info); 00208 x=mu*(ss+ss*dss*ss); 00209 DSDPLogInfo(0,2,"DOBJCone SS: %4.4e, RESIDUAL X: %4.4e\n",1.0/ss,x); 00210 if (fabs(x*ss)>1.0 && mu < 1) printf("Check Dual Min Bound\n"); 00211 info=DSDPVecAXPY(-x,K->bb,AX);DSDPCHKERR(info); 00212 DSDPFunctionReturn(0); 00213 } 00214 00215 #undef __FUNCT__ 00216 #define __FUNCT__ "DSDPComputeRLog" 00217 static int DSDPComputeRLog(void *dspcone, double *logobj, double *logdet){ 00218 BDCone *K=(BDCone*)dspcone; 00219 DSDPFunctionBegin; 00220 *logobj=0; 00221 *logdet=log(K->dss); 00222 DSDPFunctionReturn(0); 00223 } 00224 00225 #undef __FUNCT__ 00226 #define __FUNCT__ "DSDPRANorm2" 00227 static int DSDPRANorm2(void *dspcone, DSDPVec Anorm2){ 00228 DSDPFunctionBegin; 00229 DSDPFunctionReturn(0); 00230 } 00231 00232 00233 #undef __FUNCT__ 00234 #define __FUNCT__ "DSDPRMultiplyAdd" 00235 static int DSDPRMultiplyAdd(void *dspcone,double mu,DSDPVec vrow,DSDPVec vin,DSDPVec vout){ 00236 BDCone *K=(BDCone*)dspcone; 00237 DSDPVec T=K->T; 00238 int info; 00239 double dd,ss=K->dss; 00240 DSDPFunctionBegin; 00241 info=DSDPVecDot(vin,K->bb,&dd);DSDPCHKERR(info); 00242 dd=-mu*dd/(ss*ss); 00243 info=DSDPVecPointwiseMult(K->bb,vrow,T);DSDPCHKERR(info); 00244 info=DSDPVecAXPY(dd,T,vout);DSDPCHKERR(info); 00245 DSDPFunctionReturn(0); 00246 } 00247 00248 00249 #undef __FUNCT__ 00250 #define __FUNCT__ "DSDPRMonitor" 00251 static int DSDPRMonitor( void *dspcone, int tag){ 00252 DSDPFunctionBegin; 00253 DSDPFunctionReturn(0); 00254 } 00255 00256 static struct DSDPCone_Ops kops; 00257 static const char* matname="Dual Obj Cone"; 00258 00259 #undef __FUNCT__ 00260 #define __FUNCT__ "BConeOperationsInitialize" 00261 static int BConeOperationsInitialize(struct DSDPCone_Ops* coneops){ 00262 int info; 00263 if (coneops==NULL) return 0; 00264 info=DSDPConeOpsInitialize(coneops); DSDPCHKERR(info); 00265 coneops->conehessian=DSDPRHessian; 00266 coneops->conesetup=DSDPSetupBCone; 00267 coneops->conesetup2=DSDPSetupBCone2; 00268 coneops->conedestroy=DSDPDestroyBCone; 00269 coneops->conecomputes=DSDPComputeRS; 00270 coneops->coneinverts=DSDPInvertRS; 00271 coneops->conecomputex=DSDPRX; 00272 coneops->conesetxmaker=DSDPSetX; 00273 coneops->conemaxsteplength=DSDPComputeRStepLength; 00274 coneops->conelogpotential=DSDPComputeRLog; 00275 coneops->conesize=DSDPRSize; 00276 coneops->conesparsity=DSDPRSparsity; 00277 coneops->coneanorm2=DSDPRANorm2; 00278 coneops->conemonitor=DSDPRMonitor; 00279 coneops->conehmultiplyadd=DSDPRMultiplyAdd; 00280 coneops->conerhs=DSDPRRHS; 00281 coneops->id=119; 00282 coneops->name=matname; 00283 return 0; 00284 } 00285 00286 #undef __FUNCT__ 00287 #define __FUNCT__ "DSDPAddBCone" 00288 int DSDPAddBCone(DSDP dsdp, DSDPVec bb, double dmin){ 00289 BDCone *rcone; 00290 int info; 00291 DSDPFunctionBegin; 00292 info=BConeOperationsInitialize(&kops); DSDPCHKERR(info); 00293 DSDPCALLOC1(&rcone,BDCone,&info); DSDPCHKERR(info); 00294 rcone->b=bb; 00295 rcone->dmin=dmin; 00296 rcone->dsdp=dsdp; 00297 rcone->useit=DSDP_TRUE; 00298 info=DSDPAddCone(dsdp,&kops,(void*)rcone); DSDPCHKERR(info); 00299 DSDPFunctionReturn(0); 00300 } 00301 00302 #include "dsdp.h" 00303 #include "dsdp5.h" 00304 00305 #undef __FUNCT__ 00306 #define __FUNCT__ "DSDPSetDualLowerBound" 00307 int DSDPSetDualLowerBound(DSDP dsdp, double dobjmin){ 00308 int info; 00309 DSDPFunctionBegin; 00310 info = DSDPAddBCone(dsdp,dsdp->b,dobjmin);DSDPCHKERR(info); 00311 DSDPFunctionReturn(0); 00312 }