Gerschgorin.java

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

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.RealNumericalScalar;


/**
 * ゲルシュゴリンの定理を実装したクラスです。
 * 
 * <p>このクラスを使えばゲルシュゴリンの定理による 固有値の存在範囲を知ることができます。
 * 
 * ゲルシュゴリンの定理は固有値の存在範囲を大まかに 示すだけなので精度保証に用いるにはふさわしくありません。
 * 
 * @author Hiroki
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 * @param <RS> スカラーの型
 * @param <RM> 行列の型
 */
public class Gerschgorin<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 x;

  /**
   * 任意の行列を設定します。
   * 
   * @param x 任意の行列。
   */
  public Gerschgorin(final RM x) {
    this.x = x;
  }

  /**
   * メイン。
   * 
   * @param args 任意の文字列。
   */
  public static void main(final String[] args) {
    DoubleMatrix A = new DoubleMatrix(new double[][] { {1, 0, 1}, {-1, -2, -0.5}, {0.5, -1, 3}});

    // Matrix A = Matrix.random(1000, 1000);
    Gerschgorin<DoubleNumber,DoubleMatrix,DoubleComplexNumber,DoubleComplexMatrix> g = new Gerschgorin<>(A);
    DoubleMatrix[] ans = g.solve();
    ans[0].print("R"); //$NON-NLS-1$
    ans[1].print("r"); //$NON-NLS-1$

  }

  /**
   * ゲルシュゴリンの定理により、固有値の 存在範囲を示します。
   * 
   * @return 第一要素は固有値の存在範囲の中心点、第二要素は中心点からの存在範囲の半径です。
   */
  @SuppressWarnings("unchecked")
  public RM[] solve() {
    int length = this.x.length();
    RS[] r_elements = this.x.getElement(1, 1).createArray(length);
    RM r = this.x.createZero(length, 1);
    RM R = this.x.createZero(length, 1);

    for (int i = 0; i < length; i++) {
      for (int j = 0; j < length; j++) {
        if (i != j) {
          r_elements[i] = r_elements[i].add(this.x.getElement(j + 1, i + 1).abs());
          R.setElement(i + 1, 1, this.x.getElement(i + 1, i + 1).unaryMinus());
          r.setElement(i + 1, 1, r_elements[i]);
        }

      }
    }

    RM[] ans = (RM[])new RealNumericalMatrix[2];
    ans[0] = R;
    ans[1] = r;
    return ans;
  }
}