Residuals.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.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.tool.Tools;


/**
 * class <code>Residuals</code> compute the residual of primal-dual solution which is equal to the formula,<code>mu := X・Y</code>.
 * If you want to konow the more detail of implementing of algorithm, please refer to web site:
 *  <a href = "http://sdpa.sourceforge.net/" style=text-decoration:none>
 * <tt>http://sdpa.sourceforge.net/</tt></a>.<br>
 * 
 * @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 Residuals<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 s specified data type */
  private RS elementUnit;
  /**primal vector */
  private Vector<RS,RM,CS,CM> primalVec;
  /**dual matrix */
  private BlockDenseMatrix<RS,RM,CS,CM> dualMat;
  /**center norm */
  private RS centerNorm;
  /**norm of primal vector */
  private RS normPrimalVec;
  /**norm of dual matrix */
  private RS normDualMat;

  /**
   * 新しく生成された{@link Residuals}オブジェクトを初期化します。
   * 
   * @param elementUnit 単位数
   */
  public Residuals(RS elementUnit) {
    this.elementUnit = elementUnit;
    this.normPrimalVec = elementUnit.createZero();
    this.normDualMat = elementUnit.createZero();
    this.centerNorm = elementUnit.createZero();
    this.primalVec = new Vector<>(0, elementUnit.createZero());
    this.dualMat = new BlockDenseMatrix<>(elementUnit);
  }

  /**
   * constructor of class <code>Residuals</code>
   *<br> コンストラクタ
   * 
   * @param m dimension of vector
   * @param blockSize size information of block-matrices matrix
   * @param blockStruct structure information of block-matrices matrix
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code> := <code>F0</code>
   * @param A block-matrices matrix array <code>A</code> := <code>Fi</code>
   * @param currentPt current primal-dual solution
   */
  public Residuals(int m, int blockSize, int[] blockStruct, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, Solution<RS,RM,CS,CM> currentPt) {
    this(b.getElementUnit());
    initialize(m, blockSize, blockStruct, b, C, A, currentPt);
  }

  /**
   * initialize class <code>Residuals</code>. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(rsdpa_parts.cpp: SDPA 6.2.0 Line No. : 1399)</code>
   * 
   * @param m dimension of vector
   * @param blockSize size information of block-matrices matrix
   * @param blockStruct structure information of block-matrices matrix
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code> := <code>F0</code>
   * @param A block-matrices matrix array <code>A</code> := <code>Fi</code>
   * @param currentPt current primal-dual solution
   */
  private void initialize(int m, int blockSize, int[] blockStruct, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, Solution<RS,RM,CS,CM> currentPt) {
    this.primalVec.initialize(m);
    this.dualMat.initialize(blockSize, blockStruct);
    // p[k] = b[k] - A[k].X;
    for (int i = 0; i < m; ++i) {
      RS ip = Algebra.run(A[i], '.', currentPt.getXmat());
      this.primalVec.setElement(i, b.getElement(i).subtract(ip));
    }
    // D = C - Z - \sum A[k]y[k]
    this.dualMat = Algebra.let(C, '+', currentPt.getZmat(), this.elementUnit.create(-1));
    for (int i = 0; i < m; ++i) {
      this.dualMat = Algebra.let(this.dualMat, '-', A[i], currentPt.getYvec().getElement(i));
    }
    this.normPrimalVec = computeMaxNorm(this.primalVec);
    this.normDualMat = computeMaxNorm(this.dualMat);
  }

  /**
   * dualMatを返します。
   * 
   * @return dualMat
   */
  public BlockDenseMatrix<RS,RM,CS,CM> getDualMat() {
    return this.dualMat;
  }

  /**
   * primalVecを返します。
   * 
   * @return primalVec
   */
  public Vector<RS,RM,CS,CM> getPrimalVec() {
    return this.primalVec;
  }

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

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

  /**
   * compute the maximum norm of vector. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(<i>rsdpa_parts.cpp: SDPA 6.2.0 Line No. : 1430</i>)</code>
   * 
   * @param vec vector
   * @return max maximum norm of vector
   */
  private RS computeMaxNorm(Vector<RS,RM,CS,CM> vec) {
    final RS unit = vec.getElement(0);

    RS ret = unit.create(0);
    for (int i = 0; i < vec.getDimension(); ++i) {
      RS tmp = vec.getElement(i).abs();
      if (tmp.isGreaterThan(ret)) {
        ret = tmp;
      }
    }
    return ret;
  }

  /**
   * compute the maximum norm of dense block-matrices matrix. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(<i>rsdpa_parts.cpp: SDPA 6.2.0 Line No. : 1447</i>)</code>
   * 
   * @param mat block-matrices matrix
   * @return maximum norm of block-matrices matrix
   */
  private RS computeMaxNorm(BlockDenseMatrix<RS,RM,CS,CM> mat) {
    int blockSize = mat.getBlockSize();
    final RS unit = mat.getElementUnit();

    RS ret = unit.create(0);

    for (int i = 0; i < blockSize; ++i) {
      int struct = mat.getBlockStruct(i);
      if (struct < 0) {
        // case Diagonal
        final RS[] target = mat.getBlock(i).diagonalElements;
        for (int j = 0; j < -struct; ++j) {
          RS tmp = target[j].abs();
          if (tmp.isGreaterThan(ret)) {
            ret = tmp;
          }
        }
      } else {
        // case non Diagonal
        final RS[] target = mat.getBlock(i).denseElements;
        for (int j = 0; j < struct * struct; ++j) {
          RS tmp = target[j].abs();
          if (tmp.isGreaterThan(ret)) {
            ret = tmp;
          }
        }
      }
    }

    return ret;
  }

  /**
   * method <code>update</code> is to update the information of <code>mu</code> by each iteration. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(<i>rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 1498</i>)</code>
   * 
   * @param m dimension of vector
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code> := <code>F0</code>
   * @param A block-matrices matrix array <code>A</code> := <code>Fi</code>
   * @param currentPt current primal-dual solution
   * @param mu average complementarity <code>mu</code> := <code>X・Y</code>
   */
  public void update(int m, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, Solution<RS,RM,CS,CM> currentPt, AverageComplementarity<RS,RM,CS,CM> mu) {
    final RS unit = b.getElement(0);

    for (int i = 0; i < m; ++i) {
      RS ip = Algebra.run(A[i], '.', currentPt.getXmat());
      this.primalVec.setElement(i, b.getElement(i).subtract(ip));
    }
    this.normPrimalVec = computeMaxNorm(this.primalVec);
    this.dualMat = Algebra.let(C, '+', currentPt.getZmat(), unit.create(-1));
    for (int i = 0; i < m; ++i) {
      this.dualMat = Algebra.let(this.dualMat, '-', A[i], currentPt.getYvec().getElement(i));
    }
    this.normDualMat = computeMaxNorm(this.dualMat);
    this.centerNorm = unit.create(1).subtract(currentPt.xzMinEigenValue.divide(mu.getValue()));
  }

  /**
   * class <code>compute</code> calls the method <code>initialize</code>. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(<i>rsdpa_parts.cpp: SDPA 6.2.0 Line No. : 1542</i>)</code>
   * 
   * @param m dimension of vector
   * @param blockSize size information of block-matrices matrix
   * @param blockStruct structure information of block-matrices matrix
   * @param b vector <code>b</code> := <code>c</code>
   * @param C block-matrices matrix <code>C</code> := <code>F0</code>
   * @param A block-matrices matrix array <code>A</code> := <code>Fi</code>
   * @param currentPt current primal-dual solution
   */
  public void compute(int m, int blockSize, int[] blockStruct, Vector<RS,RM,CS,CM> b, BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, Solution<RS,RM,CS,CM> currentPt) {
    initialize(m, blockSize, blockStruct, b, C, A, currentPt);
  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {
    StringBuffer buff = new StringBuffer();
    buff.append(" currentRes.primalVec = \n"); //$NON-NLS-1$
    buff.append(this.primalVec.toString());
    buff.append(" currentRes.dualMat = \n"); //$NON-NLS-1$
    buff.append(this.dualMat.toString());
    buff.append(" currentRes.normPrimalVec =\n"); //$NON-NLS-1$
    buff.append(this.normPrimalVec.toString() + "\n"); //$NON-NLS-1$
    buff.append(" currentRes.normDualMat =\n"); //$NON-NLS-1$
    buff.append(this.normDualMat.toString());
    return buff.toString();
  }

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

  /**
   * create a copy of class <i><code>Residuals</code></i>
   * 
   * @return a copy
   */
  public Residuals<RS,RM,CS,CM> createClone() {
    Residuals<RS,RM,CS,CM> inst = new Residuals<>(this.elementUnit);
    inst.elementUnit = this.elementUnit == null ? null : this.elementUnit;
    inst.primalVec = this.primalVec == null ? null : (Vector<RS,RM,CS,CM>)this.primalVec.createClone();
    inst.dualMat = this.dualMat == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.dualMat.createClone();
    inst.normPrimalVec = this.normPrimalVec == null ? null : this.normPrimalVec;
    inst.normDualMat = this.normDualMat == null ? null : this.normDualMat;
    inst.centerNorm = this.centerNorm == null ? null : this.centerNorm;
    return inst;
  }
}