00001 #include "dsdpsdp.h" 00002 #include "dsdpsys.h" 00008 #undef __FUNCT__ 00009 #define __FUNCT__ "SDPConeComputeSS" 00010 00018 int SDPConeComputeSS(SDPCone sdpcone, int blockj, DSDPVec Y, DSDPVMat SS){ 00019 int info; 00020 DSDPFunctionBegin; 00021 info=DSDPVMatZeroEntries(SS); DSDPCHKBLOCKERR(blockj,info); 00022 info=DSDPBlockASum(&sdpcone->blk[blockj].ADATA,1,Y,SS); DSDPCHKBLOCKERR(blockj,info); 00023 DSDPFunctionReturn(0); 00024 } 00025 00026 #undef __FUNCT__ 00027 #define __FUNCT__ "SDPConeComputeS" 00028 00042 int SDPConeComputeS(SDPCone sdpcone, int blockj, double cc,double y[], int nvars, double r, int n, double s[], int nn){ 00043 int i,info; 00044 char UPLQ; 00045 DSDPVMat T; 00046 DSDPVec Y=sdpcone->Work; 00047 DSDPFunctionBegin; 00048 info=SDPConeCheckN(sdpcone,blockj,n);DSDPCHKBLOCKERR(blockj,info); 00049 info=SDPConeCheckM(sdpcone,nvars);DSDPCHKERR(info); 00050 if (n<1){DSDPFunctionReturn(0);} 00051 info=DSDPVecSetC(Y,-1.0*cc); 00052 info=DSDPVecSetR(Y,-r); 00053 for (i=0;i<nvars;i++){info=DSDPVecSetElement(Y,i+1,y[i]);} 00054 info=SDPConeGetStorageFormat(sdpcone,blockj,&UPLQ);DSDPCHKBLOCKERR(blockj,info); 00055 info=DSDPMakeVMatWithArray(UPLQ,s,nn,n,&T);DSDPCHKBLOCKERR(blockj,info); 00056 info=SDPConeComputeSS(sdpcone,blockj,Y,T);DSDPCHKBLOCKERR(blockj,info); 00057 info=DSDPVMatDestroy(&T);DSDPCHKBLOCKERR(blockj,info); 00058 DSDPFunctionReturn(0); 00059 } 00060 00061 #undef __FUNCT__ 00062 #define __FUNCT__ "SDPConeAddADotX" 00063 00075 int SDPConeAddADotX(SDPCone sdpcone, int blockj, double alpha,double x[], int nn, double adotx[], int m){ 00076 int info,n; 00077 char UPLQ; 00078 SDPblk *blk=sdpcone->blk; 00079 double scl=blk[blockj].ADATA.scl; 00080 DSDPVec ADOTX,YW2; 00081 DSDPVMat T; 00082 DSDPFunctionBegin; 00083 info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info); 00084 info=SDPConeCheckM(sdpcone,m-2);DSDPCHKERR(info); 00085 YW2=sdpcone->Work2; 00086 info=DSDPVecSet(alpha,YW2);DSDPCHKBLOCKERR(blockj,info); 00087 info=SDPConeGetBlockSize(sdpcone,blockj,&n);DSDPCHKBLOCKERR(blockj,info); 00088 if (n<1){DSDPFunctionReturn(0);} 00089 info=DSDPVecCreateWArray(&ADOTX,adotx,m); 00090 info=SDPConeGetStorageFormat(sdpcone,blockj,&UPLQ);DSDPCHKBLOCKERR(blockj,info); 00091 info=DSDPMakeVMatWithArray(UPLQ,x,nn,n,&T);DSDPCHKBLOCKERR(blockj,info); 00092 info=DSDPBlockADot(&blk[blockj].ADATA,1.0/scl,YW2,T,ADOTX);DSDPCHKBLOCKERR(blockj,info); 00093 info=DSDPVMatDestroy(&T);DSDPCHKBLOCKERR(blockj,info); 00094 DSDPFunctionReturn(0); 00095 } 00096 00097 #undef __FUNCT__ 00098 #define __FUNCT__ "SDPConeComputeXDot" 00099 00111 int SDPConeComputeXDot(SDPCone sdpcone, int blockj, DSDPVec Y,DSDPVMat X, DSDPVec AX, double* xtrace,double *xnorm, double *tracexs){ 00112 int info; 00113 SDPblk *blk=sdpcone->blk; 00114 DSDPVec YW2=sdpcone->Work2; 00115 double one=1.0,scl=blk[blockj].ADATA.scl; 00116 DSDPFunctionBegin; 00117 info=DSDPVecZero(YW2);DSDPCHKBLOCKERR(blockj,info); 00118 info=DSDPBlockADot(&blk[blockj].ADATA,-1.0/scl,Y,X,YW2);DSDPCHKBLOCKERR(blockj,info); 00119 info=DSDPVecGetR(YW2,xtrace);DSDPCHKBLOCKERR(blockj,info); 00120 info=DSDPVecSum(YW2,tracexs);DSDPCHKBLOCKERR(blockj,info); 00121 info=DSDPVMatNormF2(X,xnorm);DSDPCHKBLOCKERR(blockj,info); 00122 info=DSDPVecSet(one,YW2);DSDPCHKBLOCKERR(blockj,info); 00123 info=DSDPBlockADot(&blk[blockj].ADATA,1.0/scl,YW2,X,AX);DSDPCHKBLOCKERR(blockj,info); 00124 DSDPFunctionReturn(0); 00125 } 00126 00127 #undef __FUNCT__ 00128 #define __FUNCT__ "SDPConeComputeX3" 00129 00140 int SDPConeComputeX3(SDPCone sdpcone, int blockj, double mu, DSDPVec Y,DSDPVec DY,DSDPVMat X){ 00141 int info; 00142 double xshift=1e-12,xscale=1e-12; 00143 SDPblk *blk=sdpcone->blk; 00144 DSDPTruth psdefinite1=DSDP_FALSE,psdefinite2=DSDP_FALSE,full; 00145 DSDPDualMat SS; 00146 00147 DSDPFunctionBegin; 00148 SS=blk[blockj].SS; 00149 info=SDPConeComputeSS(sdpcone,blockj,Y,X);DSDPCHKBLOCKERR(blockj,info); 00150 info=DSDPDualMatSetArray(SS,X); DSDPCHKBLOCKERR(blockj,info); 00151 info=DSDPDualMatCholeskyFactor(SS,&psdefinite1); DSDPCHKBLOCKERR(blockj,info); 00152 if (psdefinite1 == DSDP_FALSE){ 00153 DSDPLogInfo(0,2,"Primal SDP Block %2.0f not PSD\n",blockj); 00154 } 00155 info=DSDPDualMatInvert(SS);DSDPCHKBLOCKERR(blockj,info); 00156 info=SDPConeComputeXX(sdpcone,blockj,DY,mu,SS,X);DSDPCHKBLOCKERR(blockj,info); 00157 info=DSDPDualMatIsFull(SS,&full);DSDPCHKBLOCKERR(blockj,info); 00158 psdefinite2=DSDP_FALSE; 00159 while (full==DSDP_TRUE && psdefinite2==DSDP_FALSE && xscale<2e-1){ 00160 info=DSDPVMatShiftDiagonal(X,xshift); DSDPCHKBLOCKERR(blockj,info); 00161 info=DSDPVMatScaleDiagonal(X,1.0+xscale); DSDPCHKBLOCKERR(blockj,info); 00162 DSDPLogInfo(0,10,"VMat: shift diagonal: %4.2e, scale diagonal: %4.2e.\n",xshift,1+xscale); 00163 info=DSDPDualMatSetArray(SS,X); DSDPCHKBLOCKERR(blockj,info); 00164 info=DSDPDualMatCholeskyFactor(SS,&psdefinite2); DSDPCHKBLOCKERR(blockj,info); 00165 xshift*=10;xscale*=10; 00166 } 00167 if (full==DSDP_FALSE){ 00168 xshift=1e-12,xscale=1e-10; 00169 info=DSDPVMatShiftDiagonal(X,xshift); DSDPCHKBLOCKERR(blockj,info); 00170 info=DSDPVMatScaleDiagonal(X,1.0+xscale); DSDPCHKBLOCKERR(blockj,info); 00171 DSDPLogInfo(0,10,"XMat: shift diagonal: %4.2e, scale diagonal: %4.2e.\n",xshift,1+xscale); 00172 } 00173 DSDPFunctionReturn(0); 00174 } 00175 00176 00177 #undef __FUNCT__ 00178 #define __FUNCT__ "SDPConeComputeX" 00179 00191 int SDPConeComputeX(SDPCone sdpcone, int blockj, int n, double x[], int nn){ 00192 int info; 00193 double mu=sdpcone->xmakermu; 00194 double xnorm,xtrace,trxs; 00195 char UPLQ; 00196 DSDPVec DY=sdpcone->DYX,Y=sdpcone->YX,AX=sdpcone->Work; 00197 DSDPVMat T; 00198 00199 DSDPFunctionBegin; 00200 info=SDPConeCheckN(sdpcone,blockj,n);DSDPCHKBLOCKERR(blockj,info); 00201 if (n<1){DSDPFunctionReturn(0);} 00202 info=SDPConeGetStorageFormat(sdpcone,blockj,&UPLQ);DSDPCHKBLOCKERR(blockj,info); 00203 info=DSDPMakeVMatWithArray(UPLQ,x,nn,n,&T);DSDPCHKBLOCKERR(blockj,info); 00204 info=SDPConeComputeX3(sdpcone,blockj,mu,Y,DY,T);DSDPCHKBLOCKERR(blockj,info); 00205 info=SDPConeComputeXDot(sdpcone,blockj,Y,T,AX,&xtrace,&xnorm,&trxs);DSDPCHKBLOCKERR(blockj,info); 00206 info=DSDPVMatDestroy(&T);DSDPCHKBLOCKERR(blockj,info); 00207 DSDPFunctionReturn(0); 00208 } 00209 00210 #undef __FUNCT__ 00211 #define __FUNCT__ "SDPConeViewX" 00212 00223 int SDPConeViewX(SDPCone sdpcone, int blockj, int n, double x[], int nn){ 00224 int info; 00225 char UPLQ; 00226 DSDPVMat T; 00227 00228 DSDPFunctionBegin; 00229 info=SDPConeCheckN(sdpcone,blockj,n);DSDPCHKBLOCKERR(blockj,info); 00230 info=SDPConeGetStorageFormat(sdpcone,blockj,&UPLQ);DSDPCHKBLOCKERR(blockj,info); 00231 info=DSDPMakeVMatWithArray(UPLQ,x,nn,n,&T);DSDPCHKBLOCKERR(blockj,info); 00232 info=DSDPVMatView(T);DSDPCHKBLOCKERR(blockj,info); 00233 info=DSDPVMatDestroy(&T);DSDPCHKBLOCKERR(blockj,info); 00234 DSDPFunctionReturn(0); 00235 } 00236 00237 #undef __FUNCT__ 00238 #define __FUNCT__ "SDPConeXVMultiply" 00239 00251 int SDPConeXVMultiply(SDPCone sdpcone, int blockj, double vin[], double vout[], int n){ 00252 int info; 00253 SDPblk *blk=sdpcone->blk; 00254 SDPConeVec V1,V2,V3,V4; 00255 DSDPDualMat S,SS; 00256 00257 DSDPFunctionBegin; 00258 info=SDPConeCheckN(sdpcone,blockj,n);DSDPCHKBLOCKERR(blockj,info); 00259 if (sdpcone->blk[blockj].n>1){ 00260 S=blk[blockj].S; SS=blk[blockj].SS; 00261 V2=blk[blockj].W; V3=blk[blockj].W2; 00262 info=SDPConeVecCreateWArray(&V1,vin,n); 00263 info=SDPConeVecCreateWArray(&V4,vout,n); 00264 if (0){ 00265 info=DSDPDualMatCholeskySolveForward(S,V1,V3);DSDPCHKERR(info); 00266 info=SDPConeVecScale(sqrt(sdpcone->xmakermu),V3);DSDPCHKERR(info); 00267 info=DSDPDualMatCholeskySolveBackward(S,V3,V2);DSDPCHKERR(info); 00268 info=DSDPDualMatCholeskyBackwardMultiply(SS,V2,V1);DSDPCHKERR(info); 00269 } 00270 info=DSDPDualMatCholeskyForwardMultiply(SS,V1,V2);DSDPCHKERR(info); 00271 info=DSDPDualMatCholeskySolveForward(S,V2,V3);DSDPCHKERR(info); 00272 info=SDPConeVecScale(sqrt(sdpcone->xmakermu),V3);DSDPCHKERR(info); 00273 info=DSDPDualMatCholeskySolveBackward(S,V3,V4);DSDPCHKERR(info); 00274 } 00275 DSDPFunctionReturn(0); 00276 } 00277 00278 #undef __FUNCT__ 00279 #define __FUNCT__ "SDPConeAddXVAV" 00280 00292 int SDPConeAddXVAV(SDPCone sdpcone, int blockj, double vin[], int n, double sum[], int mm){ 00293 int info; 00294 SDPblk *blk=sdpcone->blk; 00295 SDPConeVec V2; 00296 DSDPVec W,Wout; 00297 DSDPFunctionBegin; 00298 info=SDPConeCheckN(sdpcone,blockj,n);DSDPCHKBLOCKERR(blockj,info); 00299 W=sdpcone->Work; 00300 info=DSDPVecSet(1.0,sdpcone->Work);DSDPCHKBLOCKERR(blockj,info); 00301 if (sdpcone->blk[blockj].n>1){ 00302 info=SDPConeVecCreateWArray(&V2,vin,n);DSDPCHKERR(info); 00303 info=DSDPVecCreateWArray(&Wout,sum,mm);DSDPCHKERR(info); 00304 info=DSDPBlockvAv(&blk[blockj].ADATA,1.0,sdpcone->Work,V2,Wout);DSDPCHKBLOCKERR(blockj,info); 00305 } 00306 DSDPFunctionReturn(0); 00307 } 00308 00309 00310 00311 #undef __FUNCT__ 00312 #define __FUNCT__ "SDPConeComputeXV" 00313 00325 int SDPConeComputeXV(SDPCone sdpcone, int blockj, int *derror){ 00326 00327 int info; double rr; 00328 DSDPVec Y,DY,W; 00329 SDPblk *blk=sdpcone->blk; 00330 DSDPTruth psdefinite1=DSDP_FALSE,psdefinite2=DSDP_FALSE; 00331 DSDPDualMat S,SS; 00332 DSDPVMat T; 00333 00334 DSDPFunctionBegin; 00335 *derror=0; 00336 info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKBLOCKERR(blockj,info); 00337 if (sdpcone->blk[blockj].n>1){ 00338 Y=sdpcone->YX; DY=sdpcone->DYX; W=sdpcone->Work; 00339 T=blk[blockj].T; S=blk[blockj].S; SS=blk[blockj].SS; 00340 info=DSDPVecWAXPY(W,-1.0,DY,Y);DSDPCHKBLOCKERR(blockj,info); 00341 00342 while (psdefinite1==DSDP_FALSE){ 00343 info=DSDPVecGetR(W,&rr); 00344 info=DSDPVecSetR(W,10*rr-1e-12); 00345 info=SDPConeComputeSS(sdpcone,blockj,W,T);DSDPCHKBLOCKERR(blockj,info); 00346 info=DSDPDualMatSetArray(SS,T); DSDPCHKBLOCKERR(blockj,info); 00347 info=DSDPDualMatCholeskyFactor(SS,&psdefinite1); DSDPCHKBLOCKERR(blockj,info); 00348 } 00349 00350 while (psdefinite2==DSDP_FALSE){ 00351 info=SDPConeComputeSS(sdpcone,blockj,Y,T);DSDPCHKBLOCKERR(blockj,info); 00352 info=DSDPDualMatSetArray(S,T); DSDPCHKBLOCKERR(blockj,info); 00353 info=DSDPDualMatCholeskyFactor(S,&psdefinite2); DSDPCHKBLOCKERR(blockj,info); 00354 if (psdefinite2==DSDP_FALSE){ 00355 info=DSDPVecGetR(Y,&rr); 00356 info=DSDPVecSetR(Y,10*rr-1e-15); 00357 } 00358 } 00359 if (psdefinite1==DSDP_FALSE || psdefinite2==DSDP_FALSE) *derror=1; 00360 } 00361 DSDPFunctionReturn(0); 00362 }