YamamotoMethod.java

/*
 * Created on 2004/10/18
 *
 * $Id: YamamotoMethod.java,v 1.21 2008/02/02 03:07:13 koga Exp $
 * Copyright (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.nfc.matrix.NormType;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 線形方程式の精度保証付き解を山本の定理の方法で求めるクラスです。
 * 
 * @author hiroki
 * @version $Revision: 1.21 $. 2004/10/18
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class YamamotoMethod<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 M _A;
  /** 右辺係数ベクトル */
  private M _b;

  /** 係数行列 */
  private IM _IA;
  /** 右辺係数ベクトル */
  private IM _Ib;

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

//  /**
//   * コンストラクタ
//   * 
//   */
//  public YamamotoMethod() {
//    //
//  }

  /**
   * 新しく生成された<code>YamamotoMethod</code>オブジェクトを初期化します。
   * 
   * @param A 係数行列
   * @param b 右辺係数ベクトル
   */
  public YamamotoMethod(final IM A, final IM b) {
    this._IA = A;
    this._Ib = 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._A = A.clone();
//    this._b = b.clone();
//    solve();
//  }


//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#solve(org.mklab.cga.interval.matrix.IntervalMatrix, org.mklab.nfc.matrix.NumericalMatrix)
//   */
//  public void solve(final IM A, final M b) {
//    throw new RuntimeException("Not implemented"); //$NON-NLS-1$
//  }
//
//  /**
//   * @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) {
//    throw new RuntimeException("Not implemented"); //$NON-NLS-1$
//  }

  /**
   * @see org.mklab.cga.linear.LinearEquationVerifier#getSolution()
   */
  public IM getSolution() {
    return this.solution;
  }

  /**
   * {@inheritDoc}
   */
  public void solve() {
    //this._IA = IA;
    this._A = this._IA.getMiddle();
    //this._Ib = Ib;
    this._b = this._Ib.getMiddle();

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

//    String message = Abcdchk.abcdchk(this._A, this._b);
//    if (message.length() > 0) {
//      throw new IllegalArgumentException(message);
//    }

    M R = this._A.inverse();
    M x = R.multiply(this._b);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    int length = this._A.getColumnSize();
    M Gd = ((R.multiply(this._A).subtract(this._A.createUnit(length)))).absElementWise();
    M rd = this._A.multiply(x).subtract(this._b);

    manager.setRoundMode(RoundMode.ROUND_UP);

    M Gu = ((R.multiply(this._A).subtract(this._A.createUnit(length)))).absElementWise();
    M ru = this._A.multiply(x).subtract(this._b);
    M center = (ru.add(rd)).divide(2);
    M radius = center.subtract(rd);

    // 成分毎に値が大きな方を取得
    M G_max = (Gd).maxElementWise(Gu);
    S G_norm = ((G_max).norm(NormType.INFINITY));
    M Rru = ((R.multiply(center).add((R).absElementWise().multiply(radius)))).absElementWise();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    M Rrd = ((R.multiply(center).add((R).absElementWise().multiply(radius)))).absElementWise();
    S D = G_norm.unaryMinus().add(1);

    if (D.isLessThanOrEquals(0)) {
      manager.setRoundMode(oldRoundMode);
      throw new RuntimeException("false"); //$NON-NLS-1$
    }

    manager.setRoundMode(RoundMode.ROUND_UP);

    M Rr = ((Rru)).maxElementWise(Rrd);
    S Rr_norm = Rr.norm(NormType.INFINITY);
    S error = Rr_norm.divide(D);

    M kappa = G_max.multiply(this._A.createOnes(length, 1));
    M error2 = Rr.add(kappa.multiply(error));

    manager.setRoundMode(oldRoundMode);
    //this.solution = new RealIntervalMatrix(x, error2, "midrad"); //$NON-NLS-1$
    M inf = inf(x, error2);
    M sup = sup(x, error2);
    
    this.solution = this._IA.createInfSup(inf, sup);
  }

  /**
   * Returns supremum.
   * 
   * @param mid middle
   * @param rad radius
   * @return supremum
   */
  private M sup(M mid, M rad) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();
    
    manager.setRoundMode(RoundMode.ROUND_UP);
    
    M sup = mid.add(rad);
    
    manager.setRoundMode(oldRoundMode);
    return sup;
  }
  
  /**
   * Returns infimum.
   * 
   * @param mid middle
   * @param rad radius
   * @return infimum
   */
  private M inf(M mid, M rad) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();
    
    manager.setRoundMode(RoundMode.ROUND_DOWN);
    
    M inf = mid.subtract(rad);
    
    manager.setRoundMode(oldRoundMode);
    return inf;
  }

}