BlockDenseMatrix.java

package org.mklab.sdpj.datastructure;

import java.io.PrintStream;
import java.util.Arrays;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;

/**
 * ブロック密行列を表すクラスです。
 * 
 * @author takafumi
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
 */
public class BlockDenseMatrix<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>> {

  /** ブロック数 */
  private int blockSize;
  /** ブロック構造 */
  private int[] blockStruct;
  /** 成分 */
  private DenseMatrix<RS,RM,CS,CM>[] blocks;
  
  /** */
  private RS unit;

  /**
   * 新しく生成された{@link BlockDenseMatrix}オブジェクトを初期化します。
   * @param unit unit
   */
  public BlockDenseMatrix(RS unit) {
    this.unit = unit;
    this.blockSize = 0;
    this.blockStruct = null;
    this.blocks = null;
  }

  /**
   * 新しく生成された{@link BlockDenseMatrix}オブジェクトを初期化します。
   * 
   * @param blockSize ブロック数
   * @param blockStruct ブロック構造
   * @param unit unit
   */
  public BlockDenseMatrix(int blockSize, int[] blockStruct, RS unit) {
    this.unit = unit;
    this.initialize(blockSize, blockStruct);
  }

  /**
   * Returns unit.
   * 
   * @return unit
   */
  public RS getUnit() {
    return this.unit;
  }
  
  /**
   * 初期化します。
   * 
   * @param blockSize1 ブロック数
   * @param blockStruct1 ブロック構造
   */
  public void initialize(int blockSize1, int[] blockStruct1) {
    this.blockSize = blockSize1;

    if (blockSize1 < 0) {
      throw new RuntimeException("BlockDenseMatrix::nBlock is nonpositive"); //$NON-NLS-1$
    }

    this.blockStruct = new int[blockSize1];
    for (int i = 0; i < this.blockSize; i++) {
      this.blockStruct[i] = blockStruct1[i];
    }

    this.blocks = new DenseMatrix[blockSize1];
    for (int i = 0; i < blockSize1; i++) {
      int size = blockStruct1[i];
      if (size > 0) {
        this.blocks[i] = new DenseMatrix<>(size, size, DenseDiagonal.DENSE, this.unit);
      } else {
        size = (-size);
        this.blocks[i] = new DenseMatrix<>(size, size, DenseDiagonal.DIAGONAL, this.unit);
      }
    }
  }

  /**
   * This function is to compute the condition number, minimum singular, and
   * maximum singular value of matrix.
   * 
   * @param out Print result of condition number
   * @param min Print minimum singular value of matrix
   * @param max Print maximum singular value of matrix
   */
  public void getCondtionNumber(PrintStream out, PrintStream min, PrintStream max) {
    //RS unit = Tools.getUnitNumber();

    int rowSize = getRowSize();
    int columnSize = getColumnSize();

    RS zero = this.unit.createZero();
    RS[][] data = this.unit.createArray(rowSize, columnSize);
    for (int i = 0; i < rowSize; ++i) {
      for (int j = 0; j < columnSize; ++j) {
        data[i][j] = zero;
      }
    }

    int shift = 0;
    for (int k = 0; k < this.blockSize; ++k) {
      if (k > 0) {
        shift += this.blocks[k - 1].getColumnSize();
      }

      DenseMatrix<RS,RM,CS,CM> matrix = this.blocks[k];
      int rowNumber = matrix.getRowSize();
      int columnNumber = matrix.getColumnSize();

      switch (matrix.getDenseOrDiagonal()) {
        case DENSE:
          for (int i = 0; i < rowNumber; ++i) {
            for (int j = 0; j < columnNumber; ++j) {
              data[i + shift][j + shift] = matrix.denseElements[i + columnNumber * j];
            }
          }
          break;
        case DIAGONAL:
          for (int j = 0; j < columnNumber; ++j) {
            data[j + shift][j + shift] = matrix.diagonalElements[j];
          }
          break;
        default:
          throw new UnsupportedOperationException();
      }
    }
    //NumericalMatrix<RS> matrix = new BaseNumericalMatrix<>(columnSize, columnSize, data);
    RM matrix =  data[0][0].createGrid(data);
    out.printf("%s\n", matrix.conditionNumber().toString("%4.3e")); //$NON-NLS-1$ //$NON-NLS-2$
    min.printf("%s \n", matrix.minSingularValue().toString("%4.2e")); //$NON-NLS-1$ //$NON-NLS-2$
    max.printf("%s \n", matrix.maxSingularValue().toString("%4.2e")); //$NON-NLS-1$ //$NON-NLS-2$
  }

  /**
   * 出力します。
   * 
   * @param out 出力先
   */
  @SuppressWarnings("boxing")
  public void display(PrintStream out) {
    out.print("{\n"); //$NON-NLS-1$

    for (int k = 0; k < this.blockSize; ++k) {
      DenseMatrix<RS,RM,CS,CM> matrix = this.blocks[k];
      int rowSize = matrix.getRowSize();
      int columnSize = matrix.getColumnSize();

      switch (matrix.getDenseOrDiagonal()) {
        case DENSE:
          out.print("{"); //$NON-NLS-1$
          for (int i = 0; i < rowSize - 1; ++i) {
            if (i != 0) {
              out.print(" "); //$NON-NLS-1$
            }
            out.print(" "); //$NON-NLS-1$
            out.print("{ "); //$NON-NLS-1$
            for (int j = 0; j < columnSize - 1; ++j) {
              out.printf("%s, ", matrix.denseElements[i + columnSize * j].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
            }
            out.printf("%s },\n", matrix.denseElements[i + columnSize * (columnSize - 1)].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }
          if (rowSize > 1) {
            out.print("  { "); //$NON-NLS-1$
          }
          for (int j = 0; j < columnSize - 1; ++j) {
            out.printf("%s, ", matrix.denseElements[(rowSize - 1) + columnSize * j].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }
          out.printf("%s }", matrix.denseElements[(rowSize - 1) + columnSize * (columnSize - 1)].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$

          if (rowSize > 1) {
            out.print("   }"); //$NON-NLS-1$
          }
          out.print("\n"); //$NON-NLS-1$
          break;

        case DIAGONAL:
          out.print("{"); //$NON-NLS-1$
          for (int j = 0; j < columnSize - 1; ++j) {
            out.printf("%+4.6e, ", Double.parseDouble(matrix.diagonalElements[j].toString())); //$NON-NLS-1$
          }
          if (columnSize > 0) {
            out.printf("%+4.6e }\n", Double.parseDouble(matrix.diagonalElements[columnSize - 1].toString())); //$NON-NLS-1$
          }
          break;
        default:
          throw new UnsupportedOperationException();
      }
    }

    out.print("} \n"); //$NON-NLS-1$
  }

  /**
   * Matlab形式で出力します。
   * 
   * @param out 出力先
   */
  public void displayAsMatlabFormat(PrintStream out) {
    out.print("[\n"); //$NON-NLS-1$

    for (int k = 0; k < this.blockSize; ++k) {
      DenseMatrix<RS,RM,CS,CM> matrix = this.blocks[k];
      int rowSize = matrix.getRowSize();
      int columnSize = matrix.getColumnSize();

      switch (matrix.getDenseOrDiagonal()) {
        case DENSE:
          for (int i = 0; i < rowSize - 1; ++i) {
            for (int j = 0; j < columnSize - 1; ++j) {
              out.printf("%s;\n", matrix.denseElements[i + columnSize * j].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
            }
            out.printf("%s;\n", matrix.denseElements[i + columnSize * (columnSize - 1)].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }

          for (int j = 0; j < columnSize - 1; ++j) {
            out.printf("%s;\n", matrix.denseElements[(rowSize - 1) + columnSize * j].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }
          out.printf("%s;\n", matrix.denseElements[(rowSize - 1) + columnSize * (columnSize - 1)].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          break;
        case DIAGONAL:
          for (int j = 0; j < columnSize - 1; ++j) {
            out.printf("%s;\n", matrix.diagonalElements[j].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }
          if (columnSize > 0) {
            out.printf("%s;\n", matrix.diagonalElements[columnSize - 1].toString("%+4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
          }
          break;
        default:
          throw new UnsupportedOperationException();
      }
    }

    out.print("];\n\n"); //$NON-NLS-1$

  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {
    StringBuffer buff = new StringBuffer();
    buff.append("{\n"); //$NON-NLS-1$

    for (int i = 0; i < this.blockSize; i++) {
      buff.append(this.blocks[i].toString());
    }

    buff.append("} \n"); //$NON-NLS-1$
    return buff.toString();
  }

  /**
   * コピーします。
   * 
   * @param source コピー元
   */
  public void copyFrom(BlockDenseMatrix<RS,RM,CS,CM> source) {
    if (source.blockSize < 0) {
      throw new RuntimeException("rBlockDenseMatix:: nBlock is nonpositive"); //$NON-NLS-1$
    }

    this.blockSize = source.blockSize;
    this.blockStruct = new int[source.blockSize];

    for (int i = 0; i < source.blockSize; i++) {
      this.blockStruct[i] = source.blockStruct[i];
    }

    this.blocks = new DenseMatrix[source.blockSize];
    for (int i = 0; i < source.blockSize; i++) {
      this.blocks[i] = source.blocks[i].createClone();
    }
  }

  /**
   * ゼロ行列に設定します。
   */
  public void setZero() {
    for (int i = 0; i < this.blockSize; i++) {
      this.blocks[i].setZero();
    }
  }

  /**
   * 単位行列に設定します。
   * 
   * @param unit 単位
   */
  public void setIdentity(RS unit) {
    for (int i = 0; i < this.blockSize; i++) {
      this.blocks[i].setIdentity(unit);
    }
  }

  /**
   * ブロック数を返します。
   * 
   * @return ブロック数
   */
  public int getBlockSize() {
    return this.blockSize;
  }

  /**
   * ブロックの構造を返します。
   * 
   * @param index 指数
   * @return ブロックの構造
   */
  public int getBlockStruct(final int index) {
    return this.blockStruct[index];
  }

  /**
   * ブロックの構造を返します。
   * 
   * @return ブロックの構造
   */
  public int[] getBlockStructs() {
    return Arrays.copyOf(this.blockStruct, this.blockSize);
  }

  /**
   * ブロック成分を返します。
   * 
   * @param index 指数
   * @return　ブロック成分
   */
  public DenseMatrix<RS,RM,CS,CM> getBlock(int index) {
    return this.blocks[index];
  }

  /**
   * ブロック成分を設定します。
   * 
   * @param index 指数
   * @param block ブロック成分
   */
  public void setBlock(int index, DenseMatrix<RS,RM,CS,CM> block) {
    this.blocks[index] = block;
  }

  /**
   * 成分の単位を返します。
   * 
   * @return　成分の単位
   */
  public RS getElementUnit() {
    return this.blocks[0].getElementUnit();
  }

  /**
   * Returns a copy of instance of class <code>BlockDenseMatrix</code>.
   * 
   * @return a copy of instance
   */
  public BlockDenseMatrix<RS,RM,CS,CM> createClone() {
    BlockDenseMatrix<RS,RM,CS,CM> inst = new BlockDenseMatrix<>(this.unit);
    inst.blockSize = this.blockSize;

    if (this.blockStruct != null) {
      inst.blockStruct = new int[this.blockStruct.length];
      for (int i = 0; i < this.blockStruct.length; i++) {
        inst.blockStruct[i] = this.blockStruct[i];
      }
    } else {
      inst.blockStruct = null;
    }

    if (this.blocks != null) {
      inst.blocks = new DenseMatrix[this.blocks.length];
      for (int i = 0; i < this.blocks.length; i++) {
        inst.blocks[i] = this.blocks[i] == null ? null : (DenseMatrix<RS,RM,CS,CM>)this.blocks[i].createClone();
      }
    } else {
      inst.blocks = null;
    }
    return inst;
  }

  /**
   * 行数を返します。
   * 
   * @return 行数
   */
  public int getRowSize() {
    int rowSize = 0;
    for (DenseMatrix<RS,RM,CS,CM> matrix : this.blocks) {
      rowSize += matrix.getRowSize();
    }
    return rowSize;
  }

  /**
   * 列数を返します。
   * 
   * @return 列数
   */
  public int getColumnSize() {
    int columnSize = 0;
    for (DenseMatrix<RS,RM,CS,CM> matrix : this.blocks) {
      columnSize += matrix.getColumnSize();
    }
    return columnSize;
  }

  /**
   * Return an array of elements.
   * 
   * @return array of elements
   */
  public RS[][] getArrayOfElements() {
    final int rowSize = getRowSize();
    final int columnSize = getColumnSize();
    //final RS unit = getElementUnit();
    final RS[][] data = this.unit.createArray(rowSize, columnSize);

    RS zero = this.unit.createZero();
    for (int i = 0; i < rowSize; i++) {
      for (int j = 0; j < columnSize; j++) {
        data[i][j] = zero;
      }
    }

    int shift = 0;
    for (int k = 0; k < this.blockSize; ++k) {
      if (k > 0) {
        shift += getBlock(k - 1).getColumnSize();
      }

      final DenseMatrix<RS,RM,CS,CM> matrix = getBlock(k);
      final int rowNumber = matrix.getRowSize();
      final int columnNumber = matrix.getColumnSize();

      switch (matrix.getDenseOrDiagonal()) {
        case DENSE:
          for (int i = 0; i < rowNumber; ++i) {
            for (int j = 0; j < columnNumber; ++j) {
              data[i + shift][j + shift] = matrix.denseElements[i + columnNumber * j];
            }
          }
          break;
        case DIAGONAL:
          for (int j = 0; j < columnNumber; ++j) {
            data[j + shift][j + shift] = matrix.diagonalElements[j];
          }
          break;
        default:
          throw new UnsupportedOperationException();
      }
    }

    return data;
  }

  /**
   * Set elements
   * 
   * @param a matrix
   */
  public void setElements(RM a) {
    if (getRowSize() != a.getRowSize() || getColumnSize() != a.getColumnSize()) {
      throw new IllegalArgumentException();
    }

    int shift = 0;
    for (int k = 0; k < getBlockSize(); ++k) {
      if (k > 0) {
        shift += getBlock(k - 1).getColumnSize();
      }

      DenseMatrix<RS,RM,CS,CM> matrix = getBlock(k);
      int rowSize = matrix.getRowSize();
      int columnSize = matrix.getColumnSize();

      switch (matrix.getDenseOrDiagonal()) {
        case DENSE:
          for (int i = 0; i < rowSize; ++i) {
            for (int j = 0; j < columnSize; ++j) {
              matrix.denseElements[i + columnSize * j] = a.getElement(i + shift + 1, j + shift + 1);
            }
          }
          break;
        case DIAGONAL:
          for (int i = 0; i < columnSize; ++i) {
            matrix.diagonalElements[i] = a.getElement(i + shift, i + shift);
          }
          break;
        default:
          throw new UnsupportedOperationException();
      }
    }
  }

}