SequentialMethod.java

/*
 * Created on 2004/01/09
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.cga.eigen;

import org.mklab.cga.interval.matrix.IntervalComplexNumericalMatrix;
import org.mklab.cga.interval.matrix.IntervalRealNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalComplexNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalRealNumericalScalar;
import org.mklab.nfc.eig.EigenSolution;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 実対象行列の固有ベクトルの誤差を求めるクラスです。 未完成。
 * 
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 * @author Hiroki
 */
public class SequentialMethod<RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, 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>> {

  /** 任意の行列 */
  private RM matrix;

  /** 任意の行列 */
  private RM A;

  /** 近似固有値 */
  private CM val;

  /** 近似固有ベクトル */
  private CM vec;

  /** 近似固有値と近似固有ベクトル */
  private EigenSolution<RS,RM,CS,CM> eigen;

  /**
   * 固有値を求めたい行列を設定します。
   * 
   * @param matrix 固有値を求めたい任意の行列。
   */
  public SequentialMethod(final RM matrix) {
    this.matrix = matrix;
    this.eigen = this.A.eigenDecompose();
    this.val = this.eigen.getValue();
    this.vec = this.eigen.getVector();
  }

  /**
   * Bauer-Fike型の事前誤差評価により対角化可能行列の精度保証付き固有値を求めます。
   * 
   * @return 近似固有値と誤差をMatrix型の配列で返す 。 第一要素は近似固有値 、 第二要素は誤差 。
   */
  // public Matrix[] eigenValueError() {
  //
  // Matrix e_max, f_max, error;
  // IntervalMatrix i_vec, i_e, i_f;
  // double e_norm, f_norm, q_norm, q_inv_norm;
  //
  // /** ここから QTQ'=matrix+E, Q'Q=I+F を用いる。 **/
  // /** Qは近似固有ベクトルvec、Tは近似固有値valに相当します。 **/
  // /** 上の式は E=vec*val*vec'-matrix, F=vec'*vec-I に変形できます。 **/
  // /** さらにEとFを求めるのに区間演算の考え(半径と中心の考え)を **/
  // /** 利用して計算を行います。 **/
  //
  // i_vec = new RealIntervalMatrix(vec);
  //
  // i_e = i_vec.multiply(val).multiply(i_vec.inverse()).subtract(matrix).abs();
  // i_f =
  // i_vec.multiply(i_vec.inverse()).subtract(Matrix.I(matrix.getColSize())).abs();
  // /** EとFの計算ここまで。 **/
  //
  // /** 計算したEとFに関して上限と下限の絶対値の大きいほうを **/
  // /** 選び、EとFが最大になるようにします。 **/
  //
  // //processManager.calcPrepare();
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  // e_max = i_e.getUpper().maxElementWise(i_e.getLower());
  // f_max = i_f.getUpper().maxElementWise(i_f.getLower());
  //
  // /** ノルムの計算。それぞれの値ができるだけ誤差を大きくするように **/
  // /** 丸めモードを切り替えます。**/
  // q_norm = vec.norm(2);
  // q_inv_norm = vec.inverse().norm(2);
  // e_norm = e_max.norm(2);
  //
  // double matrix_norm = matrix.norm(2);
  //
  // FPU.setRoundMode(FPU.ROUND_DOWN);
  // f_norm = f_max.norm(2);
  // /** ノルムの計算ここまで。 **/
  //
  // /** 誤差を計算します。より大きな値になるように上向きに丸めて計算します。**/
  // FPU.setRoundMode(FPU.ROUND_UP);
  //
  // error = new RealMatrix(new double[] {((q_norm * q_inv_norm) / (1 - f_norm))
  // * e_norm });
  //
  // Matrix error2 = new RealMatrix(new double[] { e_norm + matrix_norm * f_norm
  // });
  //
  // FPU.setRoundMode(FPU.ROUND_NEAR);
  //
  // //processManager.calcComplete();
  //
  // /** 戻り値はMatrixの配列。第一要素は近似固有値val, 第二要素は誤差。**/
  // Matrix[] answer = { val, error, error2 };
  // return answer;
  // }

  /**
   * 実対称行列の全ての固有値に対する誤差を求めます。
   * 
   * @return 近似固有値と誤差をMatrix型の配列で返します。第一要素は近似固有値、第二要素は誤差。
   */
  // public Matrix[] eigenValueErrorOfRSM() {
  // Matrix e_max, f_max, preans, ans;
  // double e_norm, f_norm;
  // int matrix_size;
  //
  // /** 近似固有値と近似固有ベクトルの計算ここまで。**/
  //
  // Matrix[] ef = efCalc(val, vec);
  //
  // e_max = ef[0];
  // f_max = ef[1];
  //
  // preans = val.diagToVector();
  // ans = preans.getRealPart();
  //
  // matrix_size = matrix.getColSize();
  //
  // //processManager.calcPrepare();
  //
  // FPU.setRoundMode(FPU.ROUND_DOWN);
  //
  // f_norm = f_max.norm("Inf");
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  //
  // e_norm = e_max.norm("fro");
  //
  // Matrix error = new RealMatrix(new double[] {((1 + f_norm) / ((1 - f_norm) *
  // (1 - Math.sqrt(matrix_size) * f_norm))) * e_norm });
  //
  // Matrix[] answer = { ans, error };
  //
  // FPU.setRoundMode(FPU.ROUND_NEAR);
  //
  // //processManager.calcComplete();
  //
  // return answer;
  //
  // }

  /**
   * 実対称行列のぞれぞれの固有値の誤差を求めます。
   * 
   * @return 近似固有値と誤差をMatrix型の配列で返します。第一要素は近似固有値、第二要素は誤差。
   */
  // public Matrix[] eachEigenValueErrorOfRSM() {
  // Matrix e_max, f_max, preans, ans, error;
  // double e_norm, f_norm;
  //
  // Matrix[] ef = efCalc(val, vec);
  //
  // e_max = ef[0];
  // f_max = ef[1];
  //
  // ans = val.diagToVector().getRealPart();
  //
  // //processManager.calcPrepare();
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  // e_norm = e_max.norm("Inf");
  // f_norm = f_max.norm("Inf");
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  //
  // error =
  // val.absElementWise().multiply(f_norm).add(Matrix.I(matrix.getColSize()).multiply(e_norm)).diagToVector();
  //
  // Matrix[] answer = { ans, error };
  //
  // FPU.setRoundMode(FPU.ROUND_NEAR);
  //
  // //processManager.calcComplete();
  //
  // return answer;
  //
  // }

  /**
   * Bauer-Fike型の事後誤差評価により対角化可能行列の精度保証付き固有値を求めます。 事前に誤差を求めたい行列が正規かどうか判別できれば計算できるらしい。未完成。
   * 
   * @return 近似固有値と誤差をMatrix型の配列で返します。第一要素は近似固有値、第二要素は誤差。
   */
  // public Matrix[] eigenErrorWithBF() {
  //
  // Matrix r_max;
  // double r_norm, x_norm, error;
  // RealIntervalMatrix i_matrix, i_r, i_val;
  // Matrix matrix_ct = matrix.conjugate().transpose();
  // int matrix_rowSize = matrix.getRowSize();
  // Matrix valErrors = new ComplexMatrix(matrix_rowSize, 1);
  // Matrix val_error = null;
  // Complex rambda;
  // Matrix vector;
  // RealIntervalMatrix[] i_vector;
  //
  // if (matrix.multiply(matrix_ct).equals(matrix_ct.multiply(matrix))) {
  // i_matrix = new RealIntervalMatrix(matrix);
  //
  // for (int i = 1; i <= matrix_rowSize; i++) {
  // rambda = val.getComplexElement(i, i);
  // vector = vec.getSubMatrix(1, matrix_rowSize, i, 1);
  // //i_vector[i] = new RealIntervalMatrix(vector);
  // //未完成。
  // i_r = i_matrix.multiply(vector).subtract(vector.multiply(rambda)).abs();
  //
  // //processManager.calcPrepare();
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  // r_max = i_r.getUpper().maxElementWise(i_r.getLower());
  //
  // r_norm = r_max.norm(2);
  //
  // FPU.setRoundMode(FPU.ROUND_DOWN);
  // x_norm = vec.norm(2);
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  // error = r_norm / x_norm;
  //
  // FPU.setRoundMode(FPU.ROUND_NEAR);
  // //processManager.calcComplete();
  // valErrors.setElement(i, 1, error);
  // }
  //
  // } else {
  // System.out.println("正規行列を設定してください");
  // }
  //
  // val_error = new ComplexMatrix(new Complex[] {((ComplexMatrix)
  // valErrors).max()});
  //
  // Matrix[] answer = { val, valErrors, val_error };
  //
  // return answer;
  // }

  /**
   * 実対称行列の固有ベクトルの誤差を求めます。未完成。 よくわからない。
   * 
   * @return 近似解と誤差。
   */
  public RM[] eigenVectorErrorOfRSM() {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    RM r_down, r_up, r_max, vector;
    RS r_norm, rambda, rate;
    int matrix_rowSize;

    RM[] eigenValueError = eigenValueErrorOfRSM();

    RS valueError = (eigenValueError[1]).getElement(1, 1);
    matrix_rowSize = this.matrix.getRowSize();

    RS[] error =  valueError.createArray(matrix_rowSize);
    RM vecError =  valueError.createZeroGrid(matrix_rowSize, 1);

    for (int i = 1; i <= matrix_rowSize; i++) {
      //System.out.println(i);

      manager.setRoundMode(RoundMode.ROUND_NEAR);

      rambda = (this.val.getRealPart()).getElement(i, i);
      System.out.println(rambda);
      vector = this.vec.getRealPart().getSubMatrix(1, matrix_rowSize, i, i);

      manager.setRoundMode(RoundMode.ROUND_DOWN);

      r_down = (this.matrix.multiply(vector).subtract(vector.multiply(rambda))).absElementWise();

      manager.setRoundMode(RoundMode.ROUND_UP);

      r_up = (this.matrix.multiply(vector).subtract(vector.multiply(rambda))).absElementWise();

      r_max = (r_down).maxElementWise(r_up);

      r_norm = (r_max).norm(NormType.TWO);

      rate = r_norm.divide(valueError);

      System.out.println("r_norm/valueError" + rate); //$NON-NLS-1$
      error[i - 1] = rate.multiply(rate).add(1).sqrt().multiply(rate);

      //manager.setRoundMode(RoundMode.ROUND_NEAR);

      vecError.setElement(i, 1, error[i - 1]);
    }

    RM[] answer = (RM[])new RealNumericalMatrix[2];
    answer[0] = (this.vec).getRealPart();
    answer[1] = vecError;

    manager.setRoundMode(oldRoundMode);

    return answer;
  }

  /**
   * @return ?
   */
  public CM[] eigenValueError() {
    // TODO Auto-generated method stub
    return null;
  }

  /**
   * @return ?
   */
  public RM[] eigenValueErrorOfRSM() {
    // TODO Auto-generated method stub
    return null;
  }

  /**
   * @return ?
   */
  public CM[] eachEigenValueErrorOfRSM() {
    // TODO Auto-generated method stub
    return null;
  }

  /**
   * @return ?
   */
  public CM[] eigenErrorWithBF() {
    // TODO Auto-generated method stub
    return null;
  }

  /**
   * 実対称行列の固有値の誤差を求めるのに必要な計算を行います。
   * 
   * @param val 近似固有値。
   * @param vec 近似固有ベクトル。
   * @return Matrixの配列。第一引数はQTQ'=matrix+EのEの値、第二引数はQ'Q=I+FのFの値。
   */
  // private Matrix[] efCalc(final Matrix val, final Matrix vec) {
  // Matrix e_max, f_max;
  // RealIntervalMatrix i_vec, i_e, i_f;
  //
  // i_vec = new RealIntervalMatrix(vec);
  // i_e = i_vec.multiply(val).multiply(vec.transpose()).subtract(matrix).abs();
  // i_f =
  // i_vec.multiply(vec.transpose()).subtract(Matrix.I(matrix.getColSize())).abs();
  //
  // //processManager.calcPrepare();
  //
  // FPU.setRoundMode(FPU.ROUND_UP);
  //
  // e_max = i_e.getUpper().maxElementWise(i_e.getLower());
  // f_max = i_f.getUpper().maxElementWise(i_f.getLower());
  //
  // Matrix[] ef_max = { e_max, f_max };
  //
  // FPU.setRoundMode(FPU.ROUND_NEAR);
  //
  // //processManager.calcComplete();
  // return ef_max;
  // }
}