Rsf2csf.java

/*
 * $Id: Rsf2csf.java,v 1.27 2008/03/26 14:31:23 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.NormType;
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>Real Schur to complex Schur form conversion
 * 
 * @author koga
 * @version $Revision: 1.27 $
 */
public class Rsf2csf {
  /**
   * 実擬似上シュア形の行列を複素上シュア形に変換します。 
   * 
   * @param Ui 実シュア形式
   * @param Ti 実シュア形式
   * @return {Uo, To} 複素シュア形式 (complex schur form)
   */
  public static List<DoubleComplexMatrix> rsf2csf(DoubleMatrix Ui, DoubleMatrix Ti) {
    final int n = Ti.length();
    final DoubleComplexMatrix To = new DoubleComplexMatrix(Ti);
    final DoubleComplexMatrix Uo = new DoubleComplexMatrix(Ui);
    final DoubleNumber tolerance = To.frobNorm().multiply(To.frobNorm().getMachineEpsilon()).getRealPart();

    for (int i = n; i > 1; i--) {
      final DoubleComplexNumber ToElement1 = To.getElement(i,  i -1);
      final DoubleComplexNumber ToElement2 = To.getElement(i, i);

      if (ToElement1.abs().isGreaterThan(tolerance)) {
        final DoubleComplexMatrix mu = To.getSubMatrix(i - 1, i, i - 1, i).eigenValue().subtractElementWise(ToElement2);
        final DoubleComplexNumber[] elements = new DoubleComplexNumber[] {mu.getElement(1, 1),ToElement1};
        final DoubleComplexNumber norm = elements[0].createGrid(elements).norm(NormType.TWO);
        final DoubleComplexNumber c = mu.getElement(1, 1).divide(norm);
        final DoubleComplexNumber s = ToElement1.divide(norm);
        final DoubleComplexNumber[][] elements2 = new DoubleComplexNumber[][] { {c.conjugate(), s}, {s.unaryMinus(), c}};
        final DoubleComplexMatrix Q = elements2[0][0].createGrid(2, 2, elements2);
        To.setSubMatrix(i - 1, i, i - 1, n, Q.multiply(To.getSubMatrix(i - 1, i, i - 1, n)));
        To.setSubMatrix(1, i, i - 1, i, To.getSubMatrix(1, i, i - 1, i).multiply(Q.conjugateTranspose()));
        Uo.setSubMatrix(1, n, i - 1, i, Uo.getSubMatrix(1, n, i - 1, i).multiply(Q.conjugateTranspose()));
      }

      To.setElement(i, i - 1, 0);
    }

    return new ArrayList<>(Arrays.asList(new DoubleComplexMatrix[] {Uo, To}));
  }

  /**
   * 実擬似上シュア形の行列を複素上シュア形に変換します。 
   * 
   * @param <RS> スカラーの型
   * @param <RM> 行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param Ui 実シュア形式
   * @param Ti 実シュア形式
   * @return {Uo, To} 複素シュア形式 (complex schur form)
   */
  public static <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>> List<CM> rsf2csf(RM Ui, RM Ti) {
    final int n = Ti.length();
    final CM To = Ti.toComplex();
    final CM Uo =Ui.toComplex();
    final RS tolerance = To.frobNorm().multiply(To.frobNorm().getMachineEpsilon()).getRealPart();

    for (int i = n; i > 1; i--) {
      final CS ToElement1 = To.getElement(i,  i -1);
      final CS ToElement2 = To.getElement(i, i);

      if (ToElement1.abs().getRealPart().isGreaterThan(tolerance)) {
        final CM mu = To.getSubMatrix(i - 1, i, i - 1, i).eigenValue().subtractElementWise(ToElement2);
        final CS[] elements = ToElement1.createArray(2);
        elements[0] = mu.getElement(1, 1);
        elements[1] = ToElement1;
        
        final CS norm = elements[0].createGrid(elements).norm(NormType.TWO);
        final CS c = mu.getElement(1, 1).divide(norm);
        final CS s = ToElement1.divide(norm);
        final CS[][] elements2 = c.createArray(2, 2);
        elements2[0][0] = c.conjugate();
        elements2[0][1] = s;
        elements2[1][0] = s.unaryMinus();
        elements2[1][1] = c;
        
        final CM Q = elements2[0][0].createGrid(2, 2, elements2);
        To.setSubMatrix(i - 1, i, i - 1, n, Q.multiply(To.getSubMatrix(i - 1, i, i - 1, n)));
        To.setSubMatrix(1, i, i - 1, i, To.getSubMatrix(1, i, i - 1, i).multiply(Q.conjugateTranspose()));
        Uo.setSubMatrix(1, n, i - 1, i, Uo.getSubMatrix(1, n, i - 1, i).multiply(Q.conjugateTranspose()));
      }

      To.setElement(i, i - 1, 0);
    }
    
    List<CM> ans = new ArrayList<>();  
    ans.add(Uo);
   ans.add(To);
   return ans;
  }
}