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 } BCone;
00015
00016
00017 static int BComputeS(BCone *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 BCone *K=(BCone*)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
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 BCone *K=(BCone*)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
00072
00073
00074
00075
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(BCone *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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *K=(BCone*)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 BCone *rcone;
00290 int info;
00291 DSDPFunctionBegin;
00292 info=BConeOperationsInitialize(&kops); DSDPCHKERR(info);
00293 DSDPCALLOC1(&rcone,BCone,&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 #undef __FUNCT__
00304 #define __FUNCT__ "DSDPSetDualLowerBound"
00305 int DSDPSetDualLowerBound(DSDP dsdp, double dobjmin){
00306 int info;
00307 DSDPFunctionBegin;
00308 info = DSDPAddBCone(dsdp,dsdp->b,dobjmin);DSDPCHKERR(info);
00309 DSDPFunctionReturn(0);
00310 }