AlgebraWithEFFloat.java

/**
 * Created on 2014/01/28
 * Copyright (C) 2014 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.sdpj.algebra;

import org.mklab.mpf.ap.EFFloat;
import org.mklab.mpf.ap.EFFloatComplexMatrix;
import org.mklab.mpf.ap.EFFloatComplexNumber;
import org.mklab.mpf.ap.EFFloatMatrix;
import org.mklab.mpf.ap.MPFloat;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.BlockVector;
import org.mklab.sdpj.datastructure.DenseDiagonal;
import org.mklab.sdpj.datastructure.DenseMatrix;
import org.mklab.sdpj.datastructure.SparseDenseDiagonal;
import org.mklab.sdpj.datastructure.SparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.gpack.blaswrap.BLAS;
import org.mklab.sdpj.gpack.lapackwrap.Clapack;
import org.mklab.sdpj.tool.Tools;


/**
 * EFFloatを使って行列の演算を行います.
 * 
 * 
 * @author Hashimoto
 * @version $Revision$, 2014/01/28
 */
public class AlgebraWithEFFloat {

  /**
   * 最小固有値の計算を行います aMat is rewritten. aMat must be symmetric. eigenVec is the space of eigen values and needs memory of length aMat.nRow workVec is temporary space and needs 3*aMat.nRow-1 length
   * memory.
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param eigenValues 固有値
   * @param workVec 作業領域
   * @return 最小固有値
   */
  private 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>> RS getMinEigenValue(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> eigenValues, Vector<RS,RM,CS,CM> workVec) {
    final RS unit = a.getElementUnit();

    int n = a.getRowSize();

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        int lwork = 3 * n - 1;
        // "N" means that we need not eigen vectors
        // "L" means that we refer only lower triangular.
        int info = Clapack.dsyev("NonVectors", "Lower", n, a.denseElements, n, eigenValues.getElements(), workVec.getElements(), lwork); //$NON-NLS-1$ //$NON-NLS-2$
        if (info != 0) {
          return unit.createZero();
        }
        return eigenValues.getElements()[0];
      case DIAGONAL:
        RS minimumEigenValue = a.diagonalElements[0];
        eigenValues.getElements()[0] = a.diagonalElements[0];
        for (int i = 1; i < n; ++i) {
          eigenValues.getElements()[i] = a.diagonalElements[i];
          if (a.diagonalElements[i].isLessThan(minimumEigenValue)) {
            minimumEigenValue = a.diagonalElements[i];
          }
        }
        return minimumEigenValue;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * 最小固有値の計算を行います。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param eigenVec 固有ベクトル
   * @param workVec 作業領域
   * @return 最小固有値
   */
  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>> RS getMinEigenValue(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> eigenVec, BlockVector<RS,RM,CS,CM> workVec) {
    int n = a.getBlock(0).getRowSize();
    if (eigenVec.getBlock(0).getDimension() != n) {
      throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
    }
    if (workVec.getBlock(0).getDimension() != 3 * n - 1) {
      throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
    }
    RS min_eigen = getMinEigenValue(a.getBlock(0), eigenVec.getBlock(0), workVec.getBlock(0));

    for (int l = 1; l < a.getBlockSize(); ++l) {
      n = a.getBlock(l).getRowSize();
      if (eigenVec.getBlock(l).getDimension() != n) {
        throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
      }
      if (workVec.getBlock(l).getDimension() != 3 * n - 1) {
        throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
      }
      RS tmp_eigen = getMinEigenValue(a.getBlock(l), eigenVec.getBlock(l), workVec.getBlock(l));
      if (tmp_eigen.isLessThan(min_eigen)) {
        min_eigen = tmp_eigen;
      }
    }
    return min_eigen;
  }

  /**
   * 内積の計算を行います。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @return aVecとbVecの内積
   */
  private 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>> RS getInnerProduct(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
    int size = a.getDimension();
    if (size != b.getDimension()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }
    final RS unit = a.getElementUnit();

    RS ans = unit.create(0);

    if (ans instanceof MPFloat) {
      EFFloat result = new EFFloat(0);
      EFFloat[] effloatVectorElementsA = new EFFloat[size];
      EFFloat[] effloatVectorElementsB = new EFFloat[size];
      for (int j = 0; j < effloatVectorElementsA.length; j++) {
        effloatVectorElementsA[j] = new EFFloat((MPFloat)a.getElement(j));
        effloatVectorElementsB[j] = new EFFloat((MPFloat)b.getElement(j));
      }
      //      result = getInnerProduct(new Vector<EFFloat>(effloatVectorElementsA), new Vector<EFFloat>(effloatVectorElementsB));
      result = BLAS.<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> ddot(size, effloatVectorElementsA, 1, effloatVectorElementsB, 1);
      return (RS)result.toMPFloatWithDefaultPrecision();
    }
    return BLAS.<RS,RM,CS,CM> ddot(size, a.getElements(), 1, b.getElements(), 1);
  }

  /**
   * 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @return 最小固有値
   */
  private 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>> RS getInnerProduct(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
    // 内積＝Al.getInnerProduct(vec1,vec2)；に書き換えたほうがいいかもね！
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();

    RS ret = unit.create(0);
    if (ret instanceof MPFloat) {
      EFFloat result = new EFFloat(0);
      for (int i = 0; i < a.getBlockSize(); ++i) {
        EFFloat[] effloatBlockVectorElementsA = new EFFloat[a.getBlock(i).getDimension()];
        for (int j = 0; j < effloatBlockVectorElementsA.length; j++) {
          effloatBlockVectorElementsA[j] = new EFFloat((MPFloat)a.getBlock(i).getElement(j));
        }
        EFFloat[] effloatBlockVectorElementsB = new EFFloat[b.getBlock(i).getDimension()];
        for (int j = 0; j < effloatBlockVectorElementsB.length; j++) {
          effloatBlockVectorElementsB[j] = new EFFloat((MPFloat)b.getBlock(i).getElement(j));
        }

        result = result.add(getInnerProduct(new Vector<>(effloatBlockVectorElementsA), new Vector<>(effloatBlockVectorElementsB)));
      }
      return (RS)result.toMPFloat();
    }

    for (int i = 0; i < a.getBlockSize(); ++i) {
      ret = ret.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
    }
    return ret;
  }

  /**
   * 内積の計算を行います。 AI.getInnerProduct(ret,aVec,bVec)をAI.ret=getInnerProduct(aVec,bVec)に変更。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return aMat bMatの内積
   */
  private 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>> RS getInnerProduct(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        int length = a.getRowSize() * a.getColumnSize();
        return BLAS.<RS,RM,CS,CM> ddot(length, a.denseElements, 1, b.denseElements, 1);
      case DIAGONAL:
        return BLAS.ddot(a.getColumnSize(), a.diagonalElements, 1, b.diagonalElements, 1);
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 内積
   */
  private 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>> RS getInnerProduct(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final RS unit = b.getElementUnit();
    int index = 0;
    int counter = 0;
    RS ans;

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        // Attension: in SPARSE case, only half elements
        // are stored. And bMat must be DENSE case.
        ans = unit.create(0);
        int amari = a.nonZeroCount % 4;
        int shou = a.nonZeroCount / 4;
        for (index = 0; index < amari; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index);
          if (i == j) {
            ans = ans.add(value.multiply(b.denseElements[i + b.getRowSize() * j]));
          } else {
            ans = ans.add(value.multiply((b.denseElements[i + b.getRowSize() * j].add(b.denseElements[j + b.getRowSize() * i]))));
          }
        }

        index = amari;
        counter = 0;
        for (counter = 0; counter < shou; ++counter) {
          int i1 = a.rowIndex[index];
          int j1 = a.columnIndex[index];
          RS value1 = a.getSparseElement(index);
          RS ans1;
          if (i1 == j1) {
            ans1 = value1.multiply(b.denseElements[i1 + b.getRowSize() * j1]);
          } else {
            ans1 = value1.multiply((b.denseElements[i1 + b.getRowSize() * j1].add(b.denseElements[j1 + b.getRowSize() * i1])));
          }

          int i2 = a.rowIndex[index + 1];
          int j2 = a.columnIndex[index + 1];
          RS value2 = a.getSparseElement(index + 1);
          RS ans2;
          if (i2 == j2) {
            ans2 = value2.multiply(b.denseElements[i2 + b.getRowSize() * j2]);
          } else {
            ans2 = value2.multiply((b.denseElements[i2 + b.getRowSize() * j2].add(b.denseElements[j2 + b.getRowSize() * i2])));

          }

          int i3 = a.rowIndex[index + 2];
          int j3 = a.columnIndex[index + 2];
          RS value3 = a.getSparseElement(index + 2);
          RS ans3;
          if (i3 == j3) {
            ans3 = value3.multiply(b.denseElements[i3 + b.getRowSize() * j3]);
          } else {
            ans3 = value3.multiply((b.denseElements[i3 + b.getRowSize() * j3].add(b.denseElements[j3 + b.getRowSize() * i3])));
          }

          int i4 = a.rowIndex[index + 3];
          int j4 = a.columnIndex[index + 3];
          RS value4 = a.getSparseElement(index + 3);
          RS ans4;
          if (i4 == j4) {
            ans4 = value4.multiply(b.denseElements[i4 + b.getRowSize() * j4]);
          } else {
            ans4 = value4.multiply((b.denseElements[i4 + b.getRowSize() * j4].add(b.denseElements[j4 + b.getRowSize() * i4])));
          }

          index += 4;
          ans = ans.add(ans1.add(ans2).add(ans3).add(ans4));
          //20090120修正したけど　不具合でたらここもチェックで。
        }
        return ans;
      case DENSE:
        int length = a.getRowSize() * a.getColumnSize();
        return BLAS.ddot(length, a.denseElements, 1, b.denseElements, 1);
      case DIAGONAL:
        return BLAS.ddot(a.getColumnSize(), a.diagonalElements, 1, b.diagonalElements, 1);
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * 内積の計算を行います Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 内積
   */
  private 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>> RS getInnerProduct(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    RS ans = unit.create(0);
    if (ans instanceof MPFloat) {
      EFFloat result = new EFFloat(0);
      for (int i = 0; i < a.getBlockSize(); ++i) {
      	EFFloat[] effloatElementsA, effloatElementsB;
      	if (a.getBlock(i).isDense()){
      		effloatElementsA = new EFFloat[a.getBlock(i).denseElements.length];
      		for (int j = 0; j < effloatElementsA.length; j++) {
      			effloatElementsA[j] = new EFFloat((MPFloat)a.getBlock(i).denseElements[j]);
      		}
      		effloatElementsB = new EFFloat[b.getBlock(i).denseElements.length];
      		for (int j = 0; j < effloatElementsB.length; j++) {
      			effloatElementsB[j] = new EFFloat((MPFloat)b.getBlock(i).denseElements[j]);
      		}
      	} else {
      			effloatElementsA = new EFFloat[a.getBlock(i).diagonalElements.length];
      		for (int j = 0; j < effloatElementsA.length; j++) {
      			effloatElementsA[j] = new EFFloat((MPFloat)a.getBlock(i).diagonalElements[j]);
      		}
      			effloatElementsB = new EFFloat[b.getBlock(i).diagonalElements.length];
      		for (int j = 0; j < effloatElementsB.length; j++) {
      			effloatElementsB[j] = new EFFloat((MPFloat)b.getBlock(i).diagonalElements[j]);
      		}
      	}
        result = result.add(Algebra.<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> run(new Vector<>(effloatElementsA), '.', new Vector<>(effloatElementsB)));
      }
      return (RS)result.toMPFloatWithDefaultPrecision();
    }
    for (int i = 0; i < a.getBlockSize(); ++i) {
      ans = ans.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
    }
    return ans;
  }

  /**
   * 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 内積
   */
  private 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>> RS getInnerProduct(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    RS ans = unit.create(0);
    if (ans instanceof MPFloat) {
      EFFloat eunit = new EFFloat(1);
      
      BlockSparseMatrix<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> effloatSparseElementsA = new BlockSparseMatrix<>(a.getBlockSize(), a.getBlockStruct());
      BlockDenseMatrix<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> effloatDenseElementsB = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), eunit);

      for (int i = 0; i < a.getBlockSize(); i++) {
        SparseMatrix<RS,RM,CS,CM> sparseMatrix = a.getBlock(i);
      
        SparseMatrix<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> effloatSparseMatrix = new SparseMatrix<>(sparseMatrix.getRowSize(), sparseMatrix.getColumnSize(), sparseMatrix.getSparseOrDenseOrDiagonal(),
            sparseMatrix.getNonZeroNumber(), eunit);

        effloatSparseMatrix.nonZeroCount = sparseMatrix.nonZeroCount;
        effloatSparseMatrix.nonZeroEffect = sparseMatrix.nonZeroEffect;
        effloatSparseMatrix.rowIndex = sparseMatrix.rowIndex;
        effloatSparseMatrix.columnIndex = sparseMatrix.columnIndex;

        if (sparseMatrix.getSparseOrDenseOrDiagonal() == SparseDenseDiagonal.SPARSE) {
          effloatSparseMatrix.sparseElements = new EFFloat[sparseMatrix.sparseElements.length];
          for (int j = 0; j < sparseMatrix.sparseElements.length; j++) {
            MPFloat value = (MPFloat)sparseMatrix.sparseElements[j];
            effloatSparseMatrix.sparseElements[j] = new EFFloat(value);
          }
        }
        if (sparseMatrix.getSparseOrDenseOrDiagonal() == SparseDenseDiagonal.DENSE) {
          effloatSparseMatrix.denseElements = new EFFloat[sparseMatrix.denseElements.length];
          for (int j = 0; j < sparseMatrix.denseElements.length; j++) {
            MPFloat value = (MPFloat)sparseMatrix.denseElements[j];
            effloatSparseMatrix.denseElements[j] = new EFFloat(value);
          }
        }
        if (sparseMatrix.getSparseOrDenseOrDiagonal() == SparseDenseDiagonal.DIAGONAL) {
          effloatSparseMatrix.diagonalElements = new EFFloat[sparseMatrix.diagonalElements.length];
          for (int j = 0; j < sparseMatrix.diagonalElements.length; j++) {
            MPFloat value = (MPFloat)sparseMatrix.diagonalElements[j];
            effloatSparseMatrix.diagonalElements[j] = new EFFloat(value);
          }
        }

        effloatSparseElementsA.setBlock(i, effloatSparseMatrix);
      }

      for (int i = 0; i < b.getBlockSize(); i++) {
        DenseMatrix<RS,RM,CS,CM> denseMatrix = b.getBlock(i);
        DenseMatrix<EFFloat,EFFloatMatrix,EFFloatComplexNumber,EFFloatComplexMatrix> effloatDenseMatrix = new DenseMatrix<>(denseMatrix.getRowSize(), denseMatrix.getColumnSize(), denseMatrix.getDenseOrDiagonal(), eunit);

        if (denseMatrix.getDenseOrDiagonal() == DenseDiagonal.DENSE) {
          effloatDenseMatrix.denseElements = new EFFloat[denseMatrix.denseElements.length];
          for (int j = 0; j < denseMatrix.denseElements.length; j++) {
            MPFloat value = (MPFloat)denseMatrix.denseElements[j];
            effloatDenseMatrix.denseElements[j] = new EFFloat(value);
          }
        }
        if (denseMatrix.getDenseOrDiagonal() == DenseDiagonal.DIAGONAL) {
          effloatDenseMatrix.diagonalElements = new EFFloat[denseMatrix.diagonalElements.length];
          for (int j = 0; j < denseMatrix.diagonalElements.length; j++) {
            MPFloat value = (MPFloat)denseMatrix.diagonalElements[j];
            effloatDenseMatrix.diagonalElements[j] = new EFFloat(value);
          }
        }

        effloatDenseElementsB.setBlock(i, effloatDenseMatrix);
      }

      return (RS)Algebra.run(effloatSparseElementsA, '.', effloatDenseElementsB).toMPFloatWithDefaultPrecision();
    }

    for (int i = 0; i < a.getBlockSize(); ++i) {
      ans = ans.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
    }
    return ans;
  }

  /**
   * getCholesky(retMat,aMat)をret=getCholesky(aMat)に変更
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a 元になる行列
   * @return ret 結果retMatと同じ
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> getCholesky(DenseMatrix<RS,RM,CS,CM> a) {
    final RS unit = a.getElementUnit();

    DenseMatrix<RS,RM,CS,CM> ans;
    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        ans = a.createClone();
        int info = Clapack.dpotrf("Lower", ans.getRowSize(), ans.denseElements, ans.getRowSize()); //$NON-NLS-1$

        if (info != 0) {
          throw new RuntimeException("cannot cholesky decomposition\n" + "Could you try with smaller gammaStar?"); //$NON-NLS-1$ //$NON-NLS-2$
        }

        // Make matrix as lower triangular matrix
        for (int j = 0; j < ans.getColumnSize(); ++j) {
          int shou = j / 4;
          int amari = j % 4;
          for (int i = 0; i < amari; ++i) {
            ans.denseElements[i + ans.getColumnSize() * j] = unit.create(0);
          }
          int i;
          int count;
          for (i = amari, count = 0; count < shou; ++count, i += 4) {
            ans.denseElements[i + ans.getColumnSize() * j] = unit.create(0);
            ans.denseElements[i + 1 + ans.getColumnSize() * j] = unit.create(0);
            ans.denseElements[i + 2 + ans.getColumnSize() * j] = unit.create(0);
            ans.denseElements[i + 3 + ans.getColumnSize() * j] = unit.create(0);
          }
        }
        return ans;
      case DIAGONAL:
        ans = new DenseMatrix<>(a.getRowSize(), a.getColumnSize(), a.getDenseOrDiagonal(), unit);

        int shou = a.getColumnSize() / 4;
        int amari = a.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.diagonalElements[j].sqrt();
        }

        int j;
        int count;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.diagonalElements[j].sqrt();
          ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].sqrt();
          ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].sqrt();
          ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].sqrt();
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @return retMatのコピー
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> getInvLowTriangularMatrix(DenseMatrix<RS,RM,CS,CM> a) {
    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), a.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        ans.setIdentity(unit.create(1));

        int nRow = a.getRowSize();
        int nCol = a.getColumnSize();
        int retnRow = ans.getRowSize();
        BLAS.dtrsm("Left", "Lower", "NoTraspose", "NonUnitDiagonal", nRow, nCol, unit.create(1), a.denseElements, nRow, ans.denseElements, retnRow); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
        return ans;
      case DIAGONAL:
        int shou = a.getColumnSize() / 4;
        int amari = a.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = unit.create(1).divide(a.diagonalElements[j]);
        }
        for (int j = amari, counter = 0; counter < shou; ++counter, j += 4) {
          ans.diagonalElements[j] = unit.create(1).divide(a.diagonalElements[j]);// 1.0
          ans.diagonalElements[j + 1] = unit.create(1).divide(a.diagonalElements[j + 1]);
          ans.diagonalElements[j + 2] = unit.create(1).divide(a.diagonalElements[j + 2]);
          ans.diagonalElements[j + 3] = unit.create(1).divide(a.diagonalElements[j + 3]);
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @return コレスキー分解と逆行列
   */
  @SuppressWarnings("unchecked")
  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>> BlockDenseMatrix<RS,RM,CS,CM>[] getCholeskyAndInv(BlockDenseMatrix<RS,RM,CS,CM> a) {
    BlockDenseMatrix<RS,RM,CS,CM> choleskyMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
    BlockDenseMatrix<RS,RM,CS,CM> inverseMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> ans1 = getCholesky(a.getBlock(i));
      choleskyMat.setBlock(i, ans1);
      DenseMatrix<RS,RM,CS,CM> ans2 = getInvLowTriangularMatrix(ans1);
      inverseMat.setBlock(i, ans2);
    }

    return new BlockDenseMatrix[] {choleskyMat, inverseMat};
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   */
  private 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>> void getSymmetrize(DenseMatrix<RS,RM,CS,CM> a) {
    final RS unit = a.getElementUnit();

    //int index = 0;
    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        if (a.getRowSize() != a.getColumnSize()) {
          throw new RuntimeException("getSymmetrize:: different memory size"); //$NON-NLS-1$
        }
        for (int index = 0; index < a.getRowSize() - 1; ++index) {
          int index1 = index + index * a.getRowSize() + 1;
          int index2 = index + (index + 1) * a.getRowSize();
          int length = a.getRowSize() - 1 - index;
          // aMat.de_ele[index1] += aMat.de_ele[index2]

          RS[] aMat_index2 = unit.createArray(a.denseElements.length - index2);
          RS[] aMat_index1 = unit.createArray(a.denseElements.length - index1);

          System.arraycopy(a.denseElements, index1, aMat_index1, 0, aMat_index1.length);
          System.arraycopy(a.denseElements, index2, aMat_index2, 0, aMat_index2.length);
          //BLAS.daxpy(length, unit.create(1), aMat_index2, a.getRowSize(), aMat_index1, 1);
          BLAS.dxpy(length, aMat_index2, a.getRowSize(), aMat_index1, 1);
          System.arraycopy(aMat_index1, 0, a.denseElements, index1, aMat_index1.length);

          // aMat.de_ele[index1] /= 2.0

          RS half = unit.create(5).multiply(unit.create(10).power(-1));
          System.arraycopy(a.denseElements, index1, aMat_index1, 0, aMat_index1.length);
          BLAS.dscal(length, half, aMat_index1, 1);
          System.arraycopy(aMat_index1, 0, a.denseElements, index1, aMat_index1.length);

          // aMat.de_ele[index2] = aMat.de_ele[index1]
          RS[] A1 = unit.createArray(a.denseElements.length - index1);
          System.arraycopy(a.denseElements, index1, A1, 0, A1.length);
          RS[] A2 = unit.createArray(a.denseElements.length - index2);
          System.arraycopy(a.denseElements, index2, A2, 0, A2.length);
          BLAS.dcopy(length, A1, 1, A2, a.getRowSize());
          System.arraycopy(A2, 0, a.denseElements, index2, A2.length);
        }
        break;
      case DIAGONAL:
        // Nothing needs.
        break;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   */
  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>> void getSymmetrize(BlockDenseMatrix<RS,RM,CS,CM> a) {
    for (int i = 0; i < a.getBlockSize(); ++i) {
      getSymmetrize(a.getBlock(i));
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @return 転置行列
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> getTranspose(DenseMatrix<RS,RM,CS,CM> a) {
    if (a.getRowSize() != a.getColumnSize()) {
      throw new RuntimeException("getTranspose:: different memory size"); //$NON-NLS-1$
      // Of course, a non-symmetric matrix has its transposed matrix,
      // but in this algorithm we have to make transposed matrix only when symmetric matrix.
    }
    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        for (int i = 0; i < a.getRowSize(); ++i) {
          int shou = (i + 1) / 4;
          int amari = (i + 1) / 4;
          for (int j = 0; j < amari; ++j) {
            int index1 = i + a.getColumnSize() * j;
            int index2 = j + a.getColumnSize() * i;
            ans.denseElements[index1] = a.denseElements[index2];
            ans.denseElements[index2] = a.denseElements[index1];
          }
          int counter, j;
          for (j = amari, counter = 0; counter < shou; ++counter, j += 4) {
            int index1 = i + a.getColumnSize() * j;
            int index_1 = j + a.getColumnSize() * i;
            ans.denseElements[index1] = a.denseElements[index_1];
            ans.denseElements[index_1] = a.denseElements[index1];
            int index2 = i + a.getColumnSize() * (j + 1);
            int index_2 = (j + 1) + a.getColumnSize() * i;
            ans.denseElements[index2] = a.denseElements[index_2];
            ans.denseElements[index_2] = a.denseElements[index2];
            int index3 = i + a.getColumnSize() * (j + 2);
            int index_3 = (j + 2) + a.getColumnSize() * i;
            ans.denseElements[index3] = a.denseElements[index_3];
            ans.denseElements[index_3] = a.denseElements[index3];
            int index4 = i + a.getColumnSize() * (j + 3);
            int index_4 = (j + 3) + a.getColumnSize() * i;
            ans.denseElements[index4] = a.denseElements[index_4];
            ans.denseElements[index_4] = a.denseElements[index4];
          }
        }
        return ans;
      case DIAGONAL:
        throw new IllegalArgumentException();
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @return ブロック転置行列
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> getTranspose(BlockDenseMatrix<RS,RM,CS,CM> a) {
    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = getTranspose(a.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @return 計算が正常終了したならばtrue
   */
  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>> boolean choleskyFactorWithAdjust(DenseMatrix<RS,RM,CS,CM> a) {
    int rowSize = a.getRowSize();

    //int info = new Dpotrf<RS>().rATL_dpotrfL(rowSize, a.denseElements, rowSize);
    //a.setRowSize(rowSize);

    int info = Clapack.dpotrf("Lower", rowSize, a.denseElements, rowSize); //$NON-NLS-1$
    a.setRowSize(rowSize);

    if (info < 0) {
      System.err.println("cholesky argument is wrong " + (-info)); //$NON-NLS-1$
      return false;
    } else if (info > 0) {
      System.err.println("cholesky miss condition :: not positive definite" + " :: info = " + info); //$NON-NLS-1$ //$NON-NLS-2$
      return false;
    }
    return true;
  }

  /**
   * solve aMat * xVec = bVec aMat must be Cholesky Factorized. aMat must have done Cholesky factorized.
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b b
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> solveSystems(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getDimension() || a.getRowSize() != a.getColumnSize()) {
      throw new RuntimeException("solveSystems:: different memory size"); //$NON-NLS-1$
    }
    if (a.isDense() == false) {
      throw new RuntimeException("solveSystems:: matrix type must be DENSE"); //$NON-NLS-1$
    }
    Vector<RS,RM,CS,CM> x = b.createClone();
    BLAS.dtrsv("Lower", "NoTranspose", "NonUnit", a.getRowSize(), a.denseElements, a.getColumnSize(), x.getElements(), 1); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
    BLAS.dtrsv("Lower", "Transpose", "NonUnit", a.getRowSize(), a.denseElements, a.getColumnSize(), x.getElements(), 1); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
    return x;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getColumnSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
          ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
          ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
          ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
          ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
          ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
          ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * Dense=Sparse*Dense
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = b.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal(), unit);

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }
        for (int index = 0; index < a.nonZeroCount; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index).multiply(scalar);
          if (i != j) {
            RS[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            RS[] retMat_i = unit.createArray(ans.denseElements.length - i);
            System.arraycopy(ans.denseElements, i, retMat_i, 0, retMat_i.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_i, ans.getRowSize());
            System.arraycopy(retMat_i, 0, ans.denseElements, i, retMat_i.length);

            RS[] bMat_i = unit.createArray(b.denseElements.length - i);
            System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
            System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_i, b.getRowSize(), retMat_j, ans.getRowSize());
            System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
          } else {
            RS[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
            System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_j, ans.getRowSize());
            System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || b.isDiagonal() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = scalar.multiply(a.getDiagonalElement(j).multiply(b.diagonalElements[j]));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = scalar.multiply(a.getDiagonalElement(j)).multiply(b.diagonalElements[j]);
          ans.diagonalElements[j + 1] = scalar.multiply(a.getDiagonalElement(j + 1).multiply(b.diagonalElements[j + 1]));
          ans.diagonalElements[j + 2] = scalar.multiply(a.getDiagonalElement(j + 2).multiply(b.diagonalElements[j + 2]));
          ans.diagonalElements[j + 3] = scalar.multiply(a.getDiagonalElement(j + 3).multiply(b.diagonalElements[j + 3]));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * TODO TODO Dense=Sparse*Dense
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = b.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal(), unit);

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        for (int index = 0; index < a.nonZeroCount; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index);
          if (i != j) {
            RS[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            RS[] retMat_i = unit.createArray(ans.denseElements.length - i);
            System.arraycopy(ans.denseElements, i, retMat_i, 0, retMat_i.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_i, ans.getRowSize());
            System.arraycopy(retMat_i, 0, ans.denseElements, i, retMat_i.length);

            RS[] bMat_i = unit.createArray(b.denseElements.length - i);
            System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
            System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_i, b.getRowSize(), retMat_j, ans.getRowSize());
            System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
          } else {
            RS[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
            System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_j, ans.getRowSize());
            System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || b.isDiagonal() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.getDiagonalElement(j).multiply(b.diagonalElements[j]);
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.getDiagonalElement(j).multiply(b.diagonalElements[j]);
          ans.diagonalElements[j + 1] = a.getDiagonalElement(j + 1).multiply(b.diagonalElements[j + 1]);
          ans.diagonalElements[j + 2] = a.getDiagonalElement(j + 2).multiply(b.diagonalElements[j + 2]);
          ans.diagonalElements[j + 3] = a.getDiagonalElement(j + 3).multiply(b.diagonalElements[j + 3]);
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * TODO MULTIPLY_NON_ATLASがFalseの場合でやるのもいいかも。 falseの場合でやると、配列のコピーの回数がここの分だけ減る。 Dense=Dense*Sparse
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (b.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        for (int index = 0; index < b.nonZeroCount; ++index) {
          int i = b.rowIndex[index];
          int j = b.columnIndex[index];
          RS value = b.getSparseElement(index).multiply(scalar);
          if (i != j) {
            RS[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
            System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
            System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_i, 1, retMat_j, 1);
            System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);

            RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            RS[] retMat_i = unit.createArray(ans.denseElements.length - ans.getRowSize() * i);
            System.arraycopy(ans.denseElements, ans.getRowSize() * i, retMat_i, 0, retMat_i.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_i, 1);
            System.arraycopy(retMat_i, 0, ans.denseElements, ans.getRowSize() * i, retMat_i.length);
          } else {
            RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_j, 1);
            System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || a.isDiagonal() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.getDiagonalElement(j)));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.getDiagonalElement(j)));
          ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.getDiagonalElement(j + 1)));
          ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.getDiagonalElement(j + 2)));
          ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.getDiagonalElement(j + 3)));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * TODO MULTIPLY_NON_ATLASがFalseの場合でやるのもいいかも。 falseの場合でやると、配列のコピーの回数がここの分だけ減る。 Dense=Dense*Sparse
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (b.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }

        for (int index = 0; index < b.nonZeroCount; ++index) {
          int i = b.rowIndex[index];
          int j = b.columnIndex[index];
          RS value = b.getSparseElement(index);
          if (i != j) {
            RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            RS[] retMat_i = unit.createArray(ans.denseElements.length - ans.getRowSize() * i);
            System.arraycopy(ans.denseElements, ans.getRowSize() * i, retMat_i, 0, retMat_i.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_i, 1);
            System.arraycopy(retMat_i, 0, ans.denseElements, ans.getRowSize() * i, retMat_i.length);

            RS[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
            System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
            System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_i, 1, retMat_j, 1);
            System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
          } else {
            RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
            BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_j, 1);
            System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dgemm(
            "NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || a.isDiagonal() == false) {
          throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
        }
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.getDiagonalElement(j));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.getDiagonalElement(j));
          ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.getDiagonalElement(j + 1));
          ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.getDiagonalElement(j + 2));
          ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.getDiagonalElement(j + 3));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, RS scalar) {
    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dscal(length, scalar, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        BLAS.dscal(ans.getColumnSize(), scalar, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, RS scalar) {
    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param scalar スカラー
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> multiply(Vector<RS,RM,CS,CM> a, RS scalar) {
    Vector<RS,RM,CS,CM> ans = a.createClone();
    BLAS.dscal(ans.getDimension(), scalar, ans.getElements(), 1);
    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockVector<RS,RM,CS,CM> multiply(BlockVector<RS,RM,CS,CM> a, RS scalar) {
    BlockVector<RS,RM,CS,CM> ans = a.createClone();

    for (int i = 0; i < a.getBlockSize(); ++i) {
      Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b b
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
    if (a.getColumnSize() != b.getDimension()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    Vector<RS,RM,CS,CM> ans = b.createClone();

    final RS unit = a.getElementUnit();

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemv("NoTranspose", a.getRowSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.getElements(), 1, unit.create(0), ans.getElements(), 1); //$NON-NLS-1$
        return ans;
      case DIAGONAL:
        int shou = ans.getDimension() / 4;
        int amari = ans.getDimension() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.setElement(j, a.diagonalElements[j].multiply(b.getElement(j)));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.setElement(j, a.diagonalElements[j].multiply(b.getElement(j)));
          ans.setElement(j + 1, a.diagonalElements[j + 1].multiply(b.getElement(j + 1)));
          ans.setElement(j + 2, a.diagonalElements[j + 2].multiply(b.getElement(j + 2)));
          ans.setElement(j + 3, a.diagonalElements[j + 3].multiply(b.getElement(j + 3)));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b, RS scalar) {
    if (a.getColumnSize() != b.getDimension()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    Vector<RS,RM,CS,CM> ans = new Vector<>(a.getRowSize(), unit);

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemv("NoTranspose", a.getRowSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.getElements(), 1, unit.create(0), ans.getElements(), 1); //$NON-NLS-1$
        return ans;
      case DIAGONAL:
        int shou = ans.getDimension() / 4;
        int amari = ans.getDimension() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.setElement(j, scalar.multiply(a.diagonalElements[j].multiply(b.getElement(j))));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.setElement(j, scalar.multiply(a.diagonalElements[j].multiply(b.getElement(j))));
          ans.setElement(j + 1, scalar.multiply(a.diagonalElements[j + 1].multiply(b.getElement(j + 1))));
          ans.setElement(j + 2, scalar.multiply(a.diagonalElements[j + 2].multiply(b.getElement(j + 2))));
          ans.setElement(j + 3, scalar.multiply(a.diagonalElements[j + 3].multiply(b.getElement(j + 3))));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockVector<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

    for (int i = 0; i < a.getBlockSize(); ++i) {
      Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b b
   * @return 解
   */
  private 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>> BlockVector<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
    if (b.getBlockSize() != a.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockVector<RS,RM,CS,CM> ans = b.createClone();

    for (int i = 0; i < b.getBlockSize(); ++i) {
      Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = aMat**T * bMat
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> tran_multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemm(
            "Transpose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getColumnSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:

        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
          ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
          ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
          ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = aMat**T * bMat
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> tran_multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
    
    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        BLAS.dgemm(
            "Transpose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getColumnSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
          ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
          ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
          ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> tran_multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = tran_multiply(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> tran_multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = tran_multiply(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = aMat * bMat**T
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply_tran(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        // The Point is the first argument is "NoTranspose".
        BLAS.dgemm(
            "NoTranspose", "Transpose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getColumnSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
          ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
          ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
          ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = aMat * bMat**T
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> multiply_tran(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getColumnSize() != b.getColumnSize()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final RS unit = a.getElementUnit();
    DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal(), unit);

    switch (ans.getDenseOrDiagonal()) {
      case DENSE:
        // The Point is the first argument is "NoTranspose".
        BLAS.dgemm(
            "NoTranspose", "Transpose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getColumnSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
            ans.getRowSize());
        return ans;
      case DIAGONAL:
        int shou = ans.getColumnSize() / 4;
        int amari = ans.getColumnSize() % 4;
        for (int j = 0; j < amari; ++j) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
        }
        int count;
        int j;
        for (j = amari, count = 0; count < shou; ++count, j += 4) {
          ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
          ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
          ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
          ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
        }
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply_tran(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = multiply_tran(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> multiply_tran(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    RS scalar = null;
    BlockDenseMatrix<RS,RM,CS,CM> ans = multiply_tran(a, b, scalar);
    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> plus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b, RS scalar) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    Vector<RS,RM,CS,CM> ans = a.createClone();
    BLAS.daxpy(ans.getDimension(), scalar, b.getElements(), 1, ans.getElements(), 1);
    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> plus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    Vector<RS,RM,CS,CM> ans = a.createClone();
    BLAS.dxpy(ans.getDimension(), b.getElements(), 1, ans.getElements(), 1);
    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @return 解
   */
  private 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>> Vector<RS,RM,CS,CM> minus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    Vector<RS,RM,CS,CM> ans = a.createClone();
    BLAS.dxmy(ans.getDimension(), b.getElements(), 1, ans.getElements(), 1);
    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        int length = a.getRowSize() * a.getColumnSize();
        BLAS.daxpy(length, scalar, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        BLAS.daxpy(ans.getColumnSize(), scalar, b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    final DenseMatrix<RS,RM,CS,CM> ans;

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        int length = a.getRowSize() * a.getColumnSize();
        ans = a.createClone();
        BLAS.dxpy(length, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        ans = a.createClone();
        BLAS.dxpy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> minus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    final DenseMatrix<RS,RM,CS,CM> ans;

    switch (a.getDenseOrDiagonal()) {
      case DENSE:
        int length = a.getRowSize() * a.getColumnSize();
        ans = a.createClone();
        BLAS.dxmy(length, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        ans = a.createClone();
        BLAS.dxmy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ret = (*scalar) * b
    DenseMatrix<RS,RM,CS,CM> ans = multiply(b, scalar);

    // ret += a

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }

        int shou = a.nonZeroCount / 4;
        int amari = a.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].add(value);
            ans.denseElements[ji] = ans.denseElements[ji].add(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].add(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = a.rowIndex[index];
          int j1 = a.columnIndex[index];
          RS value1 = a.getSparseElement(index);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].add(value1);
            ans.denseElements[ji] = ans.denseElements[ji].add(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].add(value1);
          }
          int i2 = a.rowIndex[index + 1];
          int j2 = a.columnIndex[index + 1];
          RS value2 = a.getSparseElement(index + 1);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].add(value2);
            ans.denseElements[ji] = ans.denseElements[ji].add(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].add(value2);
          }
          int i3 = a.rowIndex[index + 2];
          int j3 = a.columnIndex[index + 2];
          RS value3 = a.getSparseElement(index + 2);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].add(value3);
            ans.denseElements[ji] = ans.denseElements[ji].add(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].add(value3);
          }
          int i4 = a.rowIndex[index + 3];
          int j4 = a.columnIndex[index + 3];
          RS value4 = a.getSparseElement(index + 3);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].add(value4);
            ans.denseElements[ji] = ans.denseElements[ji].add(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].add(value4);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || b.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dxpy(length, a.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || b.isDiagonal() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dxpy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<RS,RM,CS,CM> ans = b.createClone();

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (b.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }

        int shou = a.nonZeroCount / 4;
        int amari = a.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].add(value);
            ans.denseElements[ji] = ans.denseElements[ji].add(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].add(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = a.rowIndex[index];
          int j1 = a.columnIndex[index];
          RS value1 = a.getSparseElement(index);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].add(value1);
            ans.denseElements[ji] = ans.denseElements[ji].add(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].add(value1);
          }
          int i2 = a.rowIndex[index + 1];
          int j2 = a.columnIndex[index + 1];
          RS value2 = a.getSparseElement(index + 1);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].add(value2);
            ans.denseElements[ji] = ans.denseElements[ji].add(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].add(value2);
          }
          int i3 = a.rowIndex[index + 2];
          int j3 = a.columnIndex[index + 2];
          RS value3 = a.getSparseElement(index + 2);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].add(value3);
            ans.denseElements[ji] = ans.denseElements[ji].add(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].add(value3);
          }
          int i4 = a.rowIndex[index + 3];
          int j4 = a.columnIndex[index + 3];
          RS value4 = a.getSparseElement(index + 3);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].add(value4);
            ans.denseElements[ji] = ans.denseElements[ji].add(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].add(value4);
          }
        }
        return ans;
      case DENSE:
        if (b.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dxpy(length, a.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        if (b.isDiagonal() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dxpy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> minus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<RS,RM,CS,CM> ans = b.createClone();

    switch (a.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (b.isDense() == false) {
          throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
        }

        int shou = a.nonZeroCount / 4;
        int amari = a.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = a.rowIndex[index];
          int j = a.columnIndex[index];
          RS value = a.getSparseElement(index);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = a.rowIndex[index];
          int j1 = a.columnIndex[index];
          RS value1 = a.getSparseElement(index);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value1);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value1);
          }
          int i2 = a.rowIndex[index + 1];
          int j2 = a.columnIndex[index + 1];
          RS value2 = a.getSparseElement(index + 1);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value2);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value2);
          }
          int i3 = a.rowIndex[index + 2];
          int j3 = a.columnIndex[index + 2];
          RS value3 = a.getSparseElement(index + 2);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value3);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value3);
          }
          int i4 = a.rowIndex[index + 3];
          int j4 = a.columnIndex[index + 3];
          RS value4 = a.getSparseElement(index + 3);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value4);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value4);
          }
        }
        return ans;
      case DENSE:
        if (b.isDense() == false) {
          throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
        }
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dxmy(length, a.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        if (b.isDiagonal() == false) {
          throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.dxmy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ret = a
    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    //ret += (*scalar) * b

    switch (b.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        int shou = b.nonZeroCount / 4;
        int amari = b.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = b.rowIndex[index];
          int j = b.columnIndex[index];
          RS value = b.getSparseElement(index).multiply(scalar);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].add(value);
            ans.denseElements[ji] = ans.denseElements[ji].add(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].add(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = b.rowIndex[index];
          int j1 = b.columnIndex[index];
          RS value1 = b.getSparseElement(index).multiply(scalar);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].add(value1);
            ans.denseElements[ji] = ans.denseElements[ji].add(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].add(value1);
          }
          int i2 = b.rowIndex[index + 1];
          int j2 = b.columnIndex[index + 1];
          RS value2 = b.getSparseElement(index + 1).multiply(scalar);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].add(value2);
            ans.denseElements[ji] = ans.denseElements[ji].add(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].add(value2);
          }
          int i3 = b.rowIndex[index + 2];
          int j3 = b.columnIndex[index + 2];
          RS value3 = b.getSparseElement(index + 2).multiply(scalar);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].add(value3);
            ans.denseElements[ji] = ans.denseElements[ji].add(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].add(value3);
          }
          int i4 = b.rowIndex[index + 3];
          int j4 = b.columnIndex[index + 3];
          RS value4 = b.getSparseElement(index + 3).multiply(scalar);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].add(value4);
            ans.denseElements[ji] = ans.denseElements[ji].add(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].add(value4);
          }
        }
        return ans;
      case DENSE:
        if (ans.isDense() == false || a.isDense() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.daxpy(length, scalar, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        if (ans.isDiagonal() == false || a.isDiagonal() == false) {
          throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
        }
        BLAS.daxpy(ans.getColumnSize(), scalar, b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ans = a
    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    // ans += (*scalar) * b

    switch (b.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        int shou = b.nonZeroCount / 4;
        int amari = b.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = b.rowIndex[index];
          int j = b.columnIndex[index];
          RS value = b.getSparseElement(index);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].add(value);
            ans.denseElements[ji] = ans.denseElements[ji].add(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].add(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = b.rowIndex[index];
          int j1 = b.columnIndex[index];
          RS value1 = b.getSparseElement(index);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].add(value1);
            ans.denseElements[ji] = ans.denseElements[ji].add(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].add(value1);
          }
          int i2 = b.rowIndex[index + 1];
          int j2 = b.columnIndex[index + 1];
          RS value2 = b.getSparseElement(index + 1);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].add(value2);
            ans.denseElements[ji] = ans.denseElements[ji].add(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].add(value2);
          }
          int i3 = b.rowIndex[index + 2];
          int j3 = b.columnIndex[index + 2];
          RS value3 = b.getSparseElement(index + 2);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].add(value3);
            ans.denseElements[ji] = ans.denseElements[ji].add(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].add(value3);
          }
          int i4 = b.rowIndex[index + 3];
          int j4 = b.columnIndex[index + 3];
          RS value4 = b.getSparseElement(index + 3);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].add(value4);
            ans.denseElements[ji] = ans.denseElements[ji].add(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].add(value4);
          }
        }
        return ans;
      case DENSE:
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dxpy(length, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        BLAS.dxpy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> DenseMatrix<RS,RM,CS,CM> minus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<RS,RM,CS,CM> ans = a.createClone();

    switch (b.getSparseOrDenseOrDiagonal()) {
      case SPARSE:
        int shou = b.nonZeroCount / 4;
        int amari = b.nonZeroCount % 4;
        for (int index = 0; index < amari; ++index) {
          int i = b.rowIndex[index];
          int j = b.columnIndex[index];
          RS value = b.getSparseElement(index);
          if (i != j) {
            int ij = i + ans.getColumnSize() * j;
            int ji = j + ans.getColumnSize() * i;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value);
          } else {
            int ii = i + ans.getColumnSize() * i;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value);
          }
        }

        int index;
        int counter;
        for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
          int i1 = b.rowIndex[index];
          int j1 = b.columnIndex[index];
          RS value1 = b.getSparseElement(index);
          if (i1 != j1) {
            int ij = i1 + ans.getColumnSize() * j1;
            int ji = j1 + ans.getColumnSize() * i1;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value1);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value1);
          } else {
            int ii = i1 + ans.getColumnSize() * i1;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value1);
          }
          int i2 = b.rowIndex[index + 1];
          int j2 = b.columnIndex[index + 1];
          RS value2 = b.getSparseElement(index + 1);
          if (i2 != j2) {
            int ij = i2 + ans.getColumnSize() * j2;
            int ji = j2 + ans.getColumnSize() * i2;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value2);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value2);
          } else {
            int ii = i2 + ans.getColumnSize() * i2;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value2);
          }
          int i3 = b.rowIndex[index + 2];
          int j3 = b.columnIndex[index + 2];
          RS value3 = b.getSparseElement(index + 2);
          if (i3 != j3) {
            int ij = i3 + ans.getColumnSize() * j3;
            int ji = j3 + ans.getColumnSize() * i3;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value3);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value3);
          } else {
            int ii = i3 + ans.getColumnSize() * i3;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value3);
          }
          int i4 = b.rowIndex[index + 3];
          int j4 = b.columnIndex[index + 3];
          RS value4 = b.getSparseElement(index + 3);
          if (i4 != j4) {
            int ij = i4 + ans.getColumnSize() * j4;
            int ji = j4 + ans.getColumnSize() * i4;
            ans.denseElements[ij] = ans.denseElements[ij].subtract(value4);
            ans.denseElements[ji] = ans.denseElements[ji].subtract(value4);
          } else {
            int ii = i4 + ans.getColumnSize() * i4;
            ans.denseElements[ii] = ans.denseElements[ii].subtract(value4);
          }
        }
        return ans;
      case DENSE:
        int length = ans.getRowSize() * ans.getColumnSize();
        BLAS.dxmy(length, b.denseElements, 1, ans.denseElements, 1);
        return ans;
      case DIAGONAL:
        BLAS.dxmy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
        return ans;
      default:
        throw new IllegalArgumentException();
    }
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockVector<RS,RM,CS,CM> plus(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    RS unit = a.getElementUnit();
    BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

    for (int i = 0; i < a.getBlockSize(); ++i) {
      Vector<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param b b
   * @return 解
   */
  private 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>> BlockVector<RS,RM,CS,CM> plus(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    RS unit = a.getElementUnit();
    BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

    for (int i = 0; i < a.getBlockSize(); ++i) {
      Vector<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans + a * b*scalar
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param b B
   * @return 解
   */
  private 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>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());

    for (int i = 0; i < a.getBlockSize(); ++i) {
      DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
      ans.setBlock(i, blockAns);
    }

    return ans;
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  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>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, RS scalar) {
    switch (operator) {
      case '*':
        Vector<RS,RM,CS,CM> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  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>> BlockVector<RS,RM,CS,CM> let(BlockVector<RS,RM,CS,CM> a, final char operator, RS scalar) {
    switch (operator) {
      case '*':
        BlockVector<RS,RM,CS,CM> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, RS scalar) {
    switch (operator) {
      case '*':
        BlockDenseMatrix<RS,RM,CS,CM> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  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>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b, RS scalar) {
    Vector<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  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>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
    Vector<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' 't' 'T' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    DenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      case 't':
        // ret = aMat**T * bMat
        ans = tran_multiply(a, b, scalar);
        return ans;
      case 'T':
        // ret = aMat * bMat**T
        ans = multiply_tran(a, b, scalar);
        return ans;
      default:
        throw new IllegalArgumentException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' 't' 'T' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
    DenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      case 't':
        // ret = aMat**T * bMat
        ans = tran_multiply(a, b);
        return ans;
      case 'T':
        // ret = aMat * bMat**T
        ans = multiply_tran(a, b);
        return ans;
      default:
        throw new IllegalArgumentException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    DenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
    DenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    DenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b) {
    DenseMatrix<RS,RM,CS,CM> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' 't' 'T' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      case 't':
        // ret = aMat**T * bMat
        ans = tran_multiply(a, b, scalar);
        return ans;
      case 'T':
        // ret = aMat * bMat**T
        ans = multiply_tran(a, b, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' 't' 'T' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      case 't':
        // ret = aMat**T * bMat
        ans = tran_multiply(a, b);
        return ans;
      case 'T':
        // ret = aMat * bMat**T
        ans = multiply_tran(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final RS minus_scalar = scalar.unaryMinus();
        ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' '*' b*(*scalar)
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b) {
    BlockDenseMatrix<RS,RM,CS,CM> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final RS unit = a.getElementUnit();
        //final RS minus_scalar = unit.create(-1);
        //ans = plus(a, b, minus_scalar);
        return ans;
      case '*':
        ans = multiply(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = aMat '*' '/' bVec
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  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>> Vector<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
    Vector<RS,RM,CS,CM> ans;
    switch (operator) {
      case '*':
        ans = multiply(a, b);
        return ans;
      case '/':
        // ret = aMat^{-1} * bVec;
        // aMat is positive definite and already cholesky factorized.
        ans = solveSystems(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  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>> RS run(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(b, a);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  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>> RS run(BlockVector<RS,RM,CS,CM> a, final char operator, BlockVector<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  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>> RS run(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b) {
    switch (operator) {
      case '.':
        RS ans = getInnerProduct(b, a);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }
}