DSDP
|
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