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

int DSDPComputeDualStepDirections ( DSDP  dsdp  ) 

Compute the step direction by computing a linear system and solving it.

Parameters:
dsdp the solver
DSDP first attempts unpreconditioned CG to the matrix. Once the number of iterations becomes too large, it swithes a CG preconditioned by the Cholesky factorization. Usually only one iteration of the preconditioned CG is necessary, but solutions with large norms and very precise solutions may require additional iterations.

Definition at line 370 of file dualalg.c.

References DSDP_FALSE, DSDP_INDEFINITE_SCHUR_MATRIX, DSDP_TRUE, DSDPCGSolve(), DSDPComputeG(), DSDPComputeHessian(), DSDPGetMaxYElement(), DSDPInvertS(), DSDPSchurMatFactor(), DSDPSchurMatShiftDiagonal(), DSDPSchurMatView(), and DSDPSetConvergenceFlag().

Referenced by DSDPSolveDynamicRho().

                                            {
  int info,computem=1;
  double madd,ymax,cgtol=1e-7;
  DSDPTruth cg1,cg2,psdefinite;
  DSDPFunctionBegin;

  if (dsdp->itnow>30) dsdp->slestype=3;
  if (dsdp->rgap<1e-3) dsdp->slestype=3;
  if (dsdp->m<40) dsdp->slestype=3;
  if (0 && dsdp->itnow>20 && dsdp->m<500) dsdp->slestype=3;
  info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
  if (dsdp->slestype==1){
    cg1=DSDP_TRUE; cg2=DSDP_TRUE;
    info=DSDPInvertS(dsdp);DSDPCHKERR(info);
    info=DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
    info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
    if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
    if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=2;
  }
  if (dsdp->slestype==2){
    cg1=DSDP_TRUE; cg2=DSDP_TRUE;
    DSDPLogInfo(0,9,"Compute Hessian\n");
    info=DSDPInvertS(dsdp);DSDPCHKERR(info);
    info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
    computem=0;
    DSDPLogInfo(0,9,"Apply CG\n");
    info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
    if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
    if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=3;
    
  }
  if (dsdp->slestype==3){
    DSDPLogInfo(0,9,"Factor Hessian\n");
    psdefinite=DSDP_FALSE;
    if (dsdp->Mshift < 1e-12 || dsdp->rgap<0.1 || dsdp->Mshift > 1e-6){
      madd=dsdp->Mshift;
    } else {
      madd=1e-13;
    }
    if (computem){
      info=DSDPInvertS(dsdp);DSDPCHKERR(info);
    }
    while (psdefinite==DSDP_FALSE){
      if (0==1 && dsdp->Mshift>dsdp->maxschurshift){ 
      info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);  
      break;
      }
      if (0 && dsdp->Mshift*ymax>dsdp->pinfeastol/10){ 
      info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);  
      break;
      }
      if (madd*ymax>dsdp->pinfeastol*1000){ 
      info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);  
      break;
      }
      if (computem){
      info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);}
      if (0==1){info=DSDPSchurMatView(dsdp->M);DSDPCHKERR(info);}
      info = DSDPSchurMatShiftDiagonal(dsdp->M,madd);DSDPCHKERR(info);
      info = DSDPSchurMatFactor(dsdp->M,&psdefinite); DSDPCHKERR(info);
      computem=1;
      if (psdefinite==DSDP_FALSE){ 
      madd=madd*4 + 1.0e-13;
      }
    }
    dsdp->Mshift=madd;
    if (psdefinite==DSDP_TRUE){
      info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
      info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);
    }
  }
  
  DSDPFunctionReturn(0);
}


Generated by  Doxygen 1.6.0   Back to index