sdpcompute.c

Go to the documentation of this file.
00001 #include "dsdpsdp.h"
00002 #include "dsdpsys.h"
00008 /*
00009 static int tt1=0,tt2=0;
00010       DSDPEventLogBegin(tt1);
00011       DSDPEventLogBegin(tt2);
00012       DSDPEventLogEnd(tt2);
00013       DSDPEventLogEnd(tt1);
00014   if (tt1==0){DSDPEventLogRegister("SDP T1 Event",&tt1);}
00015   if (tt2==0){DSDPEventLogRegister("SDP T2 Event",&tt2);}
00016 */
00017 #undef __FUNCT__  
00018 #define __FUNCT__ "SDPConeComputeHessian"
00019 
00030 int SDPConeComputeHessian( SDPCone sdpcone, double mu, DSDPSchurMat M,  DSDPVec vrhs1, DSDPVec vrhs2){
00031 
00032   int i,k,kt,kk,m,n,rank,info;
00033   int ncols,ii,idA;
00034   double rtemp,ack,ggamma,bmu,scl;
00035   double rhs1i,rhs2i;
00036   DSDPTruth method1;
00037   SDPConeVec W,W2;
00038   DSDPVec MRowI=sdpcone->Work, Select=sdpcone->Work2;
00039   DSDPDataMat AA;
00040   DSDPDualMat S;
00041   DSDPVMat T;
00042   DSDPDataTranspose ATranspose=sdpcone->ATR;
00043   SDPblk *blk=sdpcone->blk;
00044   DSDPIndex IS;
00045 
00046   /* Evaluate M */
00047   DSDPFunctionBegin;
00048   SDPConeValid(sdpcone);
00049   info=DSDPVecGetSize(vrhs1,&m);DSDPCHKERR(info);
00050 
00051   for (i=0; i<m; i++){  /* One row at a time */
00052 
00053     /* Which Coluns */
00054     rhs1i=0;rhs2i=0;
00055     info=DSDPVecZero(MRowI);DSDPCHKERR(info);
00056     info=DSDPSchurMatRowColumnScaling(M,i,Select,&ncols); DSDPCHKERR(info); 
00057     if (ncols==0){continue;}
00058  
00059     for (kt=0; kt<ATranspose.nnzblocks[i]; kt++){ /* Loop over all blocks */
00060       kk=ATranspose.nzblocks[i][kt];
00061       idA=ATranspose.idA[i][kt];
00062       info=DSDPBlockGetMatrix(&blk[kk].ADATA,idA,&ii,&scl,&AA);DSDPCHKBLOCKERR(kk,info);
00063       if (ii!=i){DSDPSETERR2(8,"Data Transpose Error: var %d does not equal %d.\n",i,ii);}
00064       info = DSDPDataMatGetRank(AA,&rank,blk[kk].n);DSDPCHKBLOCKERR(kk,info);
00065       if (rank==0) continue;
00066 
00067       T=blk[kk].T; S=blk[kk].S; W=blk[kk].W; W2=blk[kk].W2;
00068       n=blk[kk].n; ggamma=blk[kk].gammamu; bmu=blk[kk].bmu; IS=blk[kk].IS;
00069 
00070       method1=DSDP_TRUE;  /* Simple heuristic */
00071       if (rank==1) method1=DSDP_FALSE;
00072       if (rank==2 && ncols<=n) method1=DSDP_FALSE;
00073       if (rank*rank*ncols<=n+1)method1=DSDP_FALSE;
00074       if (ncols*blk[kk].nnz < n*n/10) method1=DSDP_FALSE;
00075       if (ncols==1 && i==m-1)method1=DSDP_FALSE;
00076       if (n<5) method1=DSDP_TRUE;
00077       if (0==1) method1=DSDP_FALSE;
00078       if (method1==DSDP_TRUE){info=DSDPVMatZeroEntries(T);DSDPCHKBLOCKERR(kk,info);}
00079       for (k=0; k<rank; k++){
00080         
00081         info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKBLOCKERR(kk,info);
00082         if (ack==0.0) continue;
00083         ack*=scl;
00084         info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKBLOCKERR(kk,info);
00085 
00086         /* RHS terms */
00087         info = SDPConeVecDot(W,W2,&rtemp); DSDPCHKBLOCKERR(kk,info);
00088         if (rtemp==0.0) continue;
00089         rhs1i+=rtemp*ack*bmu; rhs2i+=rtemp*ack*ggamma*mu;
00090         ack*=(ggamma+bmu);
00091 
00092         if (method1==DSDP_TRUE){
00093           info=DSDPVMatAddOuterProduct(T,ack*mu,W2);DSDPCHKBLOCKERR(kk,info);
00094         } else {
00095           info=DSDPBlockvAv(&blk[kk].ADATA,ack*mu,Select,W2,MRowI);DSDPCHKBLOCKERR(kk,info);
00096         } /* End row computations for rank kk of block kk */
00097  
00098       }   /* End row computations for all of block kk     */
00099 
00100       if (method1==DSDP_TRUE){
00101         info=DSDPBlockADot(&blk[kk].ADATA,1.0,Select,T,MRowI);DSDPCHKBLOCKERR(kk,info);
00102       }   /* End row computations for all of block ll     */
00103     }     /* End row computations for all blocks          */
00104     info=DSDPVecAddElement(vrhs1,i,rhs1i);DSDPCHKERR(info);
00105     info=DSDPVecAddElement(vrhs2,i,rhs2i);DSDPCHKERR(info);
00106     info=DSDPSchurMatAddRow(M,i,1.0,MRowI);DSDPCHKERR(info);
00107   }
00108 
00109   DSDPFunctionReturn(0);
00110 }
00111 
00112 #undef __FUNCT__  
00113 #define __FUNCT__ "SDPConeComputeRHS"
00114 
00125 int SDPConeComputeRHS( SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vrhs1, DSDPVec vrhs2){
00126 
00127   int info,i,ii,k,rank,nnzmats;
00128   double dtmp,dyiscale=1,ack,scl,rtemp;
00129   SDPblk *sdp=&sdpcone->blk[blockj];
00130   SDPConeVec W=sdp->W,W2=sdp->W2;
00131   DSDPDataMat AA;
00132   DSDPVMat T=sdp->T;
00133   DSDPDualMat S=sdp->S;
00134   DSDPTruth method1;
00135   DSDPIndex IS=sdp->IS;
00136 
00137   DSDPFunctionBegin;
00138 
00139   info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
00140   method1=DSDP_TRUE;
00141   method1=DSDP_FALSE;
00142 
00143   if (method1==DSDP_TRUE){
00144     info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
00145     for (i=0;i<nnzmats;i++){
00146       info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
00147       info=DSDPVecGetElement(vrow,ii,&dyiscale);DSDPCHKVARERR(ii,info);
00148       if (dyiscale==0) continue;
00149       info=DSDPDataMatGetRank(AA,&rank,sdp->n);DSDPCHKVARERR(ii,info);
00150       for (k=0; k<rank; k++){
00151         info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKVARERR(ii,info);
00152         if (ack==0) continue;
00153         info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
00154         info=SDPConeVecDot(W,W2,&rtemp); DSDPCHKVARERR(ii,info);
00155         dtmp=rtemp*ack*mu*dyiscale*scl;
00156         info=DSDPVecAddElement(vrhs2,ii,dtmp);DSDPCHKVARERR(ii,info);
00157       }
00158     }
00159     
00160   } else {
00161     info=DSDPVMatZeroEntries(T);DSDPCHKERR(info);
00162     info=DSDPDualMatInverseAdd(S,mu,T);DSDPCHKERR(info);
00163     info=DSDPBlockADot(&sdp->ADATA,1.0,vrow,T,vrhs2);DSDPCHKERR(info);
00164   }
00165   
00166   DSDPFunctionReturn(0);
00167 }
00168 
00169 #undef __FUNCT__  
00170 #define __FUNCT__ "SDPConeMultiply"
00171 
00182 int SDPConeMultiply(SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout){
00183 
00184   int info,i,ii,k,rank,nnzmats;
00185   double dd2,dtmp,dyiscale,ack,scl,vv;
00186   SDPblk *sdp=&sdpcone->blk[blockj];
00187   SDPConeVec W=sdp->W,W2=sdp->W2;
00188   DSDPDataMat AA;
00189   DSDPVMat T=sdp->T;
00190   DSDPDSMat DS=sdp->DS;
00191   DSDPIndex IS=sdp->IS;
00192   DSDPDualMat S=sdp->S;
00193 
00194   DSDPFunctionBegin;
00195   info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
00196   info=DSDPVMatZeroEntries(T); DSDPCHKERR(info);
00197   info=DSDPBlockASum(&sdp->ADATA,-1.0,vin,T); DSDPCHKERR(info);
00198   info=DSDPDSMatSetArray(DS,T); DSDPCHKERR(info);
00199   info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
00200   for (i=0;i<nnzmats;i++){
00201     info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
00202     info=DSDPVecGetElement(vrow,ii,&dyiscale);DSDPCHKVARERR(ii,info);
00203     if (dyiscale==0) continue;
00204 
00205     info=DSDPDataMatGetRank(AA,&rank,sdp->n);DSDPCHKVARERR(ii,info);
00206     for (dd2=0,k=0; k<rank; k++){
00207       info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKVARERR(ii,info);
00208       if (ack==0) continue;
00209       info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
00210       info=DSDPDSMatVecVec(DS,W2,&vv);DSDPCHKVARERR(ii,info);
00211       dd2+=vv*ack;
00212     }
00213     dtmp = dd2 * dyiscale *mu *scl;
00214     info=DSDPVecAddElement(vout,ii,dtmp);DSDPCHKVARERR(ii,info);
00215   }
00216 
00217   DSDPFunctionReturn(0);
00218 }
00219 
00220 
00221 
00222 #undef __FUNCT__  
00223 #define __FUNCT__ "SDPConeComputeXX"
00224 
00235 int SDPConeComputeXX(SDPCone sdpcone, int blockj, DSDPVec DY, double mu,DSDPDualMat S, DSDPVMat X){
00236   
00237   int info, i, ii,k, rank, n, nnzmats;
00238   double dtmp,dyiscale,ack,scl;
00239   SDPblk *sdp=&sdpcone->blk[blockj];
00240   SDPConeVec W=sdp->W,W2=sdp->W2;
00241   DSDPDataMat AA;
00242   DSDPIndex IS=sdp->IS;
00243 
00244   DSDPFunctionBegin;
00245   info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
00246   mu=mu*sdp->gammamu; n=sdp->n;
00247   info=DSDPVMatZeroEntries(X);DSDPCHKERR(info);
00248   info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
00249   for (i=0;i<nnzmats;i++){
00250     info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKVARERR(ii,info);
00251     info=DSDPVecGetElement(DY,ii,&dyiscale);DSDPCHKVARERR(ii,info);
00252     if (dyiscale==0) continue;
00253     info=DSDPDataMatGetRank(AA,&rank,n);DSDPCHKVARERR(ii,info);
00254     for (k=0; k<rank; k++){
00255       info = DSDPDataMatGetEig(AA,k,W,IS,&ack);DSDPCHKVARERR(ii,info);
00256       if (ack==0) continue;
00257       info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
00258       dtmp = ack * dyiscale * mu * scl;
00259       info=DSDPVMatAddOuterProduct(X,dtmp,W2);DSDPCHKVARERR(ii,info);
00260     }
00261   }
00262 
00263   info=DSDPDualMatInverseAdd(S,mu,X);DSDPCHKERR(info);
00264   DSDPFunctionReturn(0);
00265 }
00266 

Generated on Thu May 22 09:42:47 2008 for DSDP by  doxygen 1.5.5