MekekeMethod.java

package org.mklab.cga.eigen;

import org.mklab.cga.interval.matrix.DoubleIntervalMatrix;
import org.mklab.cga.round.FPURoundModeSelector;
import org.mklab.nfc.eig.DoubleEigenSolution;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/*
 * Created on 2004/12/07
 *
 * TODO To change the template for this generated file go to
 * Window - Preferences - Java - Code Style - Code Templates
 */

/**
 * 大石の方法で対角化可能行列の固有値の精度保証を行うクラスです。
 * 
 * <p><code>BMFMethod</code>クラスとは、誤差評価の式が異なっています。
 * 
 * @author hiroki
 */
public class MekekeMethod {

  /**
   * @param args コマンドライン引数
   */
  public static void main(String[] args) {
    RoundModeManager manager = RoundModeManager.getManager();
    manager.add(new FPURoundModeSelector());

    DoubleMatrix A = new DoubleMatrix(new double[][] { {1, 2, 1, 2}, {2, 2, -1, 1}, {1, -1, 1, 1}, {2, 1, 1, 1}});
    DoubleEigenSolution PD = A.eigenDecompose();

    DoubleMatrix eigenVector = PD.getVector().getRealPart();
    DoubleMatrix eigenValue = PD.getValue().getRealPart();

    DoubleMatrix D = eigenValue.diagonalToVector().sort().getMatrix().transpose().vectorToDiagonal();
    DoubleMatrix P = new DoubleMatrix(eigenValue.getRowSize(), eigenValue.getColumnSize());
    P.setColumnVectors(eigenValue.diagonalToVector().sort().getIndices().transpose(), eigenVector);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    DoubleMatrix v_down = P.multiply(D);

    manager.setRoundMode(RoundMode.ROUND_UP);

    DoubleMatrix v_up = P.multiply(D);
    DoubleMatrix mV = v_down.add((v_up.subtract(v_down)).multiply(0.5));
    DoubleMatrix rV = mV.subtract(v_down);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    DoubleMatrix e_down = mV.multiply(P.transpose()).add(rV.multiply(-1).multiply(P.transpose().absElementWise())).subtract(A);
    DoubleMatrix f_down = P.multiply(P.transpose()).subtract(DoubleMatrix.unit(A.getColumnSize()));

    manager.setRoundMode(RoundMode.ROUND_UP);

    DoubleMatrix e_up = mV.multiply(P.transpose()).add(rV.multiply(P.transpose().absElementWise())).subtract(A);
    DoubleMatrix f_up = P.multiply(P.transpose()).subtract(DoubleMatrix.unit(A.getColumnSize()));

    DoubleMatrix e_max = e_down.absElementWise().maxElementWise(e_up.absElementWise());
    DoubleMatrix f_max = f_down.absElementWise().maxElementWise(f_up.absElementWise());

    double e_norm = e_max.norm(NormType.INFINITY).doubleValue();
    double f_norm = f_max.norm(NormType.INFINITY).doubleValue();

    DoubleMatrix approximateSolution = D.diagonalToVector();
    DoubleMatrix error = D.absElementWise().multiply(f_norm).add(DoubleMatrix.unit(A.getColumnSize()).multiply(e_norm)).diagonalToVector();

    approximateSolution.print(Messages.getString("MekekeMethod.0")); //$NON-NLS-1$
    error.print(Messages.getString("MekekeMethod.1")); //$NON-NLS-1$
    DoubleIntervalMatrix L =  DoubleIntervalMatrix.createMiddleRadius(D.diagonalToVector(), error);
    L.printMidRad("L"); //$NON-NLS-1$
    L.printInfSup("L"); //$NON-NLS-1$

  }
}