HsolveMethod.java

/*
 * Created on 2004/10/18
 *
 * $Id: HsolveMethod.java,v 1.23 2008/02/02 03:07:13 koga Exp $
 * Copyringht (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.cga.linear;

import org.mklab.cga.interval.matrix.IntervalNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalNumericalScalar;
import org.mklab.cga.util.IntervalUtil;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 線形方程式の精度保証付き解をHansen, Bliek, Rohn, Ning, Kearfott, Neumaierの方法で求めるクラスです。
 * 
 * @author hiroki
 * @version $Revision: 1.23 $. 2004/10/18
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class HsolveMethod<IS extends IntervalNumericalScalar<IS,IM,S,M>, IM extends IntervalNumericalMatrix<IS,IM,S,M>, S extends NumericalScalar<S,M>, M extends NumericalMatrix<S,M>> implements LinearEquationVerifier<IS,IM,S,M> {

  /** 係数区間行列 */
  private IM intA;

  /** 区間右辺ベクトル */
  private IM intB;

  /** 精度保証付き解 */
  private IM solution;

//  /**
//   * 新しく生成された<code>HsolveMethod</code>オブジェクトを初期化します。
//   */
//  public HsolveMethod() {
//    //
//  }
//
//  /**
//   * 新しく生成された<code>HsolveMethod</code>オブジェクトを初期化します。
//   * 
//   * @param A 係数行列
//   * @param b 右辺ベクトル
//   */
//  public HsolveMethod(final M A, final M b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = IntervalMatrixFactory.toInterval(b);
//  }
//
//  /**
//   * 新しく生成された<code>HsolveMethod</code>オブジェクトを初期化します。
//   * 
//   * @param A 係数行列
//   * @param b 右辺区間ベクトル
//   */
//  public HsolveMethod(final M A, final IM b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = b.clone();
//  }
//
//  /**
//   * 新しく生成された<code>HsolveMethod</code>オブジェクトを初期化します。
//   * 
//   * @param A 係数区間行列
//   * @param b 右辺ベクトル
//   */
//  public HsolveMethod(final IM A, final M b) {
//    this.intA = A.clone();
//    this.intB = IntervalMatrixFactory.toInterval(b);
//  }


//
//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#solve(org.mklab.nfc.matrix.NumericalMatrix, org.mklab.nfc.matrix.NumericalMatrix)
//   */
//  public void solve(final M A, final M b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = IntervalMatrixFactory.toInterval(b);
//    solve();
//  }
//
//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#solve(org.mklab.nfc.matrix.NumericalMatrix, org.mklab.cga.interval.matrix.IntervalMatrix)
//   */
//  public void solve(final M A, final IM b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = b.clone();
//    solve();
//  }
//
//  /**
//   * {@inheritDoc}
//   */
//  public void solve(final IM A, final M b) {
//    this.intA = A.clone();
//    this.intB = IntervalMatrixFactory.toInterval(b);
//    solve();
//  }
  
/**
 * 新しく生成された<code>HsolveMethod</code>オブジェクトを初期化します。
 * 
 * @param A 係数区間行列
 * @param b 右辺区間ベクトル
 */
public HsolveMethod(final IM A, final IM b) {
  this.intA = A;
  this.intB = b;
}

  /**
   * Creates interval matrix
   * 
   * @param matrix matrix
   * @return interval matrix
   */
  private IM createInterval(M matrix) {
    return this.intA.create(matrix);
  }

  /**
   * Creates interval with middle and radius.
   * 
   * @param mid middle
   * @param rad radius
   * @return interval matrxi
   */
  private IM createMidRad(M mid, M rad) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    M inf = mid.subtract(rad);
    
    manager.setRoundMode(RoundMode.ROUND_UP);
    M sup = mid.add(rad);
    
    IM ans = this.intA.createInfSup(inf, sup);
    
    manager.setRoundMode(oldRoundMode);
    
    return ans;
  }
  
  /**
   * {@inheritDoc}
   */
  public void solve() {
//    this.intA = A;
//    this.intB = b;

    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    M realMatrix = this.intA.getInfimum();
    int n = this.intA.getColumnSize();

    M midA = this.intA.getMiddle();

    IM C = createInterval(midA.inverse());
    IM AA = C.multiply(this.intA);// C*A
    IM bb = C.multiply(this.intB);// C*b
    IM dA = AA.diagonalToVector();// diag(A)
    M compA = IntervalUtil.comparisonMatrix(AA);// compmat(A)
    M B = compA.inverse();// A.inverse()
    M v = (B.multiply(realMatrix.createOnes(n, 1))).absElementWise();// abss(B*ones(n,1))

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    M u = compA.multiply(v);// A*v
    M uMinRowWise = (u).minRowWise();

    boolean isAaHmatrix = uMinRowWise.compareElementWise(".>", realMatrix.createZero(uMinRowWise.getRowSize(), uMinRowWise.getColumnSize())).allTrue(); //$NON-NLS-1$

    if (isAaHmatrix == false) {
      throw new RuntimeException("A is not an H-matrix"); //$NON-NLS-1$
    }

    M dAc = compA.diagonalToVector();// diag(A)
    M AAA = compA.multiply(B).subtract(realMatrix.createUnit(n));// A*B-I

    manager.setRoundMode(RoundMode.ROUND_UP);

    M w = realMatrix.createZero(1, n);// zeros(1, n)

    for (int i = 0; i < n; i++) {
      w = (w).maxElementWise(AAA.getColumnVector(i + 1).transpose().unaryMinus().divide((u).getElement(i + 1)));
    }

    M dlow = v.multiplyElementWise(w.conjugateTranspose()).subtract(B.diagonalToVector());
    dlow = dlow.unaryMinus();
    B = B.add(v.multiply(w));
    u = B.multiply(bb.abssElementWise());
    M d = B.diagonalToVector();

    M alpha = dAc.add(realMatrix.createOnes(d.getRowSize(), d.getColumnSize()).unaryMinus().divideElementWise(d));
    M beta = u.divideElementWise(dlow).subtract(bb.abssElementWise());

    M mid1 = realMatrix.createZero(beta.getRowSize(), beta.getColumnSize());
    M mid2 = realMatrix.createZero(alpha.getRowSize(), alpha.getColumnSize());
    IM x = (bb.add(createMidRad(mid1, beta))).divideElementWise(dA.add(createMidRad(mid2, alpha)));

    this.solution = x;

    manager.setRoundMode(oldRoundMode);
  }

  /**
   * {@inheritDoc}
   */
  public IM getSolution() {
    return this.solution;
  }

//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#getSolution()
//   */
//  public LinearEquationVerifiedSolution<E> getSolution() {
//    return new LinearEquationVerifiedSolution(this.solution);
//  }
}