Schord.java

/*
 * $Id: Schord.java,v 1.21 2008/07/16 15:40:00 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.IntMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.matrix.misc.GivensMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;


/**
 * 順序付きシュア分解を求めるクラスです。
 * 
 * <p>Schur decomposition with ordering
 * 
 * @author koga
 * @version $Revision: 1.21 $
 */
public class Schord {

  /**
   * <code>Ti</code>の固有値が指数<code>idx</code>の成分の昇順に並ぶようにするための直交行列<code>Q</code>を返します。
   * 
   * <p> <code>Ti</code>は上三角行列、<code>Qi</code>は直交行列です。
   * 
   * @param Qi Q行列
   * @param Ti T行列
   * @param idx_ 成分の順番
   * @return {Qo, To} 順序付きシュア分解 (schur decomposition)
   */
  public static List<DoubleComplexMatrix> schord(DoubleMatrix Qi, DoubleMatrix Ti, IntMatrix idx_) {
    if (Ti.getRowSize() != Ti.getColumnSize()) {
      throw new IllegalArgumentException(Messages.getString("Schord.0")); //$NON-NLS-1$
    }
    int size = Ti.getRowSize();

    IntMatrix idx = idx_;
    DoubleComplexMatrix To = new DoubleComplexMatrix(Ti);
    DoubleComplexMatrix Qo = new DoubleComplexMatrix(Qi);

    for (int i = 1; i <= size - 1; i++) {
      boolean flag = false;
      int k = i;
      for (int j = i + 1; j <= size; j++) {
        if (idx.getIntElement(j) < idx.getIntElement(k)) {
          k = j;
          flag = true;
        }
      }

      if (flag) {
        for (int p = k; p >= i + 1; p = p - 1) {
          //Matrix G = Givens.givens((DoubleNumber)To.getElement(p - 1, p - 1).subtract(To.getElement(p, p)), To.getElement(p - 1, p)).flipUpDown();
          DoubleComplexMatrix G = GivensMatrix.create(To.getElement(p - 1, p - 1).subtract(To.getElement(p, p)), To.getElement(p - 1, p)).flipUpDown();

          // To.setSubMatrix(1, To.getRowSize(), p - 1, p, To.getSubMatrix(1, To.getRowSize(), p - 1, p).multiply(G));
          To.setColumnVectors(p - 1, p, To.getColumnVectors(p - 1, p).multiply(G));

          // To.setSubMatrix(p - 1, p, 1, To.getColumnSize(), G.conjugateTranspose().multiply(To.getSubMatrix(p - 1, p, 1, To.getColumnSize())));
          To.setRowVectors(p - 1, p, G.conjugateTranspose().multiply(To.getRowVectors(p - 1, p)));

          // Qo.setSubMatrix(1, Qo.getRowSize(), p - 1, p, Qo.getSubMatrix(1, Qo.getRowSize(), p - 1, p).multiply(G));
          Qo.setColumnVectors(p - 1, p, Qo.getColumnVectors(p - 1, p).multiply(G));

          int ix = idx.getIntElement(p - 1);
          idx.setElement(1, p - 1, idx.getIntElement(p));
          idx.setElement(1, p, ix);
        }
      }
    }

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

  }

  /**
   * <code>Ti</code>の固有値が指数<code>idx</code>の成分の昇順に並ぶようにするための直交行列<code>Q</code>を返します。
   * 
   * <p> <code>Ti</code>は上三角行列、<code>Qi</code>は直交行列です。
   * 
   * @param <RS> スカラーの型
   * @param <RM> 行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param Qi Q行列
   * @param Ti T行列
   * @param idx_ 成分の順番
   * @return {Qo, To} 順序付きシュア分解 (schur decomposition)
   */
  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> schord(RM Qi, RM Ti, IntMatrix idx_) {
    if (Ti.getRowSize() != Ti.getColumnSize()) {
      throw new IllegalArgumentException(Messages.getString("Schord.0")); //$NON-NLS-1$
    }
    int size = Ti.getRowSize();

    IntMatrix idx = idx_;
    CM To = Ti.toComplex();
    CM Qo =Qi.toComplex();

    for (int i = 1; i <= size - 1; i++) {
      boolean flag = false;
      int k = i;
      for (int j = i + 1; j <= size; j++) {
        if (idx.getIntElement(j) < idx.getIntElement(k)) {
          k = j;
          flag = true;
        }
      }

      if (flag) {
        for (int p = k; p >= i + 1; p = p - 1) {
          //Matrix G = Givens.givens((DoubleNumber)To.getElement(p - 1, p - 1).subtract(To.getElement(p, p)), To.getElement(p - 1, p)).flipUpDown();
          CM G = GivensMatrix.create(To.getElement(p - 1, p - 1).subtract(To.getElement(p, p)), To.getElement(p - 1, p)).flipUpDown();

          // To.setSubMatrix(1, To.getRowSize(), p - 1, p, To.getSubMatrix(1, To.getRowSize(), p - 1, p).multiply(G));
          To.setColumnVectors(p - 1, p, To.getColumnVectors(p - 1, p).multiply(G));

          // To.setSubMatrix(p - 1, p, 1, To.getColumnSize(), G.conjugateTranspose().multiply(To.getSubMatrix(p - 1, p, 1, To.getColumnSize())));
          To.setRowVectors(p - 1, p, G.conjugateTranspose().multiply(To.getRowVectors(p - 1, p)));

          // Qo.setSubMatrix(1, Qo.getRowSize(), p - 1, p, Qo.getSubMatrix(1, Qo.getRowSize(), p - 1, p).multiply(G));
          Qo.setColumnVectors(p - 1, p, Qo.getColumnVectors(p - 1, p).multiply(G));

          int ix = idx.getIntElement(p - 1);
          idx.setElement(1, p - 1, idx.getIntElement(p));
          idx.setElement(1, p, ix);
        }
      }
    }

    List<CM> ans = new ArrayList<>();
    ans.add(Qo);
    ans.add(To);
    return ans;
  }

}