sdpcone.c

Go to the documentation of this file.
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 }

Generated on Wed Nov 11 20:41:02 2009 for DSDP by  doxygen 1.6.1