SolveInfo.java

package org.mklab.sdpj.algorithm;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.algebra.Algebra;
import org.mklab.sdpj.algebra.AlgebraWithEFFloat;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.tool.Tools;


/**
 * class <code>SolveInfo</code> save the solution information of solving SDP problem. 
 * 
 * @author koga
 * @version $Revision$, 2009/04/24
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
 */
public class SolveInfo<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {

  /**unit based on a generic data type */
  private RS unit;
  /**important variable <code>rho</code> */
  private RS rho;
  /**theta of primal problem */
  private RS etaPrimal;
  /**theta of dual problem */
  private RS etaDual;
  /**intermediate primal solution */
  private RS objValPrimal;
  /**intermediate dual solution */
  private RS objValDual;

  /**
   * constructor of class <code>SolveInfo</code><br>
   * コンストラクタ
   * 
   * @param nDim dimension of block-matrices matrix
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code>
   * @param A block-matrices matrix array <code>A</code>
   * @param initPt initial primal-dual solution
   * @param mu0 <code>mu</code> := <code>X・Y</code>
   * @param omegaStar value of parameter <code>omegaStar</code>
   */
  public SolveInfo(int nDim, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, Solution<RS,RM,CS,CM> initPt, RS mu0, RS omegaStar) {
    this.unit = mu0.create(1);

    this.rho = create(0);
    this.etaPrimal = create(0);
    this.etaDual = create(0);
    this.objValPrimal = create(0);
    this.objValDual = create(0);

    this.rho = create(1);
    this.etaPrimal = omegaStar.multiply(nDim).multiply(mu0);
    this.etaDual = omegaStar.multiply(nDim).multiply(mu0);
    this.objValPrimal = AlgebraWithEFFloat.run(C, '.', initPt.getXmat());
    this.objValDual = AlgebraWithEFFloat.run(b, '.', initPt.getYvec());
//    this.objValPrimal = Algebra.run(C, '.', initPt.getXmat());
//    this.objValDual = Algebra.run(b, '.', initPt.getYvec());
  }

  /**
   * rhoを返します。
   * 
   * @return rho
   */
  public RS getRho() {
    return this.rho;
  }

  /**
   * objValPrimalを返します。
   * 
   * @return objValPrimal
   */
  public RS getObjValPrimal() {
    return this.objValPrimal; 
  }

  /**
   * objValDualを返します。
   * 
   * @return objValDual
   */
  public RS getObjValDual() {
    return this.objValDual; 
  }

  /**
   * update the information of computing each step, which is ready for next iteration of computing.
   * 
   * @param nDim dimension of block-matrices matrix
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code>
   * @param initPt initial primal-dual solution
   * @param currentPt current primal-dual solution
   * @param currentRes current primal-dual residual error
   * @param mu <code>mu</code> := <code>X・Y</code>
   * @param theta value of theta
   * @param param setting parameter of solving SDP problem
   */
  public void update(int nDim, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, Solution<RS,RM,CS,CM> initPt, Solution<RS,RM,CS,CM> currentPt, Residuals<RS,RM,CS,CM> currentRes, AverageComplementarity<RS,RM,CS,CM> mu, RatioInitResCurrentRes<RS,RM,CS,CM> theta,
      Parameter<RS,RM,CS,CM> param) {
    this.objValPrimal = AlgebraWithEFFloat.run(C, '.', currentPt.getXmat());
    this.objValDual = AlgebraWithEFFloat.run(b, '.', currentPt.getYvec());
//    this.objValPrimal = Algebra.run(C, '.', currentPt.getXmat());
//    this.objValDual = Algebra.run(b, '.', currentPt.getYvec());

    RS primal = theta.getPrimal();
    RS dual = theta.getDual();
    RS omega = param.omegaStar;
    this.rho = createZero();// 0.0;0.0;
    RS x0z0 = mu.getInitialValue().multiply(nDim);
    RS xMatzMat = mu.getValue().multiply(nDim);
    RS x0zMat = Algebra.run(initPt.getXmat(), '.', currentPt.getZmat());
    RS xMatz0 = Algebra.run(currentPt.getXmat(), '.', initPt.getZmat());

    RS accuracy = param.epsilonDash;

    if (currentRes.getNormPrimalVec().isLessThanOrEquals(accuracy)) {
      // rMessage("primal accuracy");
      if (xMatz0.isLessThan(this.etaPrimal)) {
        this.etaPrimal = xMatz0;
      }
    }
    if (currentRes.getNormDualMat().isLessThanOrEquals(accuracy)) {
      // rMessage("dual accuracy");
      if (x0zMat.isLessThan(this.etaDual)) {
        this.etaDual = x0zMat;
      }
    }

    // primal is infeasible and dual is feasible
    if (currentRes.getNormPrimalVec().isGreaterThan(accuracy) && currentRes.getNormDualMat().isLessThanOrEquals(accuracy)) {
      this.rho = (primal.multiply(x0zMat)).divide((((primal.add((createUnit().subtract(primal)).multiply(omega))).multiply(this.etaDual)).add(xMatzMat)));
    }

    // primal is feasible and dual is infeasible
    if (currentRes.getNormPrimalVec().isLessThanOrEquals(accuracy) && currentRes.getNormDualMat().isGreaterThan(accuracy)) {
      this.rho = (dual.multiply(xMatz0)).divide((((dual.add((createUnit().subtract(dual)).multiply(omega))).multiply(this.etaPrimal)).add(xMatzMat)));
    }

    // primal and dual are infeasible
    if (currentRes.getNormPrimalVec().isGreaterThan(accuracy) && currentRes.getNormDualMat().isGreaterThan(accuracy)) {
      this.rho = ((dual.multiply(xMatz0)).add(primal.multiply(x0zMat))).divide(((((primal.multiply(dual)).add(omega.multiply(((primal.multiply((createUnit().subtract(dual)))).add((createUnit()
          .subtract(primal)).multiply(dual)))))).multiply(x0z0)).add(xMatzMat)));
    }
  }

  /**
   * show messages.
   */
  public void display() {
    Tools.message(this.toString());
  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {
    StringBuffer buff = new StringBuffer();
    buff.append("rSolveInfo.rho          = " + this.rho.toString() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("rSolveInfo.etaPrimal    = " + this.etaPrimal.toString() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("rSolveInfo.etaDual      = " + this.etaDual.toString() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("rSolveInfo.objValPrimal = " + this.objValPrimal.toString() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("rSolveInfo.objValDual   = " + this.objValDual.toString() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    return buff.toString();
  }

  /**
   * 単位数を返します。
   * 
   * @return 単位数
   */
  private RS createUnit() {
    return this.unit.createUnit();
  }

  /**
   * 零を返します。
   * 
   * @return 零を
   */
  private RS createZero() {
    return this.unit.createZero();
  }

  /**
   * int型に対応する値を返します。
   * 
   * @param value int型の値
   * @return E型の値
   */
  private RS create(int value) {
    return this.unit.create(value);
  }
}