RealOhishiMethod.java

/*
 * Created on 2004/11/12
 * 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;


/**
 * 実対称行列の固有値の精度保証を行うクラスです。
 * 
 * @author hiroki
 * @version $Revision: 1.27 $.2004/11/12
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class RealOhishiMethod<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>> implements RealEigenVerifier<RIS,RIM,CIS,CIM,RS,RM,CS,CM> {

  /** 対象となる行列 */
  private RM A;
  /** 対象となる区間行列 */
  private RIM IA;
  /** 精度保証付き固有値 */
  private CIM verifiedValue;

  //  /** 精度保証付き固有ベクトル */
  //  private CIM verifiedVector;

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

//  /**
//   * 新しく生成された<code>OhishiMethod</code>オブジェクトを初期化します。
//   * 
//   * @param A matrix
//   * @param IA interval matrix
//   */
//  public RealOhishiMethod(RM A, RIM IA) {
//    this.A = A;
//    this.IA = IA;
//  }

  /**
   * 新しく生成された<code>OhishiMethod</code>オブジェクトを初期化します。
   * 
   * @param IA 対象となる区間行列
   */
  public RealOhishiMethod(RIM IA) {
    this.IA = IA;
    this.A = IA.getMiddle();
  }
//
//  /**
//   * {@inheritDoc}
//   */
//  public void solve() {    EigenSolution<RS,RM,CS,CM> eigen = this.A.eigenDecompose();
//
//  CM P1 = eigen.getVectors();
//  CM D1 = eigen.getValues();
//
//  CM value = (D1.diagonalToVector()).sort().getMatrix().transpose().vectorToDiagonal();
//  CM vector = this.A.createZero(D1.getRowSize(), D1.getColumnSize()).toComplex();
//  vector.setColumnVectors((D1.diagonalToVector()).sort().getIndices().transpose(), P1);
//
//  this.verifiedValue = solve(value, vector);
//
//  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solveAsSparse() {
//    if (this.A != null) {
//      solveAsSparse(this.A);
//    } else if (this.IA != null) {
//      solveAsSparse(this.IA);
//    } else {
//      throw new RuntimeException("Target matrix is not given."); //$NON-NLS-1$
//    }
//  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solve(final RM a) {
//    this.A = a;
//
//    EigenSolution<RS,RM,CS,CM> eigen = (this.A).eigenDecompose();
//
//    CM P1 = eigen.getVectors();
//    CM D1 = eigen.getValues();
//
//    CM value = (D1.diagonalToVector()).sort().getMatrix().transpose().vectorToDiagonal();
//    CM vector = this.A.createZero(D1.getRowSize(), D1.getColumnSize()).toComplex();
//    vector.setColumnVectors((D1.diagonalToVector()).sort().getIndices().transpose(), P1);
//
//    this.verifiedValue = solve(value, vector);
//  }

//  /**
//   * @see org.mklab.cga.eigen.EigenVerifier#getSolution()
//   */
//  public EigenVerifiedSolution<E> getSolution() {
//    return new EigenVerifiedSolution<>(this.verifiedValue);
//  }

  /**
   * {@inheritDoc}
   */
  public CIM getEigenValue() {
    return this.verifiedValue;
  }
  
  /**
   * {@inheritDoc}
   */
  public CIM getEigenVector() {
    throw new UnsupportedOperationException();
  }
  
  /**
   * {@inheritDoc}
   */
  public void solve() {
    EigenSolution<RS,RM,CS,CM> eigen = this.A.eigenDecompose();

    CM P1 = eigen.getVector();
    CM D1 = eigen.getValue();

    CM value = (D1.diagonalToVector()).sort().getMatrix().transpose().vectorToDiagonal();
    CM vector = this.A.createZero(D1.getRowSize(), D1.getColumnSize()).toComplex();
    vector.setColumnVectors((D1.diagonalToVector()).sort().getIndices().transpose(), P1);

    this.verifiedValue = solve(value, vector);
  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solveAsSparse(RM a) {
//    throw new UnsupportedOperationException();
//  }

  /**
   * 精度保証付き固有値を返します。
   * 
   * @param value 近似固有値
   * @param vector 近似固有ベクトル
   * 
   * @return 精度保証付き固有値
   */
  private CIM solve(CM value, CM vector) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    CM v_down = vector.multiply(value);

    manager.setRoundMode(RoundMode.ROUND_UP);

    CM v_up = vector.multiply(value);
    CM mV = v_down.add((v_up.subtract(v_down)).divide(2));
    CM rV = mV.subtract(v_down);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    CM e_down = mV.multiply(vector.transpose()).add(rV.multiply(-1).multiply((vector.transpose()).absElementWise())).subtract(this.A.toComplex());
    CM f_down = vector.multiply(vector.transpose()).subtract(this.A.createUnit().toComplex());

    manager.setRoundMode(RoundMode.ROUND_UP);

    CM e_up = mV.multiply(vector.transpose()).add(rV.multiply((vector.transpose()).absElementWise())).subtract(this.A.toComplex());
    CM f_up = vector.multiply(vector.transpose()).subtract(this.A.createUnit().toComplex());

    CM e_max = (e_down).absElementWise().maxElementWise((e_up).absElementWise());
    CM f_max = ((f_down).absElementWise()).maxElementWise((f_up).absElementWise());

    CS e_norm = ((e_max).norm(NormType.INFINITY));
    CS f_norm = ((f_max).norm(NormType.INFINITY));

    //NumericalMatrixOperator<?> approximateSolution = (NumericalMatrixOperator<?>)value.diagonalToVector();
    RM error = (value).absElementWise().getRealPart().multiply(f_norm.getRealPart()).add(this.A.createUnit().multiply(e_norm.getRealPart())).diagonalToVector();

    CIM ans2 =  this.IA.toComplex().createMidRad(value.diagonalToVector(), error);

    manager.setRoundMode(oldRoundMode);
    return ans2;
  }

  /**
   * {@inheritDoc}
   */
  public void solveAsSparse() {
    throw new UnsupportedOperationException();
  }

//  /**
//   * {@inheritDoc}
//   */
//  public VerifiedSolution<RIS, RIM, CIS, CIM, RS, RM, CS, CM> getSolution() {
//    // TODO stub automaticaly generated
//    return null;
//  }
}