Logo Search packages:      
Sourcecode: dsdp version File versions  Download package

int SDPConeComputeHessian ( SDPCone  sdpcone,
double  mu,
DSDPSchurMat  M,
DSDPVec  vrhs1,
DSDPVec  vrhs2 
)

Compute the Hessian to the barrier term.

Parameters:
sdpcone cone
mu barrier parameter
M Schur matrix to insert elements.
vrhs1 dual objectvive gradient.
vrhs2 barrier gradient

Definition at line 30 of file sdpcompute.c.

References SDPCone_C::ATR, SDPCone_C::blk, SDPblk::bmu, DSDP_FALSE, DSDP_TRUE, DSDPBlockADot(), DSDPBlockGetMatrix(), DSDPBlockvAv(), DSDPDataMatGetEig(), DSDPDataMatGetRank(), DSDPDualMatInverseMultiply(), DSDPSchurMatAddRow(), DSDPSchurMatRowColumnScaling(), DSDPVMatAddOuterProduct(), DSDPVMatZeroEntries(), SDPblk::gammamu, DSDPDataTranspose::idA, SDPblk::IS, SDPblk::n, DSDPDataTranspose::nnzblocks, DSDPDataTranspose::nzblocks, SDPblk::S, SDPConeVecDot(), SDPblk::T, SDPblk::W, SDPblk::W2, SDPCone_C::Work, and SDPCone_C::Work2.

                                                                                                     {

  int i,k,kt,kk,m,n,rank,info;
  int ncols,ii,idA;
  double rtemp,ack,ggamma,bmu,scl;
  double rhs1i,rhs2i;
  DSDPTruth method1;
  SDPConeVec W,W2;
  DSDPVec MRowI=sdpcone->Work, Select=sdpcone->Work2;
  DSDPDataMat AA;
  DSDPDualMat S;
  DSDPVMat T;
  DSDPDataTranspose ATranspose=sdpcone->ATR;
  SDPblk *blk=sdpcone->blk;
  DSDPIndex IS;

  /* Evaluate M */
  DSDPFunctionBegin;
  SDPConeValid(sdpcone);
  info=DSDPVecGetSize(vrhs1,&m);DSDPCHKERR(info);

  for (i=0; i<m; i++){  /* One row at a time */

    /* Which Coluns */
    rhs1i=0;rhs2i=0;
    info=DSDPVecZero(MRowI);DSDPCHKERR(info);
    info=DSDPSchurMatRowColumnScaling(M,i,Select,&ncols); DSDPCHKERR(info); 
    if (ncols==0){continue;}
 
    for (kt=0; kt<ATranspose.nnzblocks[i]; kt++){ /* Loop over all blocks */
      kk=ATranspose.nzblocks[i][kt];
      idA=ATranspose.idA[i][kt];
      info=DSDPBlockGetMatrix(&blk[kk].ADATA,idA,&ii,&scl,&AA);DSDPCHKBLOCKERR(kk,info);
      if (ii!=i){DSDPSETERR2(8,"Data Transpose Error: var %d does not equal %d.\n",i,ii);}
      info = DSDPDataMatGetRank(AA,&rank,blk[kk].n);DSDPCHKBLOCKERR(kk,info);
      if (rank==0) continue;

      T=blk[kk].T; S=blk[kk].S; W=blk[kk].W; W2=blk[kk].W2;
      n=blk[kk].n; ggamma=blk[kk].gammamu; bmu=blk[kk].bmu; IS=blk[kk].IS;

      method1=DSDP_TRUE;  /* Simple heuristic */
      if (rank==1) method1=DSDP_FALSE;
      if (rank==2 && ncols<=n) method1=DSDP_FALSE;
      if (rank*rank*ncols<=n+1)method1=DSDP_FALSE;
      if (ncols*blk[kk].nnz < n*n/10) method1=DSDP_FALSE;
      if (ncols==1 && i==m-1)method1=DSDP_FALSE;
      if (n<5) method1=DSDP_TRUE;
      if (0==1) method1=DSDP_FALSE;
      if (method1==DSDP_TRUE){info=DSDPVMatZeroEntries(T);DSDPCHKBLOCKERR(kk,info);}
      for (k=0; k<rank; k++){
      
      info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKBLOCKERR(kk,info);
      if (ack==0.0) continue;
      ack*=scl;
      info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKBLOCKERR(kk,info);

      /* RHS terms */
      info = SDPConeVecDot(W,W2,&rtemp); DSDPCHKBLOCKERR(kk,info);
      if (rtemp==0.0) continue;
      rhs1i+=rtemp*ack*bmu; rhs2i+=rtemp*ack*ggamma*mu;
      ack*=(ggamma+bmu);

      if (method1==DSDP_TRUE){
        info=DSDPVMatAddOuterProduct(T,ack*mu,W2);DSDPCHKBLOCKERR(kk,info);
      } else {
        info=DSDPBlockvAv(&blk[kk].ADATA,ack*mu,Select,W2,MRowI);DSDPCHKBLOCKERR(kk,info);
      } /* End row computations for rank kk of block kk */
 
      }   /* End row computations for all of block kk     */

      if (method1==DSDP_TRUE){
      info=DSDPBlockADot(&blk[kk].ADATA,1.0,Select,T,MRowI);DSDPCHKBLOCKERR(kk,info);
      }   /* End row computations for all of block ll     */
    }     /* End row computations for all blocks          */
    info=DSDPVecAddElement(vrhs1,i,rhs1i);DSDPCHKERR(info);
    info=DSDPVecAddElement(vrhs2,i,rhs2i);DSDPCHKERR(info);
    info=DSDPSchurMatAddRow(M,i,1.0,MRowI);DSDPCHKERR(info);
  }

  DSDPFunctionReturn(0);
}


Generated by  Doxygen 1.6.0   Back to index