DenseMatrix.java

package org.mklab.sdpj.datastructure;

import java.io.PrintStream;
import java.nio.channels.UnsupportedAddressTypeException;

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.gpack.atlas.ATLAS;
import org.mklab.sdpj.gpack.blaswrap.BLAS;


/**
 * Class <code>DenseMarix</code> is a data structure of dense matrix, 
 * which is saving the data of dense-matrix data type.  
 * <br>密行列を表すクラスです。
 * 
 * @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 DenseMatrix<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 rowSize;
  /** 列数 */
  private int columnSize;
  /** enum data type, if matrix is dense, value is DENSE; else DIAGONAL*/
  private DenseDiagonal denseOrDiagonal;
  /** vector for saving data of a dense matrix*/
  public RS[] denseElements;
  /** vector for saving data of a diagonal matrix*/
  public RS[] diagonalElements;
  
  /** */
  private RS unit;

  /**
   * 新しく生成された{@link DenseMatrix}オブジェクトを初期化します。
   * @param unit unit
   */
  private DenseMatrix(RS unit) {
    this.unit = unit;
    this.columnSize = 0;
    this.rowSize = 0;
    this.denseOrDiagonal = DenseDiagonal.DENSE;
  }

  /**
   * 新しく生成された<code>DenseMatrix</code>オブジェクトを初期化します。
   * 
   * @param rowSize 行数
   * @param columnSize 列数
   * @param denseOrDiagonal 密行列または対角行列
   * @param unit unit
   */
  public DenseMatrix(int rowSize, int columnSize, DenseDiagonal denseOrDiagonal, RS unit) {
    this.unit = unit;
    if (rowSize < 0 || columnSize < 0) {
      throw new RuntimeException("rDenseMatrix:: Dimensions are nonpositive"); //$NON-NLS-1$
    }

    //final RS unit = Tools.getUnitNumber();

    this.rowSize = rowSize;
    this.columnSize = columnSize;
    this.denseOrDiagonal = denseOrDiagonal;

    switch (denseOrDiagonal) {
      case DENSE:
        final int length = rowSize * columnSize;
        this.denseElements = unit.createArray(length);
        ATLAS.setValue(length, unit.createZero(), this.denseElements, 1);
        break;
      case DIAGONAL:
        if (rowSize != columnSize) {
          throw new RuntimeException("rDenseMatrix:: Diagonal must be Square matrix"); //$NON-NLS-1$
        }
        this.diagonalElements = unit.createArray(columnSize);
        ATLAS.setValue(columnSize, unit.createZero(), this.diagonalElements, 1);
        break;
      default:
      throw new UnsupportedOperationException();  
    }
  }

  /**
   * 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 getConditionNumber(PrintStream out, PrintStream min, PrintStream max) {
    //RS unit = getElementUnit();
    RS[][] elements = this.unit.createArray(this.rowSize, this.columnSize);

    switch (this.denseOrDiagonal) {
      case DENSE:
        for (int i = 0; i < this.rowSize; ++i) {
          for (int j = 0; j < this.columnSize; ++j) {
            elements[i][j] = this.denseElements[i + this.columnSize * j];
          }
        }
        break;
      case DIAGONAL:
        //RS zero = 
        for (int i = 0; i < this.rowSize; i++) {
          elements[i][i] = this.diagonalElements[i];
          for (int j = 0; j < this.columnSize; ++j) {
            if (i != j) {
              elements[i][j] = this.unit.createZero();
            }
          }
        }
        break;
      default:
        throw new UnsupportedOperationException();
    }

    //NumericalMatrix<RS> cond = new BaseNumericalMatrix<>(this.rowSize, this.columnSize, elements);
    RM cond = elements[0][0].createGrid(elements);
    out.printf("%s\n", cond.conditionNumber().toString("%4.3e")); //$NON-NLS-1$//$NON-NLS-2$
    min.printf("%s\n", cond.minSingularValue().toString("%4.3e")); //$NON-NLS-1$//$NON-NLS-2$
    max.printf("%s\n", cond.maxSingularValue().toString("%4.3e")); //$NON-NLS-1$//$NON-NLS-2$
  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {
    final StringBuffer buff = new StringBuffer();

    switch (this.denseOrDiagonal) {
      case DENSE:
        buff.append("{"); //$NON-NLS-1$
        for (int i = 0; i < this.rowSize - 1; ++i) {
          if (i == 0) {
            buff.append(" "); //$NON-NLS-1$
          } else {
            buff.append("  "); //$NON-NLS-1$
          }
          buff.append("{"); //$NON-NLS-1$
          for (int j = 0; j < this.columnSize - 1; ++j) {
            buff.append(this.denseElements[i + this.columnSize * j] + ","); //$NON-NLS-1$
          }
          buff.append(this.denseElements[i + this.columnSize * (getColumnSize() - 1)] + " },\n"); //$NON-NLS-1$
        }
        if (getRowSize() > 1) {
          buff.append("  {"); //$NON-NLS-1$
        }
        for (int j = 0; j < this.columnSize - 1; ++j) {
          buff.append(this.denseElements[(getRowSize() - 1) + this.columnSize * j] + ","); //$NON-NLS-1$
        }
        buff.append(this.denseElements[(getRowSize() - 1) + this.columnSize * (getColumnSize() - 1)] + " }"); //$NON-NLS-1$

        if (getRowSize() > 1) {
          buff.append("   }\n"); //$NON-NLS-1$
        } else {
          buff.append("\n"); //$NON-NLS-1$
        }
        break;
      case DIAGONAL:
        buff.append("{"); //$NON-NLS-1$
        for (int j = 0; j < this.columnSize - 1; ++j) {
          buff.append(this.diagonalElements[j] + ","); //$NON-NLS-1$
        }
        if (getColumnSize() > 0) {
          buff.append(this.diagonalElements[getColumnSize() - 1] + "}\n"); //$NON-NLS-1$
        }
        break;
      default:
        throw new UnsupportedOperationException();
    }
    return buff.toString();
  }

  /**
   * コピー元の行列をコピーします。
   * 
   * @param source コピー元
   */
  public void copyFrom(DenseMatrix<RS,RM,CS,CM> source) {
    if (this == source) {
      return;
    }

    //final RS unit = source.getElementUnit();

    switch (source.denseOrDiagonal) {
      case DENSE:
        this.denseOrDiagonal = DenseDiagonal.DENSE;
        this.rowSize = source.getRowSize();
        this.columnSize = source.getColumnSize();
        int length = this.rowSize * this.columnSize;
        this.denseElements = this.unit.createArray(length);
        this.diagonalElements = null;
        BLAS.dcopy(length, source.denseElements, 1, this.denseElements, 1);
        break;
      case DIAGONAL:
        this.denseOrDiagonal = DenseDiagonal.DIAGONAL;
        this.rowSize = source.getRowSize();
        this.columnSize = source.getColumnSize();
        this.denseElements = null;
        this.diagonalElements = this.unit.createArray(getColumnSize());
        BLAS.dcopy(getColumnSize(), source.diagonalElements, 1, this.diagonalElements, 1);
        break;
      default:
        throw new UnsupportedAddressTypeException();
    }
  }

  /**
   * ゼロ行列にします。
   */
  public void setZero() {
    //final RS unit = getElementUnit();

    switch (this.denseOrDiagonal) {
      case DENSE:
        ATLAS.setValue(getRowSize() * this.columnSize, this.unit.createZero(), this.denseElements, 1);
        break;
      case DIAGONAL:
        ATLAS.setValue(getColumnSize(), this.unit.createZero(), this.diagonalElements, 1);
        break;
      default:
        throw new UnsupportedOperationException();
    }
  }

  /**
   * 単位行列にします。
   * 
   * @param value スカラー
   */
  public void setIdentity(RS value) {
    //final RS unit = value.createUnit();

    if (getRowSize() != this.columnSize) {
      throw new RuntimeException("rSparseMatrix:: Identity matrix must be square matrix"); //$NON-NLS-1$
    }

    switch (this.denseOrDiagonal) {
      case DENSE:
        final int length = this.rowSize * this.columnSize;
        ATLAS.setValue(length, this.unit.createZero(), this.denseElements, 1);
        final int step = this.columnSize + 1;
        ATLAS.setValue(getColumnSize(), value, this.denseElements, step);
        break;
      case DIAGONAL:
        ATLAS.setValue(getColumnSize(), value, this.diagonalElements, 1);
        break;
      default:
        throw new UnsupportedOperationException();
    }
  }

  /**
   * 密行列であるか判定します。
   * 
   * @return 密行列ならばtrue
   */
  public boolean isDense() {
    return this.denseOrDiagonal == DenseDiagonal.DENSE;
  }

  /**
   * 対角行列であるか判定します。
   * 
   * @return 対角行列ならばtrue
   */
  public boolean isDiagonal() {
    return this.denseOrDiagonal == DenseDiagonal.DIAGONAL;
  }

  /**
   * 密行列また対角行列の種類を返します。
   * 
   * @return 密行列また対角行列の種類
   */
  public DenseDiagonal getDenseOrDiagonal() {
    return this.denseOrDiagonal;
  }

  /**
   * 密行列の成分を返します。
   * 
   * @param index 成分の指数
   * @return 密行列の成分を返します。
   */
  public RS getDenseElement(final int index) {
    return this.denseElements[index];
  }

  /**
   * 密行列の成分を設定します。
   * 
   * @param index 成分の指数
   * @param value 値
   */
  public void setDenseElement(final int index, final RS value) {
    this.denseElements[index] = value;
  }

  /**
   * 対角行列の成分を返します。
   * 
   * @param index 成分の指数
   * @return 対角行列の成分を返します。
   */
  public RS getDiagonalElement(final int index) {
    return this.diagonalElements[index];
  }

  /**
   * 成分の単位を返します。
   * 
   * @return　成分の単位
   */
  public RS getElementUnit() {
    if (isDense()) {
      return getDenseElement(0).createUnit();
    }
    return getDiagonalElement(0).createUnit();
  }

  /**
   * 行数を設定します。
   * 
   * @param rowSize 行数
   */
  public void setRowSize(int rowSize) {
    this.rowSize = rowSize;
  }

  /**
   * 行数を返します。
   * 
   * @return 行数
   */
  public int getRowSize() {
    return this.rowSize;
  }

  /**
   * 列数を設定します。
   * 
   * @param columnSize 列数
   */
  public void setColunSize(int columnSize) {
    this.columnSize = columnSize;
  }

  /**
   * 列数を返します。
   * 
   * @return 列数
   */
  public int getColumnSize() {
    return this.columnSize;
  }

  /**
   * Returns clone.
   * 
   * @return clone
   */
  public DenseMatrix<RS,RM,CS,CM> createClone() {
    final DenseMatrix<RS,RM,CS,CM> inst = new DenseMatrix<>(this.unit);
    inst.rowSize = this.rowSize;
    inst.columnSize = this.columnSize;
    inst.denseOrDiagonal = this.denseOrDiagonal == null ? null : this.denseOrDiagonal;

    if (this.denseElements != null) {
      inst.denseElements = this.denseElements[0].createArray(this.denseElements.length);
      System.arraycopy(this.denseElements, 0, inst.denseElements, 0, this.denseElements.length);
    } else {
      inst.denseElements = null;
    }

    if (this.diagonalElements != null) {
      // inst.denseElements = this.diagonalElements[0].createArray(this.diagonalElements.length);
      inst.diagonalElements = this.diagonalElements[0].createArray(this.diagonalElements.length);
      System.arraycopy(this.diagonalElements, 0, inst.diagonalElements, 0, this.diagonalElements.length);
    } else {
      inst.diagonalElements = null;
    }
    return inst;
  }

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

    switch (getDenseOrDiagonal()) {
      case DENSE:
        for (int i = 0; i < this.rowSize; ++i) {
          for (int j = 0; j < this.columnSize; ++j) {
            data[i][j] = this.denseElements[i + this.columnSize * j];
          }
        }
        break;
      case DIAGONAL:
        //final RS zero = 
        for (int i = 0; i < this.rowSize; i++) {
          data[i][i] = this.diagonalElements[i];
          for (int j = 0; j < this.columnSize; j++) {
            if (i != j) {
              data[i][j] = this.unit.createZero();
            }
          }
        }
        break;
      default:
        throw new UnsupportedOperationException();
    }

    return data;
  }

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

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