CPD Results

The following document contains the results of PMD's CPD 6.13.0.

Duplications

File Line
org/mklab/sdpj/algebra/Algebra.java 294
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 416
    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 <T> 成分の型
   * @param a 元になる行列
   * @return ret 結果retMatと同じ
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> getCholesky(DenseMatrix<T> a) {
    final T unit = a.getElementUnit();

    DenseMatrix<T> 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());

        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 <T> 成分の型
   * @param a A
   * @return retMatのコピー
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> getInvLowTriangularMatrix(DenseMatrix<T> a) {
    final T unit = a.getElementUnit();

    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), a.getColumnSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @return コレスキー分解と逆行列
   */
  @SuppressWarnings("unchecked")
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T>[] getCholeskyAndInv(BlockDenseMatrix<T> a) {
    BlockDenseMatrix<T> choleskyMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());
    BlockDenseMatrix<T> inverseMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

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

  /**
   * @param <T> 成分の型
   * @param a A
   */
  private static <T extends NumericalScalar<T>> void getSymmetrize(DenseMatrix<T> a) {
    final T 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]

          T[] aMat_index2 = unit.createArray(a.denseElements.length - index2);
          T[] 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

          T 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]
          T[] A1 = unit.createArray(a.denseElements.length - index1);
          System.arraycopy(a.denseElements, index1, A1, 0, A1.length);
          T[] 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 <T> 成分の型
   * @param a A
   */
  public static <T extends NumericalScalar<T>> void getSymmetrize(BlockDenseMatrix<T> a) {
    for (int i = 0; i < a.getBlockSize(); ++i) {
      getSymmetrize(a.getBlock(i));
    }
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @return 転置行列
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> getTranspose(DenseMatrix<T> 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<T> 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 <T> 成分の型
   * @param a A
   * @return ブロック転置行列
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> getTranspose(BlockDenseMatrix<T> a) {
    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @return 計算が正常終了したならばtrue
   */
  public static <T extends NumericalScalar<T>> boolean choleskyFactorWithAdjust(DenseMatrix<T> a) {
    int rowSize = a.getRowSize();

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

    int info = Clapack.dpotrf("Lower", rowSize, a.denseElements, rowSize);
    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 <T> 成分の型
   * @param a A
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> solveSystems(DenseMatrix<T> a, Vector<T> 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<T> 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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(SparseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = b.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal());

    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];
          T value = a.getSparseElement(index).multiply(scalar);
          if (i != j) {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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);

            T[] bMat_i = unit.createArray(b.denseElements.length - i);
            System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
            T[] 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 {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(SparseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = b.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal());

    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];
          T value = a.getSparseElement(index);
          if (i != j) {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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);

            T[] bMat_i = unit.createArray(b.denseElements.length - i);
            System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
            T[] 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 {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, SparseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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];
          T value = b.getSparseElement(index).multiply(scalar);
          if (i != j) {
            T[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
            System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
            T[] 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);

            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            T[] 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 {
            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            T[] 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, SparseMatrix<T> b) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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];
          T value = b.getSparseElement(index);
          if (i != j) {
            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            T[] 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);

            T[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
            System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
            T[] 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 {
            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            T[] 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 <T> 成分の型
   * @param a A
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, T scalar) {
    DenseMatrix<T> 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 <T> 成分の型
   * @param a A
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, T scalar) {
    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a a
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> multiply(Vector<T> a, T scalar) {
    Vector<T> ans = a.createClone();
    BLAS.dscal(ans.getDimension(), scalar, ans.getElements(), 1);
    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a a
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockVector<T> multiply(BlockVector<T> a, T scalar) {
    BlockVector<T> ans = a.createClone();

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> multiply(DenseMatrix<T> a, Vector<T> b) {
    if (a.getColumnSize() != b.getDimension()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    Vector<T> ans = b.createClone();

    final T 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 <T> 成分の型
   * @param a A
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> multiply(DenseMatrix<T> a, Vector<T> b, T scalar) {
    if (a.getColumnSize() != b.getDimension()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    Vector<T> 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 <T> 成分の型
   * @param a A
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockVector<T> multiply(BlockDenseMatrix<T> a, BlockVector<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    BlockVector<T> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockVector<T> multiply(BlockDenseMatrix<T> a, BlockVector<T> b) {
    if (b.getBlockSize() != a.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockVector<T> ans = b.createClone();

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = a.createClone();

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockSparseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockSparseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = aMat**T * bMat
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> tran_multiply(DenseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> tran_multiply(DenseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<T> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal());
    final T unit = a.getElementUnit();

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> tran_multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> tran_multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = aMat * bMat**T
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply_tran(DenseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply_tran(DenseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getColumnSize() != b.getColumnSize()) {
      Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal());

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply_tran(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply_tran(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    T scalar = null;
    BlockDenseMatrix<T> ans = multiply_tran(a, b, scalar);
    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <T> 成分の型
   * @param a a
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> plus(Vector<T> a, Vector<T> b, T scalar) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

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

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a a
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> plus(Vector<T> a, Vector<T> b) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

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

  /**
   * ans = a - b
   * 
   * @param <T> 成分の型
   * @param a a
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> Vector<T> minus(Vector<T> a, Vector<T> b) {
    if (a.getDimension() != b.getDimension()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

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

  /**
   * ans = a + b*scalar
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(DenseMatrix<T> a, DenseMatrix<T> b, T 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<T> 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(DenseMatrix<T> a, DenseMatrix<T> 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<T> 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> minus(DenseMatrix<T> a, DenseMatrix<T> 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<T> 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 <T> 成分の型 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(SparseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ret = (*scalar) * b
    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(SparseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> minus(SparseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(DenseMatrix<T> a, SparseMatrix<T> b, T scalar) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ret = a
    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> plus(DenseMatrix<T> a, SparseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ans = a
    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> minus(DenseMatrix<T> a, SparseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
    }

    DenseMatrix<T> 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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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 <T> 成分の型
   * @param a a
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockVector<T> plus(BlockVector<T> a, BlockVector<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    T unit = a.getElementUnit();
    BlockVector<T> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

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

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a a
   * @param b b
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockVector<T> plus(BlockVector<T> a, BlockVector<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    T unit = a.getElementUnit();
    BlockVector<T> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);

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

    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

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

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> minus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

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

    return ans;
  }

  /**
   * ans = a + b*scalar
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> minus(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs());

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

    return ans;
  }

  /**
   * ans + a * b*scalar
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockSparseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockSparseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> minus(BlockDenseMatrix<T> a, BlockSparseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> Vector<T> let(Vector<T> a, final char operator, T scalar) {
    switch (operator) {
      case '*':
        Vector<T> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockVector<T> let(BlockVector<T> a, final char operator, T scalar) {
    switch (operator) {
      case '*':
        BlockVector<T> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '*' (*scalar)
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockDenseMatrix<T> a, final char operator, T scalar) {
    switch (operator) {
      case '*':
        BlockDenseMatrix<T> ans = multiply(a, scalar);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = a '+' '-' b*(*scalar)
   * 
   * @param <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param b b
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> Vector<T> let(Vector<T> a, final char operator, Vector<T> b, T scalar) {
    Vector<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  public static <T extends NumericalScalar<T>> Vector<T> let(Vector<T> a, final char operator, Vector<T> b) {
    Vector<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(DenseMatrix<T> a, final char operator, DenseMatrix<T> b, T scalar) {
    DenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(DenseMatrix<T> a, final char operator, DenseMatrix<T> b) {
    DenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(SparseMatrix<T> a, final char operator, DenseMatrix<T> b, T scalar) {
    DenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(SparseMatrix<T> a, final char operator, DenseMatrix<T> b) {
    DenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(DenseMatrix<T> a, final char operator, SparseMatrix<T> b, T scalar) {
    DenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> DenseMatrix<T> let(DenseMatrix<T> a, final char operator, SparseMatrix<T> b) {
    DenseMatrix<T> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockDenseMatrix<T> a, final char operator, BlockDenseMatrix<T> b, T scalar) {
    BlockDenseMatrix<T> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockDenseMatrix<T> a, final char operator, BlockDenseMatrix<T> b) {
    BlockDenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockSparseMatrix<T> a, final char operator, BlockDenseMatrix<T> b, T scalar) {
    BlockDenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockSparseMatrix<T> a, final char operator, BlockDenseMatrix<T> b) {
    BlockDenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockDenseMatrix<T> a, final char operator, BlockSparseMatrix<T> b, T scalar) {
    BlockDenseMatrix<T> ans;

    switch (operator) {
      case '+':
        ans = plus(a, b, scalar);
        return ans;
      case '-':
        final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> BlockDenseMatrix<T> let(BlockDenseMatrix<T> a, final char operator, BlockSparseMatrix<T> b) {
    BlockDenseMatrix<T> ans;
    switch (operator) {
      case '+':
        ans = plus(a, b);
        return ans;
      case '-':
        ans = minus(a, b);
        //final T unit = a.getElementUnit();
        //final T 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 <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  public static <T extends NumericalScalar<T>> Vector<T> let(DenseMatrix<T> a, final char operator, Vector<T> b) {
    Vector<T> 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 <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(Vector<T> a, final char operator, Vector<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(DenseMatrix<T> a, final char operator, DenseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(DenseMatrix<T> a, final char operator, SparseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(b, a);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(SparseMatrix<T> a, final char operator, DenseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a a
   * @param operator 演算子
   * @param b b
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(BlockVector<T> a, final char operator, BlockVector<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(BlockDenseMatrix<T> a, final char operator, BlockDenseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(BlockSparseMatrix<T> a, final char operator, BlockDenseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(a, b);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }

  /**
   * ans = inner_product(a,b) // op = '.'
   * 
   * @param <T> 成分の型
   * @param a A
   * @param operator 演算子
   * @param b B
   * @return 解
   */
  public static <T extends NumericalScalar<T>> T run(BlockDenseMatrix<T> a, final char operator, BlockSparseMatrix<T> b) {
    switch (operator) {
      case '.':
        T ans = getInnerProduct(b, a);
        return ans;
      default:
        throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
    }
  }
}
File Line
org/mklab/sdpj/main/SdpjMain.java 123
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 127
    buff.append("SdpjMain start at " + new java.util.Date() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$

    ParameterType parameterType = ParameterType.PARAMETER_DEFAULT;

    if (this.argc == 1) {
      showMessage(arguments[0]);
    }

    if (arguments[0].charAt(0) == '-') {
      for (int i = 0; i < this.argc; ++i) {
        String target = arguments[i];
        if (target.equals("-dd") && i + 1 < this.argc) { //$NON-NLS-1$
          this.dataFile = arguments[i + 1];
          i++;
          continue;
        }
        if (target.equals("-ds") && i + 1 < this.argc) { //$NON-NLS-1$
          this.dataFile = arguments[i + 1];
          i++;
          this.isDataSparse = true;
          continue;
        }
        if (target.equals("-id") && i + 1 < this.argc) { //$NON-NLS-1$
          this.initFile = arguments[i + 1];
          i++;
          this.isInitFile = true;
          continue;
        }
        if (target.equals("-is") && i + 1 < this.argc) { //$NON-NLS-1$
          this.initFile = arguments[i + 1];
          i++;
          this.isInitFile = true;
          this.isInitSparse = true;
          continue;
        }
        if (target.equals("-o") && i + 1 < this.argc) { //$NON-NLS-1$
          this.outFile = arguments[i + 1];
          i++;
          continue;
        }
        if (target.equals("-p") && i + 1 < this.argc) { //$NON-NLS-1$
          this.paraFile = arguments[i + 1];
          i++;
          this.isParameter = true;
          continue;
        }
        if (target.equals("-k") && i + 1 < this.argc) { //$NON-NLS-1$
          this.kappaString = arguments[i + 1];
          Tools.message("Kappa = " + this.kappaString); //$NON-NLS-1$
          i++;
          continue;
        }

        if (target.equals("-pt") && i + 1 < this.argc) { //$NON-NLS-1$
          int tmp = atoi(arguments[i + 1]);
          switch (tmp) {
            case 0:
              parameterType = ParameterType.PARAMETER_DEFAULT;
              break;
            case 1:
              parameterType = ParameterType.PARAMETER_AGGRESSIVE;
              break;
            case 2:
              parameterType = ParameterType.PARAMETER_STABLE;
              break;
            case 3:
              parameterType = ParameterType.PARAMETER_MP115_DEFAULT;
              break;
            case 4:
              parameterType = ParameterType.PARAMETER_MP115_STABLE;
              break;
            default:
              parameterType = ParameterType.PARAMETER_DEFAULT;
          }
          if (tmp == 3 || tmp == 4) {
            this.dataType = "mpfloat"; //$NON-NLS-1$
          } else {
            this.dataType = "double"; //$NON-NLS-1$
          }
          i++;
          this.paraFile = null;
          this.isParameter = false;
          continue;
        }
      }
    } else { // SDPA argument
      this.dataFile = arguments[0];
      int len = this.dataFile.length();
      if (this.dataFile.charAt(len - 1) == 's' && this.dataFile.charAt(len - 2) == '-') {
        this.isDataSparse = true;
      }

      this.outFile = arguments[1];

      for (int i = 2; i < arguments.length; ++i) {
        if (arguments[i].equals("-pt") && i + 1 < this.argc) { //$NON-NLS-1$
          int tmp = atoi(arguments[i + 1]);
          switch (tmp) {
            case 0:
              parameterType = ParameterType.PARAMETER_DEFAULT;
              break;
            case 1:
              parameterType = ParameterType.PARAMETER_AGGRESSIVE;
              break;
            case 2:
              parameterType = ParameterType.PARAMETER_STABLE;
              break;
            case 3:
              parameterType = ParameterType.PARAMETER_MP115_DEFAULT;
              break;
            case 4:
              parameterType = ParameterType.PARAMETER_MP115_STABLE;
              break;
            default:
              parameterType = ParameterType.PARAMETER_DEFAULT;
          }
          if (tmp == 3 || tmp == 4) {
            this.dataType = "mpfloat"; //$NON-NLS-1$
          } else {
            this.dataType = "double"; //$NON-NLS-1$
          }
          i++;
          this.paraFile = null;
          this.isParameter = false;
        } else {
          if (arguments[i].equals("-p") && i + 1 < this.argc) { //$NON-NLS-1$
            this.paraFile = arguments[i + 1];
            i++;
            this.isParameter = true;
          } else {
            this.initFile = arguments[i];
            this.isInitFile = true;
            len = this.initFile.length();
            if (this.initFile.charAt(len - 1) == 's' && this.initFile.charAt(len - 2) == '-') {
              this.isInitSparse = true;
            }
          }
        }
      }
    }

    if (this.dataType.equals("")) { //$NON-NLS-1$
      if (this.paraFile == null) {
        this.paraFile = "param/param.sdpj"; //$NON-NLS-1$
      }
      creatDataType(this.paraFile);
      this.isParameter = true;
    }

    if (this.dataFile == null || this.outFile == null) {
      showMessage(arguments[0]);
    }

    buff.append("data      is " + this.dataFile); //$NON-NLS-1$
    if (this.isDataSparse) {
      buff.append(" : sparse\n"); //$NON-NLS-1$
    } else {
      buff.append(" : dense\n"); //$NON-NLS-1$
    }
    if (this.paraFile != null) {
      buff.append("parameter is " + this.paraFile + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    }
    if (this.outFile != null) {
      buff.append("out       is " + this.outFile + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
    }
    if (this.initFile != null) {
      buff.append("initial   is " + this.initFile); //$NON-NLS-1$
    }
    if (this.isInitFile) {
      if (this.isInitSparse) {
        buff.append(" : sparse\n"); //$NON-NLS-1$
      } else {
        buff.append(" : dense\n"); //$NON-NLS-1$
      }
    } else {
      buff.append("\n"); //$NON-NLS-1$
    }
    if (this.paraFile == null) {
      if (parameterType == ParameterType.PARAMETER_DEFAULT) {
        buff.append("set       is Double DEFAULT\n");// << endl; //$NON-NLS-1$
      } else if (parameterType == ParameterType.PARAMETER_AGGRESSIVE) {
        buff.append("set       is Double AGGRESSIVE\n");// << endl; //$NON-NLS-1$
      } else if (parameterType == ParameterType.PARAMETER_STABLE) {
        buff.append("set       is Double STABLE\n");// << endl; //$NON-NLS-1$
      } else if (parameterType == ParameterType.PARAMETER_MP115_DEFAULT) {
        buff.append("set       is MPFloat DEFAULT\n");// << endl; //$NON-NLS-1$
      } else if (parameterType == ParameterType.PARAMETER_MP115_STABLE) {
        buff.append("set       is MPFloat STABLE\n");// << endl; //$NON-NLS-1$
      }
    }

    PrintStream fpOut = new PrintStream(new FileOutputStream(new File(this.outFile)));

    if (this.dataType.equals("mpfloat")) { //$NON-NLS-1$
      MPFloat.setDefaultPrecisionDigits(this.precision);
      this.solver = new Solver<>(new MPFloat(1));
File Line
org/mklab/sdpj/algebra/Algebra.java 131
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 176
    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return aMat bMatの内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(DenseMatrix<T> a, DenseMatrix<T> 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.<T> 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(SparseMatrix<T> a, DenseMatrix<T> b) {
    if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final T unit = b.getElementUnit();
    int index = 0;
    int counter = 0;
    T 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];
          T 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];
          T value1 = a.getSparseElement(index);
          T 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];
          T value2 = a.getSparseElement(index + 1);
          T 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];
          T value3 = a.getSparseElement(index + 2);
          T 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];
          T value4 = a.getSparseElement(index + 3);
          T 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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    T ans = unit.create(0);
File Line
org/mklab/sdpj/algebra/Algebra.java 1736
org/mklab/sdpj/algebra/Algebra.java 1845
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1858
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1967
        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];
          T 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];
          T 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];
          T 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];
          T 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];
          T 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) {
File Line
org/mklab/sdpj/main/SdpjMain.java 328
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 334
      System.out.println("SdpjMain approximate significant digits : 16"); //$NON-NLS-1$
      if (this.kappaString != null) {
        this.solver.setKappa(this.<DoubleNumber> atof(this.kappaString));
      }
    } else {
      System.err.println("Please decide a data type correctly in your parameter file!"); //$NON-NLS-1$
      System.exit(0);
    }

    //set some parameters of solver before running 
    this.solver.setDataFile(this.dataFile);
    this.solver.setInitFile(this.initFile);
    this.solver.setOutFile(this.outFile);
    this.solver.setParaFile(this.paraFile);
    this.solver.setIsInitFile(this.isInitFile);
    this.solver.setIsInitSparse(this.isInitSparse);
    this.solver.setIsDataSparse(this.isDataSparse);
    this.solver.setIsParameter(this.isParameter);
    this.solver.setParameterType(parameterType);
    this.solver.setPrintStreamForDisplay(display);
    this.solver.setPrintStreamForSave(fpOut);

    Tools.message(buff.toString(), display);
    Tools.message(buff.toString(), fpOut);

    final double startTime = System.nanoTime();
    this.solver.run();
    final double endTime = System.nanoTime();
    System.out.printf("The elapsed time of solving this SDP problem is : %.3f seconds\n\n", new Double((endTime - startTime) / 1e9)); //$NON-NLS-1$
    this.dataType = ""; //$NON-NLS-1$
  }

  /**
   * the purpose of this method is show how to use <code>SdpjMain</code> 
   * when users make a mistake in using <code>SdpjMain</code>.
   * 
   * @param argv0 a string, which is "<code>sdpj</code>"
   */
  private void showMessage(String argv0) {
    StringBuffer buff = new StringBuffer();
    buff.append("\n"); //$NON-NLS-1$
    buff.append("************* Please assign data file and output file.*************\n"); //$NON-NLS-1$
    buff.append("\n"); //$NON-NLS-1$
    buff.append("------------------- option type 1 --------------------\n"); //$NON-NLS-1$
    buff.append(argv0 + " DataFile OutputFile [InitialPtFile]" + " [-pt parameters] [-p InitialParamFile]\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("parameters = 0 default, 1 aggressive," + " 2 stable, 3 MPFloat default, 4 MPFlaot stable\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example1-1: " + argv0 + " example.dat example.result\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example1-2: " + argv0 + " example.dat-s example.result\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example1-3: " + argv0 + " example.dat example.result example.ini\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example1-4: " + argv0 + " example.dat example.result -pt 2\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example1-5: " + argv0 + " example.dat example.result -p param.sdpj\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("\n"); //$NON-NLS-1$

    buff.append("------------------- option type 2 --------------------\n"); //$NON-NLS-1$
    buff.append(argv0 + " [option filename]+ \n"); //$NON-NLS-1$
    buff.append("  -dd : data dense :: -ds : data sparse     \n"); //$NON-NLS-1$
    buff.append("  -id : init dense :: -is : init sparse     \n"); //$NON-NLS-1$
    buff.append("  -o  : output     :: -p  : parameter       \n"); //$NON-NLS-1$
    buff.append("  -pt : parameters , 0 default, 1 aggressive, 2 stable,\n"); //$NON-NLS-1$
    buff.append("                   3 MPFloat default, 4 MPFlaot stable\n"); //$NON-NLS-1$
    buff.append("example2-1: " + argv0 + " -o example.result -dd example.dat\n"); //$NON-NLS-1$ //$NON-NLS-2$
    buff.append("example2-2: " + argv0 + " -ds example.dat-s -o example.result " + "-p param.sdpj\n"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
    buff.append("example2-3: " + argv0 + " -ds example.dat-s -o example.result " + "-pt 2\n"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$

    buff.append("\nNote:if you run program by command line mode,please refer to option type1\n"); //$NON-NLS-1$
    buff.append("or type2.Otherwise,you must give parameter args some information about \n"); //$NON-NLS-1$
    buff.append("DataFile & OutputFile,etc.If you do not input the InitialParamFile in Type1,\n"); //$NON-NLS-1$
    buff.append("the datatype will be used by default. And Type2 is also the same. Otherwise,\n"); //$NON-NLS-1$
    buff.append("the datatype is appointed based on the InitialParamFile. To run the program \n"); //$NON-NLS-1$
    buff.append("correctly,you must appoint datafile and outputfile.\n"); //$NON-NLS-1$
    buff.append("For example: args[0] = \"example.dat-s\", args[1] = \"example.result\"\n"); //$NON-NLS-1$
    buff.append("Recommedation:if you want to appoint a datatype, input a InitialParamFile.\n"); //$NON-NLS-1$
    buff.append("\n******************************************************************\n"); //$NON-NLS-1$

    System.out.println(new String(buff));
  }

  /**
   * The program is to change an inputting style of solver <code>sedumi</code> into the style of solver <code>sdpj</code>
   *  
   * @param arguments a command line string
   */
  private void sedumi(String[] arguments) {
    final Sedumi sedumi = new Sedumi();
    final String writeFile = arguments[1];
    final String fileFormat = writeFile.substring(writeFile.lastIndexOf(".")); //$NON-NLS-1$
    final String datsFile = writeFile.replaceAll(Pattern.quote(fileFormat), ".dat-s"); //$NON-NLS-1$ 

    System.out.println(datsFile + " is building..."); //$NON-NLS-1$
    sedumi.SedumiLoad(writeFile, datsFile);
    System.out.println("Done"); //$NON-NLS //$NON-NLS-1$
    System.out.println("please input command:% sdpj " + datsFile + " outputfile" + " [options]\n"); //$NON-NLS //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
    System.exit(0);
  }

  /**
   * This program is to change a <code>string</code> into a <code>integer</code>, which can only change a <code>char</code> of
   * the string. And if the char is not a <code>number</code> char, maybe it occurs an exception. 
   * <br>文字列の最初に数字きたらその整数を返します。 本来は文字列の先頭から走査していき、数字以外が来たらそこまでの整数を返すというもの。
   * 今回は引数の判別で、一桁のものしかないため、最初の数字しか考えていない。
   * 
   * @param str a integer string
   * @return a integer
   */
  private int atoi(String str) {

    switch (str.charAt(0)) {
      case '0':
        return 0;
      case '1':
        return 1;
      case '2':
        return 2;
      case '3':
        return 3;
      case '4':
        return 4;
      case '5':
        return 5;
      case '6':
        return 6;
      case '7':
        return 7;
      case '8':
        return 8;
      case '9':
        return 9;
      case '-':
        return -atoi(str.substring(1));
    }
    return 0;
  }

  /**
   * Change a data string into a data based on the data type, which is used in program. 
   * 
   * @param <E> a generic data type
   * @param string a data string
   * @return a data based on a data type
   */
  private <E extends NumericalScalar<E>> E atof(String string) {
    if (this.dataType.equals("double")) { //$NON-NLS-1$
      return (E)new DoubleNumber(1).valueOf(string);
    }
    if (this.dataType.equals("mpfloat")) { //$NON-NLS-1$
      return (E)new MPFloat(1).valueOf(string);
File Line
org/mklab/sdpj/algebra/Algebra.java 21
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 33
public class Algebra<E extends NumericalScalar<E>> {

  /**
   * 最小固有値の計算を行います 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 <T> 成分の型
   * @param a A
   * @param eigenValues 固有値
   * @param workVec 作業領域
   * @return 最小固有値
   */
  private static <T extends NumericalScalar<T>> T getMinEigenValue(DenseMatrix<T> a, Vector<T> eigenValues, Vector<T> workVec) {
    final T 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:
        T 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 <T> 成分の型
   * @param a A
   * @param eigenVec 固有ベクトル
   * @param workVec 作業領域
   * @return 最小固有値
   */
  static <T extends NumericalScalar<T>> T getMinEigenValue(BlockDenseMatrix<T> a, BlockVector<T> eigenVec, BlockVector<T> 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$
    }
    T 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$
      }
      T 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 <T> 成分の型
   * @param a a
   * @param b b
   * @return aVecとbVecの内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(Vector<T> a, Vector<T> b) {
    int size = a.getDimension();
    if (size != b.getDimension()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }
File Line
org/mklab/sdpj/algebra/Algebra.java 698
org/mklab/sdpj/algebra/Algebra.java 780
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 820
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 902
          T value = a.getSparseElement(index).multiply(scalar);
          if (i != j) {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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);

            T[] bMat_i = unit.createArray(b.denseElements.length - i);
            System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
            T[] 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 {
            T[] bMat_j = unit.createArray(b.denseElements.length - j);
            System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
            T[] 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$
File Line
org/mklab/sdpj/gpack/blaswrap/Dtrmm.java 10
org/mklab/sdpj/gpack/blaswrap/Dtrsm.java 11
public class Dtrmm<E extends NumericalScalar<E>> {

  /*  Purpose   
  =======   
  DTRMM  performs one of the matrix-matrix operations   
     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),   
  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or   
  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of   
     op( A ) = A   or   op( A ) = A'.   
  Parameters   
  ==========   
  SIDE   - CHARACTER*1.   
           On entry,  SIDE specifies whether  op( A ) multiplies B from   
           the left or right as follows:   
              SIDE = 'L' or 'l'   B := alpha*op( A )*B.   
              SIDE = 'R' or 'r'   B := alpha*B*op( A ).   
           Unchanged on exit.   
  UPLO   - CHARACTER*1.   
           On entry, UPLO specifies whether the matrix A is an upper or   
           lower triangular matrix as follows:   
              UPLO = 'U' or 'u'   A is an upper triangular matrix.   
              UPLO = 'L' or 'l'   A is a lower triangular matrix.   
           Unchanged on exit.   
  TRANSA - CHARACTER*1.   
           On entry, TRANSA specifies the form of op( A ) to be used in   
           the matrix multiplication as follows:   
              TRANSA = 'N' or 'n'   op( A ) = A.   
              TRANSA = 'T' or 't'   op( A ) = A'.   
              TRANSA = 'C' or 'c'   op( A ) = A'.   
           Unchanged on exit.   
  DIAG   - CHARACTER*1.   
           On entry, DIAG specifies whether or not A is unit triangular   
           as follows:   
              DIAG = 'U' or 'u'   A is assumed to be unit triangular.   
              DIAG = 'N' or 'n'   A is not assumed to be unit   
                                  triangular.   
           Unchanged on exit.   
  M      - INTEGER.   
           On entry, M specifies the number of rows of B. M must be at   
           least zero.   
           Unchanged on exit.   
  N      - INTEGER.   
           On entry, N specifies the number of columns of B.  N must be   
           at least zero.   
           Unchanged on exit.   
  ALPHA  - DOUBLE PRECISION.   
           On entry,  ALPHA specifies the scalar  alpha. When  alpha is   
           zero then  A is not referenced and  B need not be set before   
           entry.   
           Unchanged on exit.   
  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m   
           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.   
           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k   
           upper triangular part of the array  A must contain the upper   
           triangular matrix  and the strictly lower triangular part of   
           A is not referenced.   
           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k   
           lower triangular part of the array  A must contain the lower   
           triangular matrix  and the strictly upper triangular part of   
           A is not referenced.   
           Note that when  DIAG = 'U' or 'u',  the diagonal elements of   
           A  are not referenced either,  but are assumed to be  unity.   
           Unchanged on exit.   
  LDA    - INTEGER.   
           On entry, LDA specifies the first dimension of A as declared   
           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then   
           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'   
           then LDA must be at least max( 1, n ).   
           Unchanged on exit.   
  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).   
           Before entry,  the leading  m by n part of the array  B must   
           contain the matrix  B,  and  on exit  is overwritten  by the   
           transformed matrix.   
  LDB    - INTEGER.   
           On entry, LDB specifies the first dimension of B as declared   
           in  the  calling  (sub)  program.   LDB  must  be  at  least   
           max( 1, m ).   
           Unchanged on exit.   
  Level 3 Blas routine.   
  -- Written on 8-February-1989.   
     Jack Dongarra, Argonne National Laboratory.   
     Iain Duff, AERE Harwell.   
     Jeremy Du Croz, Numerical Algorithms Group Ltd.   
     Sven Hammarling, Numerical Algorithms Group Ltd.   
     Test the input parameters.   
     Parameter adjustments */

  /**
   * @param side side
   * @param uplo uplo
   * @param transa transa
   * @param diag diag
   * @param m m
   * @param n n
   * @param alpha alpha
   * @param a a
   * @param lda lda
   * @param b b
   * @param ldb ldb
   * @return result
   */
  public int execute(String side, String uplo, String transa, String diag, int m, int n, E alpha, E[] a, int lda, E[] b, int ldb) {
    final E unit = a[0].createUnit();

    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    int b_dim1 = ldb;
    int b_offset = 1 + b_dim1 * 1;
    int pointer_b = -b_offset;

    boolean lside = BLAS.lsame(side, "L"); //$NON-NLS-1$
    int nrowa;
    if (lside) {
      nrowa = m;
    } else {
      nrowa = n;
    }
    boolean nounit = BLAS.lsame(diag, "N"); //$NON-NLS-1$
    boolean upper = BLAS.lsame(uplo, "U"); //$NON-NLS-1$
    int info = 0;
    if (!lside && !BLAS.lsame(side, "R")) { //$NON-NLS-1$
      info = 1;
    } else if (!upper && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$
      info = 2;
    } else if (!BLAS.lsame(transa, "N") && !BLAS.lsame(transa, "T") && !BLAS.lsame(transa, "C")) { //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
      info = 3;
    } else if (!BLAS.lsame(diag, "U") && !BLAS.lsame(diag, "N")) { //$NON-NLS-1$ //$NON-NLS-2$
      info = 4;
    } else if (m < 0) {
      info = 5;
    } else if (n < 0) {
      info = 6;
    } else if (lda < Math.max(1, nrowa)) {
      info = 9;
    } else if (ldb < Math.max(1, m)) {
      info = 11;
    }
    if (info != 0) {
      BLAS.xerbla("DTRMM ", info); //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 171
org/mklab/sdpj/gpack/blaswrap/Dsyrk.java 145
    if (n == 0 || (alpha.isZero() || k == 0) && beta.equals(unit.create(1))) {
      return 0;
    }
    // And when  alpha.eq.zero.
    if (alpha.isZero()) {
      if (upper) {
        if (beta.isZero()) {
          for (int j = 1; j <= n; ++j) {
            for (int i = 1; i <= j; ++i) {
              c[j * c_dim1 + i + pointer_c] = unit.createZero();
            }
          }
        } else {
          for (int j = 1; j <= n; ++j) {
            for (int i = 1; i <= j; ++i) {
              int p = j * c_dim1 + i + pointer_c;
              c[p] = beta.multiply(c[p]);
            }
          }
        }
      } else {
        if (beta.isZero()) {
          for (int j = 1; j <= n; ++j) {
            for (int i = j; i <= n; ++i) {
              c[j * c_dim1 + i + pointer_c] = unit.createZero();
            }
          }
        } else {
          for (int j = 1; j <= n; ++j) {
            for (int i = j; i <= n; ++i) {
              int p = j * c_dim1 + i + pointer_c;
              c[p] = beta.multiply(c[p]);
            }
          }
        }
      }
      return 0;
    }

    if (BLAS.lsame(trans, "N")) { //$NON-NLS-1$
      // Form  C := alpha*A*B' + alpha*B*A' + C.
      if (upper) {
        for (int j = 1; j <= n; ++j) {
          if (beta.isZero()) {
            for (int i = 1; i <= j; ++i) {
              c[j * c_dim1 + i + pointer_c] = unit.createZero();
            }
          } else if (!beta.equals(unit.createUnit())) {
File Line
org/mklab/sdpj/algebra/Algebra.java 610
org/mklab/sdpj/algebra/Algebra.java 1341
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 732
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1463
            "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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, DenseMatrix<T> b) {
File Line
org/mklab/sdpj/algebra/Algebra.java 610
org/mklab/sdpj/algebra/Algebra.java 1341
org/mklab/sdpj/algebra/Algebra.java 1471
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 732
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1463
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1593
            "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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, DenseMatrix<T> b) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 350
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 362
    } else if (sname && LibF77.compare(c2, "OR", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nb = 32;
        }
      } else if (c3[0] == 'M') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nb = 32;
        }
      }
    } else if (cname && LibF77.compare(c2, "UN", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 482
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 494
    } else if (sname && LibF77.compare(c2, "OR", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nbmin = 2;
        }
      } else if (c3[0] == 'M') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nbmin = 2;
        }
      }
    } else if (cname && LibF77.compare(c2, "UN", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 650
org/mklab/sdpj/algebra/Algebra.java 1384
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 772
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1506
            "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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(SparseMatrix<T> a, DenseMatrix<T> b, T scalar) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dsteqr.java 771
org/mklab/sdpj/gpack/lapackwrap/Dsterf.java 474
  }

  /**
   * @param type type
   * @param kl kl
   * @param ku ku
   * @param cfrom cfrom
   * @param cto cto
   * @param m m
   * @param n2 n2
   * @param aIndex d
   * @param lda lda
   * @return result
   */
  private int dlascl_d(String type, int kl, int ku, E cfrom, E cto, int m, int n2, int aIndex, int lda) {
    final E unit = cfrom.createUnit();
    E[] aTemp = unit.createArray(this.d.length - aIndex);
    System.arraycopy(this.d, aIndex, aTemp, 0, aTemp.length);
    int info = Clapack.dlascl(type, kl, ku, cfrom, cto, m, n2, aTemp, lda);
    System.arraycopy(aTemp, 0, this.d, aIndex, aTemp.length);
    return info;
  }

  /**
   * @param type type
   * @param kl kl
   * @param ku ku
   * @param cfrom cfrom
   * @param cto cto
   * @param m m
   * @param n2 n2
   * @param aIndex e
   * @param lda lda
   * @return result
   */
  private int dlascl_e(String type, int kl, int ku, E cfrom, E cto, int m, int n2, int aIndex, int lda) {
    final E unit = cfrom.createUnit();
    E[] aTemp = unit.createArray(this.e.length - aIndex);
    System.arraycopy(this.e, aIndex, aTemp, 0, aTemp.length);
    int info = Clapack.dlascl(type, kl, ku, cfrom, cto, m, n2, aTemp, lda);
    System.arraycopy(aTemp, 0, this.e, aIndex, aTemp.length);
    return info;
  }
File Line
org/mklab/sdpj/algebra/Algebra.java 1384
org/mklab/sdpj/algebra/Algebra.java 1514
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1506
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1636
            "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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> tran_multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
File Line
org/mklab/sdpj/algebra/Algebra.java 650
org/mklab/sdpj/algebra/Algebra.java 1514
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 772
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1636
            "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 <T> 成分の型
   * @param a A
   * @param b B
   * @param scalar スカラー
   * @return 解
   */
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(SparseMatrix<T> a, DenseMatrix<T> b, T scalar) {
File Line
org/mklab/sdpj/algorithm/StepLength.java 128
org/mklab/sdpj/algorithm/StepLength.java 208
    this.dual = this.unit.create(9).divide(10);

    if (phase.getValue() == PhaseType.noINFO || phase.getValue() == PhaseType.dFEAS) {
      if (this.primal.isGreaterThan(1)) {// primal is infeasible
        this.primal = createUnit();
      }
    } else {// when primal is feasible,check stepP1 is effective or not.
      E incPrimalObj = Algebra.run(C, '.', newton.getDxMat());
      if (incPrimalObj.isGreaterThan(0)) {
        if (this.primal.isGreaterThan(this.dual)) {
          this.primal = this.dual;
        }
        if (this.primal.isGreaterThan(1)) {
          this.primal = createUnit();
        }
      }
    }
    if (phase.getValue() == PhaseType.noINFO || phase.getValue() == PhaseType.pFEAS) {
      if (this.dual.isGreaterThan(1)) {
        // dual is infeasible
        this.dual = createUnit();
      }
    } else {
      // when dual is feasible,check stepD1 is effective or not.
      E incDualObj = Algebra.run(b, '.', newton.getDyVec());
      if (incDualObj.isLessThan(0)) {
        if (this.dual.isGreaterThan(this.primal)) {
          this.dual = this.primal;
        }
        if (this.dual.isGreaterThan(1)) {
          this.dual = createUnit();
        }
      }
    }
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 311
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 568
          this.arraycopyWorkTempOffsetForWork();

          for (int j = 1; j <= k; ++j) {
            for (int i = 1; i <= m; ++i) {
              this.c[j * c_dim1 + i + pointer_c] = this.c[j * c_dim1 + i + pointer_c].subtract(this.work[j * work_dim1 + i + pointer_work]);
              /* L50: */
            }
            /* L60: */
          }
        }

      } else {
        /* Let  V =  ( V1 )   
                     ( V2 )    (last K rows)   
           where  V2  is unit upper triangular. */

        if (BLAS.lsame(side, "L")) { //$NON-NLS-1$

          /* Form  H * C  or  H' * C  where  C = ( C1 )   
                                                 ( C2 )   
          
             W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)   
             W := C2' */

          for (int j = 1; j <= k; ++j) {
            E[] cTemp = unit.createArray(n);
            E[] workTemp = unit.createArray(n);
            System.arraycopy(this.c, c_dim1 + m - k + j + pointer_c, cTemp, 0, n);
            System.arraycopy(this.work, j * work_dim1 + 1 + pointer_work, workTemp, 0, n);

            BLAS.dcopy(n, cTemp, ldc, workTemp, 1);

            System.arraycopy(cTemp, 0, this.c, c_dim1 + m - k + j + pointer_c, n);
            System.arraycopy(workTemp, 0, this.work, j * work_dim1 + 1 + pointer_work, n);
            /* L70: */
          }

          //配列をコピー work-->workTempOffset c-->cTempOffset
          this.arraycopyWorkForWorkTempOffset();
File Line
org/mklab/sdpj/convert/SDPAM.java 102
org/mklab/sdpj/convert/SDPAM.java 172
  public String getSDPADatSString() {
    final StringWriter sw = new StringWriter();
    final PrintWriter pw = new PrintWriter(new BufferedWriter(sw));

    pw.print(this.mDIM);
    pw.println("=mDIM"); //$NON-NLS-1$
    pw.print(this.nBLOCK);
    pw.println("=nBLOCK"); //$NON-NLS-1$

    if (this.bLOCKsTRUCT.size() > 1) {
      pw.print("("); //$NON-NLS-1$
      for (int i = 0; i < this.bLOCKsTRUCT.size(); i++) {
        pw.print(this.bLOCKsTRUCT.get(i));
        if (i == this.bLOCKsTRUCT.size() - 1) {
          pw.print(")"); //$NON-NLS-1$
        } else {
          pw.print(","); //$NON-NLS-1$
        }
      }
    } else {
      pw.print(this.bLOCKsTRUCT.get(0));
    }

    pw.println("=bLOCKsTRUCT"); //$NON-NLS-1$
    pw.print("{"); //$NON-NLS-1$

    for (int i = 0; i < this.c.length; i++) {
      pw.print(this.c[i]);
      if (i == this.c.length - 1) {
        pw.println("}"); //$NON-NLS-1$
      } else {
        pw.print(","); //$NON-NLS-1$
      }
    }

    int x = this.sdpamF[0].length;
File Line
org/mklab/sdpj/tool/Sedumi2.java 104
org/mklab/sdpj/tool/Sedumi2.java 202
    Solver<E> solver = new Solver<>((E)new MPFloat(1));

    solver.setProblem(problem);
    solver.setOutFile("out.txt"); //$NON-NLS-1$
    solver.setPrintStreamForDisplay(System.out);
    solver.run();
    //   solver.solve(problem, "out.txt", System.out, com); //$NON-NLS-1$
    Solution<E> s = solver.getSolution();

    MatxList ml = new MatxList();
    Vector<E> y = s.getYvec();
    ml.add(y.getElement(0).createGrid(y.getElements()));
    DenseMatrix<E> x = s.getXmat().getBlock(0);
    if (x.isDense()) {
      ml.add(x.denseElements[0].createGrid(x.denseElements));
    } else {
      ml.add(x.diagonalElements[0].createGrid(x.diagonalElements));
    }
    DenseMatrix<E> z = s.getZmat().getBlock(0);
    if (z.isDense()) {
      ml.add(z.denseElements[0].createGrid(z.denseElements));
    } else {
      ml.add(z.diagonalElements[0].createGrid(z.diagonalElements));
    }

    return ml;
  }
File Line
org/mklab/sdpj/gpack/blaswrap/Dtrmv.java 10
org/mklab/sdpj/gpack/blaswrap/Dtrsv.java 11
public class Dtrmv<E extends NumericalScalar<E>> {

  /*  Purpose   
  =======   
  DTRMV  performs one of the matrix-vector operations   
     x := A*x,   or   x := A'*x,   
  where x is an n element vector and  A is an n by n unit, or non-unit,   
  upper or lower triangular matrix.   
  Parameters   
  ==========   
  UPLO   - CHARACTER*1.   
           On entry, UPLO specifies whether the matrix is an upper or   
           lower triangular matrix as follows:   
              UPLO = 'U' or 'u'   A is an upper triangular matrix.   
              UPLO = 'L' or 'l'   A is a lower triangular matrix.   
           Unchanged on exit.   
  TRANS  - CHARACTER*1.   
           On entry, TRANS specifies the operation to be performed as   
           follows:   
              TRANS = 'N' or 'n'   x := A*x.   
              TRANS = 'T' or 't'   x := A'*x.   
              TRANS = 'C' or 'c'   x := A'*x.   
           Unchanged on exit.   
  DIAG   - CHARACTER*1.   
           On entry, DIAG specifies whether or not A is unit   
           triangular as follows:   
              DIAG = 'U' or 'u'   A is assumed to be unit triangular.   
              DIAG = 'N' or 'n'   A is not assumed to be unit   
                                  triangular.   
           Unchanged on exit.   
  N      - INTEGER.   
           On entry, N specifies the order of the matrix A.   
           N must be at least zero.   
           Unchanged on exit.   
  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).   
           Before entry with  UPLO = 'U' or 'u', the leading n by n   
           upper triangular part of the array A must contain the upper   
           triangular matrix and the strictly lower triangular part of   
           A is not referenced.   
           Before entry with UPLO = 'L' or 'l', the leading n by n   
           lower triangular part of the array A must contain the lower   
           triangular matrix and the strictly upper triangular part of   
           A is not referenced.   
           Note that when  DIAG = 'U' or 'u', the diagonal elements of   
           A are not referenced either, but are assumed to be unity.   
           Unchanged on exit.   
  LDA    - INTEGER.   
           On entry, LDA specifies the first dimension of A as declared   
           in the calling (sub) program. LDA must be at least   
           max( 1, n ).   
           Unchanged on exit.   
  X      - DOUBLE PRECISION array of dimension at least   
           ( 1 + ( n - 1 )*abs( INCX ) ).   
           Before entry, the incremented array X must contain the n   
           element vector x. On exit, X is overwritten with the   
           tranformed vector x.   
  INCX   - INTEGER.   
           On entry, INCX specifies the increment for the elements of   
           X. INCX must not be zero.   
           Unchanged on exit.   
  Level 2 Blas routine.   
  -- Written on 22-October-1986.   
     Jack Dongarra, Argonne National Lab.   
     Jeremy Du Croz, Nag Central Office.   
     Sven Hammarling, Nag Central Office.   
     Richard Hanson, Sandia National Labs.   
     Test the input parameters.   
     Parameter adjustments */

  /**
   * @param uplo uplo
   * @param trans trans
   * @param diag diag
   * @param n n
   * @param a a
   * @param lda lda
   * @param x x
   * @param incx incx
   * @return result
   */
  public int execute(String uplo, String trans, String diag, int n, E[] a, int lda, E[] x, int incx) {
    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    int info = 0;
    if (!BLAS.lsame(uplo, "U") && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$ //$NON-NLS-2$
      info = 1;
    } else if (!BLAS.lsame(trans, "N") && !BLAS.lsame(trans, "T") && !BLAS.lsame(trans, "C")) { //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
      info = 2;
    } else if (!BLAS.lsame(diag, "U") && !BLAS.lsame(diag, "N")) { //$NON-NLS-1$ //$NON-NLS-2$
      info = 3;
    } else if (n < 0) {
      info = 4;
    } else if (lda < Math.max(1, n)) {
      info = 6;
    } else if (incx == 0) {
      info = 8;
    }
    if (info != 0) {
      BLAS.xerbla("DTRMV ", info); //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Dsteqr.java 723
org/mklab/sdpj/gpack/lapackwrap/Dsterf.java 516
  }

  /**
   * @param id id
   * @param n n
   * @param dIndex dIndex
   * @return result
   */
  private int dlasrt(String id, int n, int dIndex) {
    final E unit = this.d[0].createUnit();
    E[] dTemp = unit.createArray(this.d.length - dIndex);
    System.arraycopy(this.d, dIndex, dTemp, 0, dTemp.length);
    int info = Clapack.dlasrt(id, n, dTemp);
    System.arraycopy(dTemp, 0, this.d, dIndex, dTemp.length);
    return info;
  }

  /**
   * @param norm input
   * @param n input
   * @param dIndex d input
   * @param eIndex e input
   * @return ans output
   */
  private E dlanst(String norm, int n, int dIndex, int eIndex) {
    final E unit = this.d[0].createUnit();
    E[] dTemp = unit.createArray(this.d.length - dIndex);
    E[] eTemp = unit.createArray(this.e.length - eIndex);
    System.arraycopy(this.d, dIndex, dTemp, 0, dTemp.length);
    System.arraycopy(this.e, eIndex, eTemp, 0, eTemp.length);
    return Clapack.dlanst(norm, n, dTemp, eTemp);
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasrt.java 164
org/mklab/sdpj/gpack/lapackwrap/Dlasrt.java 201
            } while (d[pointer_d + i].isGreaterThan(dmnmx));
            if (i < j) {
              E tmp = d[pointer_d + i];
              d[pointer_d + i] = d[pointer_d + j];
              d[pointer_d + j] = tmp;
            }
          } while (i < j);

          if (j - start > endd - j - 1) {
            ++stkpnt;
            stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
            stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
            ++stkpnt;
            stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
            stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
          } else {
            ++stkpnt;
            stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
            stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
            ++stkpnt;
            stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
            stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
          }
        } else {
File Line
org/mklab/sdpj/algebra/Algebra.java 877
org/mklab/sdpj/algebra/Algebra.java 960
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 999
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1082
            System.arraycopy(retMat_i, 0, ans.denseElements, ans.getRowSize() * i, retMat_i.length);
          } else {
            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            T[] 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$
File Line
org/mklab/sdpj/algorithm/StepLength.java 306
org/mklab/sdpj/algorithm/StepLength.java 350
    E[][] matrixdA = newton.getDzMat().getArrayOfElements();

    NumericalMatrix<E> A = new BaseNumericalMatrix<>(matrixA);
    NumericalMatrix<E> dA = new BaseNumericalMatrix<>(matrixdA);

    double alpha = 0.1;
    for (int i = 1; i <= 10; i++) {
      alpha *= i;
      dA = (NumericalMatrix<E>)dA.multiply(alpha);
      A = (NumericalMatrix<E>)A.add(dA);

      E minimumEigenValue;
      if (guarantee) {
        Guarantor guarantor = new Guarantor((VerifierFactory<MPFloat>)new EigenVerifierFactory<>((NumericalMatrix<E>)A.multiply(alpha)));
        guarantor.solveWithEffectiveDigits(26);
        IntervalMatrix<E> ans = (IntervalMatrix<E>)((EigenVerifiedSolution<E>)guarantor.getSolution()).getValue();
        minimumEigenValue = (E)ans.getInfimum().getRealPart().min();
      } else {
        NumericalMatrix<?> eigenValue = A.eigenValue();
File Line
org/mklab/sdpj/gpack/blaswrap/Dtrsm.java 272
org/mklab/sdpj/gpack/blaswrap/Dtrsm.java 301
            for (int k = 1; k <= j - 1; ++k) {
              int ap = j * a_dim1 + k + pointer_a;
              if (!a[ap].equals(unit.create(0), unit.getMachineEpsilon())) {
                for (int i = 1; i <= m; ++i) {
                  // b_ref(i, j) = b_ref(i, j) - a_ref(k, j) * b_ref(i,k);
                  int bp = j * b_dim1 + i + pointer_b;
                  b[bp] = b[bp].subtract(a[ap].multiply(b[k * b_dim1 + i + pointer_b]));
                }
              }
            }
            if (nounit) {
              E temp = unit.create(1).divide(a[j * a_dim1 + j + pointer_a]);

              for (int i = 1; i <= m; ++i) {
                // b_ref(i, j) = temp * b_ref(i, j);
                int p = j * b_dim1 + i + pointer_b;
                b[p] = temp.multiply(b[p]);
              }
            }
          }
        } else {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 375
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 626
          System.arraycopy(vTemp, 0, this.v, v_dim1 + m - k + 1 + pointer_v, vTemp.length);
          /* C2 := C2 - W' */

          //配列を元に戻す。
          this.arraycopyCTempOffsetForC();
          this.arraycopyWorkTempOffsetForWork();

          for (int j = 1; j <= k; ++j) {
            for (int i = 1; i <= n; ++i) {
              this.c[i * c_dim1 + m - k + j + pointer_c] = this.c[i * c_dim1 + m - k + j + pointer_c].subtract(this.work[j * work_dim1 + i + pointer_work]);
              /* L80: */
            }
            /* L90: */
          }

        } else if (BLAS.lsame(side, "R")) { //$NON-NLS-1$

          /* Form  C * H  or  C * H'  where  C = ( C1  C2 )   
             W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)   
             W := C2 */

          for (int j = 1; j <= k; ++j) {
            E[] cTemp = unit.createArray(m);
            E[] workTemp = unit.createArray(m);
            System.arraycopy(this.c, (n - k + j) * c_dim1 + 1 + pointer_c, cTemp, 0, m);
File Line
org/mklab/sdpj/main/SdpjMain.java 78
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 83
  public SdpjMain() {
    this.isInitFile = false;
    this.isInitSparse = false;
    this.isDataSparse = false;
    this.isParameter = false;
    this.dataFile = null;
    this.initFile = null;
    this.outFile = null;
    this.paraFile = null;
    this.kappaString = null;
    this.argc = 0;
    this.dataType = ""; //$NON-NLS-1$
    this.precision = 77;
  }

  /**
   * Returning a solver.
   * 
   * @return information of a solver of SdpjMain
   */
  public Solver<? extends NumericalScalar<?>> getSolver() {
    return this.solver;
  }

  /**
   * This is a solver of SdpjMain, which can compute a optimal solution of SDP problem, 
   * is an entry of solving a SDP problem.  
   * 
   * @param arguments a command line string
   * @throws IOException If I/O exception occurs
   */
  public void solver(String[] arguments) throws IOException {
    if (arguments.length < 2) {
      showMessage("sdpj"); //$NON-NLS-1$
      System.exit(0);
    }

    if (arguments[0].toLowerCase().equals("sedumi")) { //$NON-NLS-1$
      sedumi(arguments);
    }

    this.argc = arguments.length;
    PrintStream display = System.out;

    StringBuffer buff = new StringBuffer();
    buff.append("SdpjMain start at " + new java.util.Date() + "\n"); //$NON-NLS-1$ //$NON-NLS-2$
File Line
org/mklab/sdpj/main/SdpjMain.java 473
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 476
      return (E)new MPFloat(1).valueOf(string);
    }
    return null;
  }

  /**
   * Create a specific data type from a parameter file. Users can freely specify a data type, 
   * which can be used in <code>SdpjMain</code>.
   *
   * @param parameterFile Parameter file
   * @throws FileNotFoundException If no file found exception occurs
   * @throws IOException If I/O exception occurs
   */
  private void creatDataType(String parameterFile) throws FileNotFoundException, IOException {
    BufferedReader reader = new BufferedReader(new FileReader(new File(parameterFile)));

    while (this.dataType.equals("double") == false && this.dataType.equals("mpfloat") == false) { //$NON-NLS-1$//$NON-NLS-2$
      String string = reader.readLine();
      if (string == null) {
        break;
      }

      String type = string.toLowerCase();
      if (type.contains("double")) { //$NON-NLS-1$
        this.dataType = "double"; //$NON-NLS-1$
      } else if (type.contains("mpfloat")) { //$NON-NLS-1$
        this.dataType = "mpfloat"; //$NON-NLS-1$
      } else {
        this.dataType = type;
      }
    }

    if (this.dataType.equals("mpfloat")) { //$NON-NLS-1$
      this.precision = IO.readInteger(reader.readLine());
    }
  }
}
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 138
org/mklab/sdpj/gpack/blaswrap/Dsyrk.java 115
    int c_dim1 = ldc;
    int c_offset = 1 + c_dim1 * 1;
    int pointer_c = -c_offset;

    int nrowa;
    if (BLAS.lsame(trans, "N")) { //$NON-NLS-1$
      nrowa = n;
    } else {
      nrowa = k;
    }
    boolean upper = BLAS.lsame(uplo, "U"); //$NON-NLS-1$

    int info = 0;
    if (!upper && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$
      info = 1;
    } else if (!BLAS.lsame(trans, "N") && !BLAS.lsame(trans, "T") && !BLAS.lsame(trans, "C")) { //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
      info = 2;
    } else if (n < 0) {
      info = 3;
    } else if (k < 0) {
      info = 4;
    } else if (lda < Math.max(1, nrowa)) {
      info = 7;
    } else if (ldb < Math.max(1, nrowa)) {
File Line
org/mklab/sdpj/benchmark/Sdplib.java 332
org/mklab/sdpj/benchmark/Sdplib.java 446
      this.target = "0"; //$NON-NLS-1$
      String replacement = "100"; //$NON-NLS-1$
      for (Integer i = 0; i < 13; i++) {
        if (i == 0) {
          this.run(replacement);
          this.target = replacement;
        }
        if (i > 0 && i <= 4) {
          replacement = "124-" + i.toString(); //$NON-NLS-1$
          this.run(replacement);
          this.target = replacement;
        }
        if (i > 4 && i <= 8) {
          replacement = "250-" + new Integer(i - 4).toString(); //$NON-NLS-1$
          this.run(replacement);
          this.target = replacement;
        }
        if (i > 8 && i <= 12) {
          replacement = "500-" + new Integer(i - 8).toString(); //$NON-NLS-1$
          this.run(replacement);
          this.target = replacement;
        }
      }
    }
  }

  /**
   * test solver SDPJ with using problems of infp
   * 
   * @param name name of problem folder
   * @throws IOException 読み込めない場合
   */
  @SuppressWarnings("boxing")
  private void testInfp(String name) throws IOException {
File Line
org/mklab/sdpj/algebra/Algebra.java 2075
org/mklab/sdpj/algebra/Algebra.java 2183
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2197
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2305
          T 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];
          T value1 = b.getSparseElement(index).multiply(scalar);
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlatrd.java 322
org/mklab/sdpj/gpack/lapackwrap/Dlatrd.java 441
  private int dsymv(String uplo, int n, E alpha, int aIndex, int lda, int xIndex, int incx, E beta, int yIndex, int incy) {
    final E unit = alpha.createUnit();
    E[] aTemp = unit.createArray(this.a.length - aIndex);
    E[] x = unit.createArray(this.a.length - xIndex);
    E[] y = unit.createArray(this.w.length - yIndex);
    System.arraycopy(this.a, aIndex, aTemp, 0, aTemp.length);
    System.arraycopy(this.a, xIndex, x, 0, x.length);
    System.arraycopy(this.w, yIndex, y, 0, y.length);
File Line
org/mklab/sdpj/algebra/Algebra.java 682
org/mklab/sdpj/algebra/Algebra.java 763
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 804
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 885
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(SparseMatrix<T> a, DenseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = b.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal());

    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];
          T value = a.getSparseElement(index).multiply(scalar);
File Line
org/mklab/sdpj/algebra/Algebra.java 846
org/mklab/sdpj/algebra/Algebra.java 929
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 968
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1051
  private static <T extends NumericalScalar<T>> DenseMatrix<T> multiply(DenseMatrix<T> a, SparseMatrix<T> b, T scalar) {
    if (a.getColumnSize() != b.getRowSize()) {
      throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    DenseMatrix<T> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal());

    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];
          T value = b.getSparseElement(index).multiply(scalar);
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 360
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 492
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 555
        }
      }
    } else if (cname && LibF77.compare(c2, "UN", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 349
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 481
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 549
      }
    } else if (sname && LibF77.compare(c2, "OR", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algorithm/Newton.java 343
org/mklab/sdpj/algorithm/Newton.java 398
            if (i != j) {
              this.bMat.denseElements[i + this.bMat.getColumnSize() * j] = this.bMat.denseElements[i + this.bMat.getColumnSize() * j].add(value);
              this.bMat.denseElements[j + this.bMat.getColumnSize() * i] = this.bMat.denseElements[j + this.bMat.getColumnSize() * i].add(value);
            } else {
              this.bMat.denseElements[i + this.bMat.getColumnSize() * i] = this.bMat.denseElements[i + this.bMat.getColumnSize() * i].add(value);
            }
          }
        }
      } else {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 249
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 263
          for (int j = 1; j <= n - 1; ++j) {
            E ctemp = c[pointer_c + j];
            E stemp = s[pointer_s + j];
            if (!ctemp.isUnit() || !stemp.isZero()) {
              for (int i = 1; i <= m; ++i) {
                int p1 = (j) * a_dim1 + i + pointer_a;
                int p2 = (j + 1) * a_dim1 + i + pointer_a;
                E temp = a[p2];
                a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
                a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
              }
            }
          }
        } else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 865
org/mklab/sdpj/algebra/Algebra.java 955
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 987
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1077
            T[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
            System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
            T[] 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);
File Line
org/mklab/sdpj/algebra/Algebra.java 872
org/mklab/sdpj/algebra/Algebra.java 948
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 994
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1070
            T[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
            System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
            T[] 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);
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 186
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 200
          for (int j = 2; j <= m; ++j) {
            E ctemp = c[pointer_c + j - 1];
            E stemp = s[pointer_s + j - 1];
            if (!ctemp.isUnit() || !stemp.isZero()) {
              for (int i = 1; i <= n; ++i) {
                int p1 = i * a_dim1 + 1 + pointer_a;
                int p2 = i * a_dim1 + j + pointer_a;
                E temp = a[p2];
                a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
                a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
              }
            }
          }
        } else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 279
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 293
          for (int j = 2; j <= n; ++j) {
            E ctemp = c[pointer_c + j - 1];
            E stemp = s[pointer_s + j - 1];
            if (!ctemp.isUnit() || !stemp.isZero()) {
              for (int i = 1; i <= m; ++i) {
                int p1 = a_dim1 + i + pointer_a;
                int p2 = (j) * a_dim1 + i + pointer_a;
                E temp = a[p2];
                a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
                a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
              }
            }
          }
        } else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 2092
org/mklab/sdpj/algebra/Algebra.java 2200
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2214
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2322
          T 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];
          T value2 = b.getSparseElement(index + 1).multiply(scalar);
File Line
org/mklab/sdpj/algebra/Algebra.java 2104
org/mklab/sdpj/algebra/Algebra.java 2212
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2226
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2334
          T 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];
          T value3 = b.getSparseElement(index + 2).multiply(scalar);
File Line
org/mklab/sdpj/algebra/Algebra.java 2116
org/mklab/sdpj/algebra/Algebra.java 2224
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2238
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2346
          T 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];
          T value4 = b.getSparseElement(index + 3).multiply(scalar);
File Line
org/mklab/sdpj/algebra/Algebra.java 1745
org/mklab/sdpj/algebra/Algebra.java 1854
org/mklab/sdpj/algebra/Algebra.java 2183
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1867
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1976
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2305
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1963
org/mklab/sdpj/algebra/Algebra.java 2282
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2085
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2404
          T 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];
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 216
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 230
          for (int j = 1; j <= m - 1; ++j) {
            E ctemp = c[pointer_c + j];
            E stemp = s[pointer_s + j];
            if (!ctemp.isUnit() || !stemp.isZero()) {
              for (int i = 1; i <= n; ++i) {
                int p1 = i * a_dim1 + m + pointer_a;
                int p2 = i * a_dim1 + j + pointer_a;
                E temp = a[p2];
                a[p2] = (stemp.multiply(a[p1])).add(ctemp.multiply(temp));
                a[p1] = (ctemp.multiply(a[p1])).subtract(stemp.multiply(temp));
              }
            }
          }
        } else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 550
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 557
    } else if (sname && LibF77.compare(c2, "OR", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nx = 128;
        }
      }
    } else if (cname && LibF77.compare(c2, "UN", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/blaswrap/Dtrmm.java 284
org/mklab/sdpj/gpack/blaswrap/Dtrmm.java 306
            for (int j = 1; j <= k - 1; ++j) {
              if (!a[k * a_dim1 + j + pointer_a].isZero()) {
                E temp = alpha.multiply(a[k * a_dim1 + j + pointer_a]);
                for (int i = 1; i <= m; ++i) {
                  int p = j * b_dim1 + i + pointer_b;
                  b[p] = b[p].add(temp.multiply(b[k * b_dim1 + i + pointer_b]));
                }
              }
            }
            E temp = alpha;
            if (nounit) {
              temp = temp.multiply(a[k * a_dim1 + k + pointer_a]);
            }
            if (!temp.equals(unit.create(1))) {
File Line
org/mklab/sdpj/algebra/Algebra.java 1745
org/mklab/sdpj/algebra/Algebra.java 1854
org/mklab/sdpj/algebra/Algebra.java 2075
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1867
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1976
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2197
          T 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];
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2.java 148
org/mklab/sdpj/gpack/blaswrap/Dsyr2.java 179
            for (int i = 1; i <= j; ++i) {
              int p = j * a_dim1 + i + pointer_a;
              a[p] = a[p].add(x[i - 1].multiply(temp1)).add(y[i - 1].multiply(temp2));
            }
          }
        }
      } else {
        for (int j = 1; j <= n; ++j) {
          if (!x[jx - 1].isZero() || !y[jy - 1].isZero()) {
            E temp1 = alpha.multiply(y[jy - 1]);
            E temp2 = alpha.multiply(x[jx - 1]);
            int ix = kx;
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 234
org/mklab/sdpj/gpack/blaswrap/Dsyrk.java 203
                c[ci] = c[ci].add(a[ai].multiply(temp1)).add(b[bi].multiply(temp2));
              }
            }
          }
        }
      } else {
        for (int j = 1; j <= n; ++j) {
          if (beta.isZero()) {
            for (int i = j; i <= n; ++i) {
              c[j * c_dim1 + i + pointer_c] = unit.createZero();
            }
          } else if (!beta.equals(unit.createUnit())) {
            for (int i = j; i <= n; ++i) {
              int p = j * c_dim1 + i + pointer_c;
              c[p] = beta.multiply(c[p]);
            }
          }
          for (int l = 1; l <= k; ++l) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dorgtr.java 592
org/mklab/sdpj/gpack/lapackwrap/Dorgtr.java 611
    if (BLAS.lsame(side, "L")) { //$NON-NLS-1$
      /* Form  H * C */
      if (!tau.isZero()) {
        //array copy
        E[] cTemp = unit.createArray(c.length - (c_offset + pointer_c));
        E[] workTemp = unit.createArray(this.work.length - (1 + pointer_work));
        E[] vTemp = unit.createArray(v.length - (1 + pointer_v));
        System.arraycopy(c, c_offset + pointer_c, cTemp, 0, cTemp.length);
        System.arraycopy(this.work, 1 + pointer_work, workTemp, 0, workTemp.length);
        System.arraycopy(v, 1 + pointer_v, vTemp, 0, vTemp.length);

        /* w := C' * v */
        BLAS.dgemv("Transpose", m, n, c_b4, cTemp, ldc, vTemp, incv, c_b5, workTemp, 1); //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 1199
org/mklab/sdpj/algebra/Algebra.java 1285
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1321
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1407
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> multiply(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
File Line
org/mklab/sdpj/algebra/Algebra.java 2417
org/mklab/sdpj/algebra/Algebra.java 2555
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2539
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2677
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b, T scalar) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

    BlockDenseMatrix<T> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs());

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

    return ans;
  }

  /**
   * ans = a + b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 275
org/mklab/sdpj/gpack/lapackwrap/Dlarfb.java 290
          BLAS.dtrmm("Right", "Lower", "No transpose", "Unit", m, k, c_b14, this.vTempOffset, ldv, this.workTempOffset, ldwork); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$

          if (n > k) {
            /* W := W + C2 * V2 */
            E[] vTemp = unit.createArray(this.v.length - (v_dim1 + k + 1 + pointer_v));
            E[] cTemp = unit.createArray(this.c.length - ((k + 1) * c_dim1 + 1 + pointer_c));
            System.arraycopy(this.v, v_dim1 + k + 1 + pointer_v, this.v, 0, vTemp.length);
            System.arraycopy(this.c, (k + 1) * c_dim1 + 1 + pointer_c, cTemp, 0, cTemp.length);

            BLAS.dgemm("No transpose", "No transpose", m, k, n - k, c_b14, cTemp, ldc, vTemp, ldv, c_b14, this.workTempOffset, ldwork); //$NON-NLS-1$ //$NON-NLS-2$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 350
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 494
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 557
    } else if (sname && LibF77.compare(c2, "OR", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 362
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 482
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 550
    } else if (cname && LibF77.compare(c2, "UN", 2, 2) == 0) { //$NON-NLS-1$
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 2168
org/mklab/sdpj/algebra/Algebra.java 2270
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2290
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2392
      throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
    }

    // ans = a
    DenseMatrix<T> 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];
          T 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);
File Line
org/mklab/sdpj/algebra/Algebra.java 2440
org/mklab/sdpj/algebra/Algebra.java 2578
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2562
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2700
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> plus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
    }

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

    return ans;
  }

  /**
   * ans = a - b
   * 
   * @param <T> 成分の型
   * @param a A
   * @param b B
   * @return 解
   */
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> minus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyrk.java 235
org/mklab/sdpj/gpack/blaswrap/Dsyrk.java 251
          for (int i = 1; i <= j; ++i) {
            E temp = unit.createZero();
            for (int l = 1; l <= k; ++l) {
              int p = i * a_dim1 + l + pointer_a;
              temp = temp.add(a[p].multiply(a[j * a_dim1 + l + pointer_a]));
            }
            int p = j * c_dim1 + i + pointer_c;
            if (beta.isZero()) {
              c[p] = alpha.multiply(temp);
            } else {
              c[p] = (alpha.multiply(temp)).add(beta.multiply(c[p]));
            }
          }
        }
      } else {
File Line
org/mklab/sdpj/iocenter/SolverIO.java 336
org/mklab/sdpj/iocenter/SolverIO.java 342
            .printf(
                "%2d %s %s %s %s %s %s %s %s\n", this.pIteration, this.mu.getValue().toString("%4.2e"), this.theta.getDual().toString("%4.2e"), this.theta.getPrimal().toString("%4.2e"), this.solveInfo.getObjValDual().unaryMinus() //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
                    .toString("%+7.2e"), this.solveInfo.getObjValPrimal().unaryMinus().toString("%+7.2e"), this.alpha.getDual().toString("%4.2e"), this.alpha.getPrimal().toString("%4.2e"), this.beta.getValue().toString("%4.2e")); //$NON-NLS-1$//$NON-NLS-2$ //$NON-NLS-3$//$NON-NLS-4$//$NON-NLS-5$
      }
File Line
org/mklab/sdpj/main/SdpjMain.java 176
org/mklab/sdpj/main/SdpjMain.java 218
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 180
org/mklab/sdpj/main/SdpjMainWithEFFloat.java 222
        if (target.equals("-pt") && i + 1 < this.argc) { //$NON-NLS-1$
          int tmp = atoi(arguments[i + 1]);
          switch (tmp) {
            case 0:
              parameterType = ParameterType.PARAMETER_DEFAULT;
              break;
            case 1:
              parameterType = ParameterType.PARAMETER_AGGRESSIVE;
              break;
            case 2:
              parameterType = ParameterType.PARAMETER_STABLE;
              break;
            case 3:
              parameterType = ParameterType.PARAMETER_MP115_DEFAULT;
              break;
            case 4:
              parameterType = ParameterType.PARAMETER_MP115_STABLE;
              break;
            default:
              parameterType = ParameterType.PARAMETER_DEFAULT;
          }
          if (tmp == 3 || tmp == 4) {
            this.dataType = "mpfloat"; //$NON-NLS-1$
          } else {
            this.dataType = "double"; //$NON-NLS-1$
          }
          i++;
          this.paraFile = null;
          this.isParameter = false;
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 355
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 367
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 487
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 499
        }
      } else if (c3[0] == 'M') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 271
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 288
          for (int i = 1; i <= j; ++i) {
            E temp1 = unit.createZero();
            E temp2 = unit.createZero();
            for (int l = 1; l <= k; ++l) {
              temp1 = temp1.add(a[i * a_dim1 + l + pointer_a].multiply(b[j * b_dim1 + l + pointer_b]));
              temp2 = temp2.add(b[i * b_dim1 + l + pointer_b].multiply(a[j * a_dim1 + l + pointer_a]));
            }
            int p = j * c_dim1 + i + pointer_c;
            if (beta.isZero()) {
              c[p] = (alpha.multiply(temp1)).add(alpha.multiply(temp2));
File Line
org/mklab/sdpj/algebra/Algebra.java 1798
org/mklab/sdpj/algebra/Algebra.java 2128
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1920
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2250
          T 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) {
File Line
org/mklab/sdpj/datastructure/BlockDenseMatrix.java 106
org/mklab/sdpj/datastructure/BlockDenseMatrix.java 427
      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();
      }
    }
File Line
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 219
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 246
            for (int i = 1; i <= j; ++i) {
              int p = j * c_dim1 + i + pointer_c;
              c[p] = beta.multiply(c[p]);
            }
          }
          for (int l = 1; l <= k; ++l) {
            int aj = l * a_dim1 + j + pointer_a;
            int bj = l * b_dim1 + j + pointer_b;
            if (!a[aj].isZero() || !b[bj].isZero()) {
              E temp1 = alpha.multiply(b[bj]);
              E temp2 = alpha.multiply(a[aj]);
              for (int i = 1; i <= j; ++i) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 351
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 356
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 363
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 368
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nb = 32;
        }
      } else if (c3[0] == 'M') {
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 483
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 488
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 495
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 500
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
          nbmin = 2;
        }
      } else if (c3[0] == 'M') {
File Line
org/mklab/sdpj/algebra/Algebra.java 1798
org/mklab/sdpj/algebra/Algebra.java 1907
org/mklab/sdpj/algebra/Algebra.java 2236
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1920
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2029
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2358
          T 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:
File Line
org/mklab/sdpj/algebra/Algebra.java 2016
org/mklab/sdpj/algebra/Algebra.java 2335
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2138
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2457
          T 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:
File Line
org/mklab/sdpj/gpack/blaswrap/Dcopy.java 21
org/mklab/sdpj/gpack/blaswrap/Dcopy.java 51
  void dcopy(int size, double[] x, int incx, double[] y, int incy) {
    if (incx == 1 && incy == 1) {
      System.arraycopy(x, 0, y, 0, size);
      return;
    }

    int ix = 0;
    int iy = 0;
    if (incx < 0) {
      ix = (-(size) + 1) * incx;
    }
    if (incy < 0) {
      iy = (-(size) + 1) * incy;
    }
    for (int i = 1; i <= size; ++i) {
      y[iy] = x[ix];
      ix += incx;
      iy += incy;
    }
  }
File Line
org/mklab/sdpj/datastructure/DenseMatrix.java 127
org/mklab/sdpj/datastructure/SparseMatrix.java 130
        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) {
File Line
org/mklab/sdpj/algebra/Algebra.java 1774
org/mklab/sdpj/algebra/Algebra.java 1883
org/mklab/sdpj/algebra/Algebra.java 2212
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1896
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2005
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2334
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1786
org/mklab/sdpj/algebra/Algebra.java 1895
org/mklab/sdpj/algebra/Algebra.java 2224
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1908
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2017
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2346
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1992
org/mklab/sdpj/algebra/Algebra.java 2311
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2114
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2433
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 2004
org/mklab/sdpj/algebra/Algebra.java 2323
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2126
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2445
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 272
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 331
    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();

    T ans = unit.create(0);
File Line
org/mklab/sdpj/algebra/Algebra.java 1907
org/mklab/sdpj/algebra/Algebra.java 2128
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2029
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2250
          T 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) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 351
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 363
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 488
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 500
      if (c3[0] == 'G') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 356
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 368
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 483
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 495
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 551
org/mklab/sdpj/gpack/lapackwrap/Ilaenv.java 558
      } else if (c3[0] == 'M') {
        if (LibF77.compare(c4, "QR", 2, 2) == 0 || LibF77.compare(c4, "RQ", 2, 2) == 0 || LibF77.compare(c4, "LQ", 2, 2) == 0 || LibF77.compare(c4, "QL", 2, 2) == 0 || LibF77.compare(c4, "HR", 2, 2) == 0 || LibF77.compare(c4, "TR", 2, 2) == 0 //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
            || LibF77.compare(c4, "BR", 2, 2) == 0) { //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 728
org/mklab/sdpj/algebra/Algebra.java 810
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 850
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 932
            "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]));
File Line
org/mklab/sdpj/algebra/Algebra.java 894
org/mklab/sdpj/algebra/Algebra.java 976
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1016
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1098
            "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)));
File Line
org/mklab/sdpj/algebra/Algebra.java 1762
org/mklab/sdpj/algebra/Algebra.java 1871
org/mklab/sdpj/algebra/Algebra.java 2200
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1884
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1993
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2322
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1980
org/mklab/sdpj/algebra/Algebra.java 2299
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2102
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2421
          T 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];
File Line
org/mklab/sdpj/problem/ProblemCreater.java 225
org/mklab/sdpj/problem/ProblemCreaterOnMemory.java 161
            if (i <= j && !tmp.isZero()) {
              SparseMatrix<E> target = this.C.getBlock(k);
              int count = target.nonZeroCount;
              target.rowIndex[count] = i;
              target.columnIndex[count] = j;
              // C must be opposite sign.
              target.sparseElements[count] = tmp.unaryMinus();
              target.nonZeroCount++;
              if (i == j) {
                target.nonZeroEffect++;
              } else {
                target.nonZeroEffect += 2;
              }
            }
          }
        }
      } else {
        for (int j = 0; j < -size; ++j) {
          E tmp = unit.valueOf(F[0][counter]);
File Line
org/mklab/sdpj/algebra/Algebra.java 2128
org/mklab/sdpj/algebra/Algebra.java 2236
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2250
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2358
          T 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:
File Line
org/mklab/sdpj/problem/ProblemCreater.java 260
org/mklab/sdpj/problem/ProblemCreaterOnMemory.java 195
              if (i <= j && !tmp.isZero()) {
                SparseMatrix<E> target = this.A[k].getBlock(p);
                int count = target.nonZeroCount;
                target.rowIndex[count] = i;
                target.columnIndex[count] = j;
                target.sparseElements[count] = tmp;
                target.nonZeroCount++;
                if (i == j) {
                  target.nonZeroEffect++;
                } else {
                  target.nonZeroEffect += 2;
                }
              }
            }
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            E tmp = unit.valueOf(F[k + 1][counter]);
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlatrd.java 334
org/mklab/sdpj/gpack/lapackwrap/Dsytd2.java 326
  }

  /**
   * @param n n
   * @param alpha a[?]
   * @param xIndex a
   * @param incx indx
   * @param tau tau
   * @return result
   */
  private E[] dlarfg(int n, E alpha, int xIndex, int incx, E tau) {
    final E unit = alpha.createUnit();
    E[] x = unit.createArray(this.a.length - xIndex);
    System.arraycopy(this.a, xIndex, x, 0, x.length);
    E[] ans = Clapack.dlarfg(n, alpha, x, incx, tau);
    System.arraycopy(x, 0, this.a, xIndex, x.length);
    return ans;
  }
File Line
org/mklab/sdpj/algebra/Algebra.java 264
org/mklab/sdpj/algebra/Algebra.java 286
  private static <T extends NumericalScalar<T>> T getInnerProduct(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
    }

    final T unit = a.getElementUnit();
    T ans = unit.create(0);

    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 <T> 成分の型
   * @param a A
   * @param b B
   * @return 内積
   */
  private static <T extends NumericalScalar<T>> T getInnerProduct(BlockSparseMatrix<T> a, BlockDenseMatrix<T> b) {
File Line
org/mklab/sdpj/algebra/Algebra.java 1762
org/mklab/sdpj/algebra/Algebra.java 1871
org/mklab/sdpj/algebra/Algebra.java 2092
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1884
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1993
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2214
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1774
org/mklab/sdpj/algebra/Algebra.java 1883
org/mklab/sdpj/algebra/Algebra.java 2104
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1896
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2005
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2226
          T 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];
File Line
org/mklab/sdpj/algebra/Algebra.java 1786
org/mklab/sdpj/algebra/Algebra.java 1895
org/mklab/sdpj/algebra/Algebra.java 2116
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1908
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2017
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2238
          T 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];
File Line
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 159
org/mklab/sdpj/gpack/lapackwrap/Dlasr.java 173
            if (!ctemp.isUnit() || stemp.isZero()) {
              for (int i = 1; i <= n; ++i) {
                int p1 = i * a_dim1 + j + pointer_a;
                int p2 = i * a_dim1 + j + 1 + pointer_a;
                E temp = a[p2];
                a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
                a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
              }
            }
          }
        } else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
File Line
org/mklab/sdpj/gpack/lapackwrap/Dpotf2.java 82
org/mklab/sdpj/gpack/lapackwrap/Dpotrf.java 79
    E c_b12 = unit.create(1);

    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    /* Function Body */
    int info = 0;
    boolean upper = BLAS.lsame(uplo, "U"); //$NON-NLS-1$
    if (!upper && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$
      info = -1;
    } else if (n < 0) {
      info = -2;
    } else if (lda < Math.max(1, n)) {
      info = -4;
    }

    if (info != 0) {
      BLAS.xerbla("DPOTF2", -info); //$NON-NLS-1$
File Line
org/mklab/sdpj/algebra/Algebra.java 1737
org/mklab/sdpj/algebra/Algebra.java 1846
org/mklab/sdpj/algebra/Algebra.java 1955
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1859
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1968
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2077
          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];
          T 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);
File Line
org/mklab/sdpj/algebra/Algebra.java 2462
org/mklab/sdpj/algebra/Algebra.java 2601
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2584
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2723
  private static <T extends NumericalScalar<T>> BlockDenseMatrix<T> minus(BlockDenseMatrix<T> a, BlockDenseMatrix<T> b) {
    if (a.getBlockSize() != b.getBlockSize()) {
      throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
    }

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

    return ans;
  }
File Line
org/mklab/sdpj/iocenter/SolverIO.java 351
org/mklab/sdpj/iocenter/SolverIO.java 357
                "%2d  %s  %s  %s  %s  %s  %s  %s  %s\n", this.pIteration, this.mu.getValue().toString(), this.theta.getPrimal().toString(), this.theta.getDual().toString(), this.solveInfo.getObjValPrimal().toString(), //$NON-NLS-1$
                this.solveInfo.getObjValDual().toString(), this.alpha.getPrimal().toString(), this.alpha.getDual().toString(), this.beta.getValue().toString());
      }
File Line
org/mklab/sdpj/algebra/Algebra.java 1753
org/mklab/sdpj/algebra/Algebra.java 1862
org/mklab/sdpj/algebra/Algebra.java 1971
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1875
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1984
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2093
            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];
          T 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);
File Line
org/mklab/sdpj/algebra/Algebra.java 2191
org/mklab/sdpj/algebra/Algebra.java 2290
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2313
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 2412
            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];
          T 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);
File Line
org/mklab/sdpj/gpack/blaswrap/Dgemm.java 123
org/mklab/sdpj/gpack/blaswrap/Dsyr2k.java 127
  int dgemm(String transa, String transb, int m, int n, int k, E alpha, E[] a, int lda, E[] b, int ldb, E beta, E[] c, int ldc) {
    final E unit = a[0].createUnit();

    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    int b_dim1 = ldb;
    int b_offset = 1 + b_dim1 * 1;
    int pointer_b = -b_offset;

    int c_dim1 = ldc;
    int c_offset = 1 + c_dim1 * 1;
    int pointer_c = -c_offset;
File Line
org/mklab/sdpj/algebra/Algebra.java 806
org/mklab/sdpj/algebra/Algebra.java 972
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 928
org/mklab/sdpj/algebra/AlgebraWithEFFloat.java 1094
        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) {
File Line
org/mklab/sdpj/gpack/lapackwrap/Dorgtr.java 283
org/mklab/sdpj/gpack/lapackwrap/Dorgtr.java 724
    if (m < 0) {
      info = -1;
    } else if (n < 0 || n > m) {
      info = -2;
    } else if (k < 0 || k > n) {
      info = -3;
    } else if (lda < Math.max(1, m)) {
      info = -5;
    } else if (lwork < Math.max(1, n) && !lquery) {
      info = -8;
    }
    if (info != 0) {
      BLAS.xerbla("DORGQL", -info); //$NON-NLS-1$