Funm.java

/*
 * $Id: Funm.java,v 1.33 2008/07/16 15:40:00 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import java.util.List;

import org.mklab.nfc.eig.SchurDecomposition;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleNumber;


/**
 * 行列関数の値を求めるクラスです。
 * 
 * <p>Matrix functions
 * 
 * @author koga
 * @version $Revision: 1.33 $
 */
public class Funm {

  /**
   * 複素数関数を表すクラス
   * 
   * @author koga
   */
  public static abstract class ComplexFunction {

    /**
     * 関数を評価します。
     * 
     * @param c 関数の引数
     * @return 評価の結果
     */
    public abstract DoubleComplexNumber eval(DoubleComplexNumber c);
  }

  /**
   * 行列関数を評価します。
   * 
   * @param A 対象となる行列
   * @param fun 関数
   * @return 行列関数の値 (evaluated result)
   */
  public static DoubleComplexMatrix funm(DoubleMatrix A, ComplexFunction fun) {
    int n = A.getRowSize();
    SchurDecomposition<DoubleNumber,DoubleMatrix> temp = A.schurDecompose();
    DoubleMatrix Qr = temp.getU();
    DoubleMatrix Tr = temp.getT();

    List<DoubleComplexMatrix> tmp = Rsf2csf.rsf2csf(Qr, Tr);
    DoubleComplexMatrix Q = tmp.get(0);
    DoubleComplexMatrix T = tmp.get(1);

    DoubleComplexNumber unit = fun.eval(T.getElement(1, 1));
    DoubleComplexMatrix F = unit.createZeroGrid(n, n);
    for (int i = 1; i <= n; i++) {
      F.setElement(i, i, fun.eval(T.getElement(i, i)));
    }

    DoubleComplexNumber tNorm = T.norm(NormType.ONE);
    DoubleComplexNumber tolerance = tNorm.multiply(tNorm.getMachineEpsilon());

    for (int k = 1; k <= n - 1; k++) {
      for (int i = 1; i <= n - k; i++) {
        int j = i + k;
        DoubleComplexNumber s = T.getElement(i, j).multiply(F.getElement(j, j).subtract(F.getElement(i, i)));
        if (k > 1) {
          DoubleComplexMatrix tp = T.getSubMatrix(i, i, i + 1, j - 1).multiply(F.getSubMatrix(i + 1, j - 1, j, j)).subtract(F.getSubMatrix(i, i, i + 1, j - 1).multiply(T.getSubMatrix(i + 1, j - 1, j, j)));
          s = s.add(tp.getElement(1));
        }
        DoubleComplexNumber d = T.getElement(j, j).subtract(T.getElement(i, i));

        if (d.abs().isLessThanOrEquals(tolerance)) {
          System.err.println(Messages.getString("Funm.0")); //$NON-NLS-1$
          d = d.createUnit().multiply(tolerance);
        }
        F.setElement(i, j, s.divide(d));
      }
    }

    F = Q.multiply(F).multiply(Q.conjugateTranspose());
    return F;
  }

  /**
   * 行列関数を評価します。
   * 
   * @param A 対象となる行列
   * @param fun 関数
   * @param tolerance 許容誤差
   * @return 行列関数の値 (evaluated result)
   */
  public static DoubleComplexMatrix funm(DoubleMatrix A, ComplexFunction fun, double tolerance) {
    int n = A.getRowSize();
    SchurDecomposition<DoubleNumber,DoubleMatrix> temp = A.schurDecompose();
    DoubleMatrix Qr = temp.getU();
    DoubleMatrix Tr = temp.getT();

    List<DoubleComplexMatrix> tmp = Rsf2csf.rsf2csf(Qr, Tr);
    DoubleComplexMatrix Q = tmp.get(0);
    DoubleComplexMatrix T = tmp.get(1);

    DoubleComplexNumber unit = fun.eval(T.getElement(1, 1));
    DoubleComplexMatrix F = unit.createZeroGrid(n, n);

    for (int i = 1; i <= n; i++) {
      F.setElement(i, i, fun.eval(T.getElement(i, i)));
    }

    for (int k = 1; k <= n - 1; k++) {
      for (int i = 1; i <= n - k; i++) {
        int j = i + k;
        DoubleComplexNumber s = T.getElement(i, j).multiply(F.getElement(j, j).subtract(F.getElement(i, i)));
        if (k > 1) {
          DoubleComplexMatrix tp = T.getSubMatrix(i, i, i + 1, j - 1).multiply(F.getSubMatrix(i + 1, j - 1, j, j)).subtract(F.getSubMatrix(i, i, i + 1, j - 1).multiply(T.getSubMatrix(i + 1, j - 1, j, j)));
          s = s.add(tp.getElement(1));
        }
        DoubleComplexNumber d = T.getElement(j, j).subtract(T.getElement(i, i));

        if (d.abs().isLessThanOrEquals(tolerance)) {
          System.err.println(Messages.getString("Funm.1")); //$NON-NLS-1$
          d = d.createUnit().multiply(tolerance);
        }
        F.setElement(i, j, s.divide(d));
      }
    }

    F = Q.multiply(F).multiply(Q.conjugateTranspose());
    return F;
  }
}