ComplexBFMethod.java

/*
 * Created on 2004/11/11
 * 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.20 $.2004/11/11
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class ComplexBFMethod<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 ComplexEigenVerifier<RIS,RIM,CIS,CIM,RS,RM,CS,CM> {

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

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

//  /**
//   * コンストラクタ
//   * 
//   */
  
//  /**
//   * Creates {@link ComplexBFMethod}.
//   */
//  public ComplexBFMethod() {
//    // nothing to do
//  }

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

//  /**
//   * 新しく生成された<code>BFMethod</code>オブジェクトを初期化します。
//   * 
//   * @param IA 対象となる区間行列
//   */
//  public ComplexBFMethod(CIM IA) {
//    this.IA = IA;
//    this.A = IA.getMiddle();
//  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solve() {
//    EigenSolution<RS,RM,CS,CM> eigen = this.A.eigenDecompose();
//    CM value = eigen.getValues();
//    CM vector = eigen.getVectors();
//    this.verifiedValue = solve(value, vector);
//  }

//  /**
//   * @see org.mklab.cga.eigen.EigenVerifier#solveAsSparse()
//   */
//  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 CM a) {
//    this.A = a;
//    EigenSolution<RS,RM,CS,CM> eigen = (this.A).eigenDecompose();
//    CM value = eigen.getValues();
//    CM vector = eigen.getVectors();
//    this.verifiedValue = solve(value, vector);
//  }

  /**
   * {@inheritDoc}
   */
  public void solve(final CIM a) {
    this.IA = a;
    this.A = a.getMiddle();
    EigenSolution<RS,RM,CS,CM> eigen = this.A.eigenDecompose();
    CM value = eigen.getValue();
    CM vector = eigen.getVector();
    this.verifiedValue = solve(value, vector);
  }

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

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

  /**
   * @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 v_m = v_down.add((v_up.subtract(v_down)).divide(2));
    CM v_r = v_m.subtract(v_down);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    CM e_down = v_m.multiply(vector.inverse()).add(v_r.multiply(-1).multiply((vector.inverse()).absElementWise())).subtract(this.A);
    CM f_down = vector.multiply(vector.inverse()).subtract(this.A.createUnit());

    manager.setRoundMode(RoundMode.ROUND_UP);

    CM e_up = v_m.multiply(vector.inverse()).add(v_r.multiply((vector.inverse()).absElementWise())).subtract(this.A);
    CM f_up = vector.multiply(vector.inverse()).subtract(this.A.createUnit());
    CM e_max = ((e_down).absElementWise()).maxElementWise((e_up).absElementWise());
    CM f_max = (f_down).absElementWise().maxElementWise((f_up).absElementWise());

    // TODO エラーを解消する
    CS e_norm = ((e_max).norm(NormType.TWO));
    // TODO エラーを解消する
    CS A_norm = ((this.A).norm(NormType.TWO));

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    CS f_norm = ((f_max).norm(NormType.TWO));

    manager.setRoundMode(RoundMode.ROUND_UP);

    RM error2 = this.A.getRealPart().createUnit(value.getRowSize(), value.getColumnSize()).multiply(e_norm.getRealPart().add(A_norm.getRealPart().multiply(f_norm.getRealPart()))).diagonalToVector();

    manager.setRoundMode(RoundMode.ROUND_NEAR);

    CIM L = this.IA.createMidRad(value.diagonalToVector(), error2);

    manager.setRoundMode(oldRoundMode);
    return L;
  }

  /**
   * {@inheritDoc}
   */
  public void solveAsSparse(CIM a) {
    throw new UnsupportedOperationException();
  }
  
  /**
   * {@inheritDoc}
   */
  public CIM getEigenValue() {
    return this.verifiedValue;
  }

  /**
   * {@inheritDoc}
   */
  public CIM getEigenVector() {
    throw new UnsupportedOperationException();
  }
}