Algebra.java
package org.mklab.sdpj.algebra;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.BlockVector;
import org.mklab.sdpj.datastructure.DenseMatrix;
import org.mklab.sdpj.datastructure.SparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.gpack.blaswrap.BLAS;
import org.mklab.sdpj.gpack.lapackwrap.Clapack;
import org.mklab.sdpj.tool.Tools;
/**
* TODO elements[n]で渡してる所はそのポインタってだけで、実はそれ以降の配列も渡してたりする？要確認 行列演算を行うクラスです。
*
* @author takafumi
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
*/
public class Algebra<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {
/**
* 最小固有値の計算を行います aMat is rewritten. aMat must be symmetric. eigenVec is the space of eigen values and needs memory of length aMat.nRow workVec is temporary space and needs 3*aMat.nRow-1 length
* memory.
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param eigenValues 固有値
* @param workVec 作業領域
* @return 最小固有値
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getMinEigenValue(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> eigenValues, Vector<RS,RM,CS,CM> workVec) {
final RS unit = a.getElementUnit();
int n = a.getRowSize();
switch (a.getDenseOrDiagonal()) {
case DENSE:
int lwork = 3 * n - 1;
// "N" means that we need not eigen vectors
// "L" means that we refer only lower triangular.
int info = Clapack.dsyev("NonVectors", "Lower", n, a.denseElements, n, eigenValues.getElements(), workVec.getElements(), lwork); //$NON-NLS-1$ //$NON-NLS-2$
if (info != 0) {
return unit.createZero();
}
return eigenValues.getElements()[0];
case DIAGONAL:
RS minimumEigenValue = a.diagonalElements[0];
eigenValues.getElements()[0] = a.diagonalElements[0];
for (int i = 1; i < n; ++i) {
eigenValues.getElements()[i] = a.diagonalElements[i];
if (a.diagonalElements[i].isLessThan(minimumEigenValue)) {
minimumEigenValue = a.diagonalElements[i];
}
}
return minimumEigenValue;
default:
throw new IllegalArgumentException();
}
}
/**
* 最小固有値の計算を行います。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param eigenVec 固有ベクトル
* @param workVec 作業領域
* @return 最小固有値
*/
static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getMinEigenValue(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> eigenVec, BlockVector<RS,RM,CS,CM> workVec) {
int n = a.getBlock(0).getRowSize();
if (eigenVec.getBlock(0).getDimension() != n) {
throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
}
if (workVec.getBlock(0).getDimension() != 3 * n - 1) {
throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
}
RS min_eigen = getMinEigenValue(a.getBlock(0), eigenVec.getBlock(0), workVec.getBlock(0));
for (int l = 1; l < a.getBlockSize(); ++l) {
n = a.getBlock(l).getRowSize();
if (eigenVec.getBlock(l).getDimension() != n) {
throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
}
if (workVec.getBlock(l).getDimension() != 3 * n - 1) {
throw new RuntimeException("getMinEigenValue:: different memory size"); //$NON-NLS-1$
}
RS tmp_eigen = getMinEigenValue(a.getBlock(l), eigenVec.getBlock(l), workVec.getBlock(l));
if (tmp_eigen.isLessThan(min_eigen)) {
min_eigen = tmp_eigen;
}
}
return min_eigen;
}
/**
* 内積の計算を行います。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @return aVecとbVecの内積
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
int size = a.getDimension();
if (size != b.getDimension()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
return BLAS.<RS,RM,CS,CM> ddot(size, a.getElements(), 1, b.getElements(), 1);
}
/**
* 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @return 最小固有値
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
// 内積＝Al.getInnerProduct(vec1,vec2)；に書き換えたほうがいいかもね！
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
RS ret = unit.create(0);
for (int i = 0; i < a.getBlockSize(); ++i) {
ret = ret.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
}
return ret;
}
/**
* 内積の計算を行います。 AI.getInnerProduct(ret,aVec,bVec)をAI.ret=getInnerProduct(aVec,bVec)に変更。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return aMat bMatの内積
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
switch (a.getDenseOrDiagonal()) {
case DENSE:
int length = a.getRowSize() * a.getColumnSize();
return BLAS.<RS,RM,CS,CM> ddot(length, a.denseElements, 1, b.denseElements, 1);
case DIAGONAL:
return BLAS.ddot(a.getColumnSize(), a.diagonalElements, 1, b.diagonalElements, 1);
default:
throw new IllegalArgumentException();
}
}
/**
* 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 内積
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
final RS unit = b.getElementUnit();
int index = 0;
int counter = 0;
RS ans;
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
// Attension: in SPARSE case, only half elements
// are stored. And bMat must be DENSE case.
ans = unit.create(0);
int amari = a.nonZeroCount % 4;
int shou = a.nonZeroCount / 4;
for (index = 0; index < amari; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index);
if (i == j) {
ans = ans.add(value.multiply(b.denseElements[i + b.getRowSize() * j]));
} else {
ans = ans.add(value.multiply((b.denseElements[i + b.getRowSize() * j].add(b.denseElements[j + b.getRowSize() * i]))));
}
}
index = amari;
counter = 0;
for (counter = 0; counter < shou; ++counter) {
int i1 = a.rowIndex[index];
int j1 = a.columnIndex[index];
RS value1 = a.getSparseElement(index);
RS ans1;
if (i1 == j1) {
ans1 = value1.multiply(b.denseElements[i1 + b.getRowSize() * j1]);
} else {
ans1 = value1.multiply((b.denseElements[i1 + b.getRowSize() * j1].add(b.denseElements[j1 + b.getRowSize() * i1])));
}
int i2 = a.rowIndex[index + 1];
int j2 = a.columnIndex[index + 1];
RS value2 = a.getSparseElement(index + 1);
RS ans2;
if (i2 == j2) {
ans2 = value2.multiply(b.denseElements[i2 + b.getRowSize() * j2]);
} else {
ans2 = value2.multiply((b.denseElements[i2 + b.getRowSize() * j2].add(b.denseElements[j2 + b.getRowSize() * i2])));
}
int i3 = a.rowIndex[index + 2];
int j3 = a.columnIndex[index + 2];
RS value3 = a.getSparseElement(index + 2);
RS ans3;
if (i3 == j3) {
ans3 = value3.multiply(b.denseElements[i3 + b.getRowSize() * j3]);
} else {
ans3 = value3.multiply((b.denseElements[i3 + b.getRowSize() * j3].add(b.denseElements[j3 + b.getRowSize() * i3])));
}
int i4 = a.rowIndex[index + 3];
int j4 = a.columnIndex[index + 3];
RS value4 = a.getSparseElement(index + 3);
RS ans4;
if (i4 == j4) {
ans4 = value4.multiply(b.denseElements[i4 + b.getRowSize() * j4]);
} else {
ans4 = value4.multiply((b.denseElements[i4 + b.getRowSize() * j4].add(b.denseElements[j4 + b.getRowSize() * i4])));
}
index += 4;
ans = ans.add(ans1.add(ans2).add(ans3).add(ans4));
//20090120修正したけど 不具合でたらここもチェックで。
}
return ans;
case DENSE:
int length = a.getRowSize() * a.getColumnSize();
return BLAS.ddot(length, a.denseElements, 1, b.denseElements, 1);
case DIAGONAL:
return BLAS.ddot(a.getColumnSize(), a.diagonalElements, 1, b.diagonalElements, 1);
default:
throw new IllegalArgumentException();
}
}
/**
* 内積の計算を行います Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 内積
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
RS ans = unit.create(0);
for (int i = 0; i < a.getBlockSize(); ++i) {
ans = ans.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
}
return ans;
}
/**
* 内積の計算を行います。 Al.getInnerProduct(ret,aVec,bVec)をret=Al.getInnerProduct(aVec,bVec)に変更。
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 内積
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS getInnerProduct(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("getInnerProduct:: different memory size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
RS ans = unit.create(0);
for (int i = 0; i < a.getBlockSize(); ++i) {
ans = ans.add(getInnerProduct(a.getBlock(i), b.getBlock(i)));
}
return ans;
}
/**
* getCholesky(retMat,aMat)をret=getCholesky(aMat)に変更
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a 元になる行列
* @return ret 結果retMatと同じ
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> getCholesky(DenseMatrix<RS,RM,CS,CM> a) {
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans;
switch (a.getDenseOrDiagonal()) {
case DENSE:
ans = a.createClone();
int info = Clapack.dpotrf("Lower", ans.getRowSize(), ans.denseElements, ans.getRowSize()); //$NON-NLS-1$
if (info != 0) {
throw new RuntimeException("cannot cholesky decomposition\n" + "Could you try with smaller gammaStar?"); //$NON-NLS-1$ //$NON-NLS-2$
}
// Make matrix as lower triangular matrix
for (int j = 0; j < ans.getColumnSize(); ++j) {
int shou = j / 4;
int amari = j % 4;
for (int i = 0; i < amari; ++i) {
ans.denseElements[i + ans.getColumnSize() * j] = unit.create(0);
}
int i;
int count;
for (i = amari, count = 0; count < shou; ++count, i += 4) {
ans.denseElements[i + ans.getColumnSize() * j] = unit.create(0);
ans.denseElements[i + 1 + ans.getColumnSize() * j] = unit.create(0);
ans.denseElements[i + 2 + ans.getColumnSize() * j] = unit.create(0);
ans.denseElements[i + 3 + ans.getColumnSize() * j] = unit.create(0);
}
}
return ans;
case DIAGONAL:
ans = new DenseMatrix<>(a.getRowSize(), a.getColumnSize(), a.getDenseOrDiagonal(), unit);
int shou = a.getColumnSize() / 4;
int amari = a.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.diagonalElements[j].sqrt();
}
int j;
int count;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.diagonalElements[j].sqrt();
ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].sqrt();
ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].sqrt();
ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].sqrt();
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @return retMatのコピー
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> getInvLowTriangularMatrix(DenseMatrix<RS,RM,CS,CM> a) {
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), a.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
ans.setIdentity(unit.create(1));
int nRow = a.getRowSize();
int nCol = a.getColumnSize();
int retnRow = ans.getRowSize();
BLAS.dtrsm("Left", "Lower", "NoTraspose", "NonUnitDiagonal", nRow, nCol, unit.create(1), a.denseElements, nRow, ans.denseElements, retnRow); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
return ans;
case DIAGONAL:
int shou = a.getColumnSize() / 4;
int amari = a.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = unit.create(1).divide(a.diagonalElements[j]);
}
for (int j = amari, counter = 0; counter < shou; ++counter, j += 4) {
ans.diagonalElements[j] = unit.create(1).divide(a.diagonalElements[j]);// 1.0
ans.diagonalElements[j + 1] = unit.create(1).divide(a.diagonalElements[j + 1]);
ans.diagonalElements[j + 2] = unit.create(1).divide(a.diagonalElements[j + 2]);
ans.diagonalElements[j + 3] = unit.create(1).divide(a.diagonalElements[j + 3]);
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @return コレスキー分解と逆行列
*/
@SuppressWarnings("unchecked")
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM>[] getCholeskyAndInv(BlockDenseMatrix<RS,RM,CS,CM> a) {
BlockDenseMatrix<RS,RM,CS,CM> choleskyMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
BlockDenseMatrix<RS,RM,CS,CM> inverseMat = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> ans1 = getCholesky(a.getBlock(i));
choleskyMat.setBlock(i, ans1);
DenseMatrix<RS,RM,CS,CM> ans2 = getInvLowTriangularMatrix(ans1);
inverseMat.setBlock(i, ans2);
}
return new BlockDenseMatrix[] {choleskyMat, inverseMat};
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> void getSymmetrize(DenseMatrix<RS,RM,CS,CM> a) {
final RS unit = a.getElementUnit();
//int index = 0;
switch (a.getDenseOrDiagonal()) {
case DENSE:
if (a.getRowSize() != a.getColumnSize()) {
throw new RuntimeException("getSymmetrize:: different memory size"); //$NON-NLS-1$
}
for (int index = 0; index < a.getRowSize() - 1; ++index) {
int index1 = index + index * a.getRowSize() + 1;
int index2 = index + (index + 1) * a.getRowSize();
int length = a.getRowSize() - 1 - index;
// aMat.de_ele[index1] += aMat.de_ele[index2]
RS[] aMat_index2 = unit.createArray(a.denseElements.length - index2);
RS[] aMat_index1 = unit.createArray(a.denseElements.length - index1);
System.arraycopy(a.denseElements, index1, aMat_index1, 0, aMat_index1.length);
System.arraycopy(a.denseElements, index2, aMat_index2, 0, aMat_index2.length);
//BLAS.daxpy(length, unit.create(1), aMat_index2, a.getRowSize(), aMat_index1, 1);
BLAS.dxpy(length, aMat_index2, a.getRowSize(), aMat_index1, 1);
System.arraycopy(aMat_index1, 0, a.denseElements, index1, aMat_index1.length);
// aMat.de_ele[index1] /= 2.0
RS half = unit.create(5).multiply(unit.create(10).power(-1));
System.arraycopy(a.denseElements, index1, aMat_index1, 0, aMat_index1.length);
BLAS.dscal(length, half, aMat_index1, 1);
System.arraycopy(aMat_index1, 0, a.denseElements, index1, aMat_index1.length);
// aMat.de_ele[index2] = aMat.de_ele[index1]
RS[] A1 = unit.createArray(a.denseElements.length - index1);
System.arraycopy(a.denseElements, index1, A1, 0, A1.length);
RS[] A2 = unit.createArray(a.denseElements.length - index2);
System.arraycopy(a.denseElements, index2, A2, 0, A2.length);
BLAS.dcopy(length, A1, 1, A2, a.getRowSize());
System.arraycopy(A2, 0, a.denseElements, index2, A2.length);
}
break;
case DIAGONAL:
// Nothing needs.
break;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> void getSymmetrize(BlockDenseMatrix<RS,RM,CS,CM> a) {
for (int i = 0; i < a.getBlockSize(); ++i) {
getSymmetrize(a.getBlock(i));
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @return 転置行列
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> getTranspose(DenseMatrix<RS,RM,CS,CM> a) {
if (a.getRowSize() != a.getColumnSize()) {
throw new RuntimeException("getTranspose:: different memory size"); //$NON-NLS-1$
// Of course, a non-symmetric matrix has its transposed matrix,
// but in this algorithm we have to make transposed matrix only when symmetric matrix.
}
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
switch (a.getDenseOrDiagonal()) {
case DENSE:
for (int i = 0; i < a.getRowSize(); ++i) {
int shou = (i + 1) / 4;
int amari = (i + 1) / 4;
for (int j = 0; j < amari; ++j) {
int index1 = i + a.getColumnSize() * j;
int index2 = j + a.getColumnSize() * i;
ans.denseElements[index1] = a.denseElements[index2];
ans.denseElements[index2] = a.denseElements[index1];
}
int counter, j;
for (j = amari, counter = 0; counter < shou; ++counter, j += 4) {
int index1 = i + a.getColumnSize() * j;
int index_1 = j + a.getColumnSize() * i;
ans.denseElements[index1] = a.denseElements[index_1];
ans.denseElements[index_1] = a.denseElements[index1];
int index2 = i + a.getColumnSize() * (j + 1);
int index_2 = (j + 1) + a.getColumnSize() * i;
ans.denseElements[index2] = a.denseElements[index_2];
ans.denseElements[index_2] = a.denseElements[index2];
int index3 = i + a.getColumnSize() * (j + 2);
int index_3 = (j + 2) + a.getColumnSize() * i;
ans.denseElements[index3] = a.denseElements[index_3];
ans.denseElements[index_3] = a.denseElements[index3];
int index4 = i + a.getColumnSize() * (j + 3);
int index_4 = (j + 3) + a.getColumnSize() * i;
ans.denseElements[index4] = a.denseElements[index_4];
ans.denseElements[index_4] = a.denseElements[index4];
}
}
return ans;
case DIAGONAL:
throw new IllegalArgumentException();
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @return ブロック転置行列
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> getTranspose(BlockDenseMatrix<RS,RM,CS,CM> a) {
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = getTranspose(a.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @return 計算が正常終了したならばtrue
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> boolean choleskyFactorWithAdjust(DenseMatrix<RS,RM,CS,CM> a) {
int rowSize = a.getRowSize();
//int info = new Dpotrf<RS>().rATL_dpotrfL(rowSize, a.denseElements, rowSize);
//a.setRowSize(rowSize);
int info = Clapack.dpotrf("Lower", rowSize, a.denseElements, rowSize); //$NON-NLS-1$
a.setRowSize(rowSize);
if (info < 0) {
System.err.println("cholesky argument is wrong " + (-info)); //$NON-NLS-1$
return false;
} else if (info > 0) {
System.err.println("cholesky miss condition :: not positive definite" + " :: info = " + info); //$NON-NLS-1$ //$NON-NLS-2$
return false;
}
return true;
}
/**
* solve aMat * xVec = bVec aMat must be Cholesky Factorized. aMat must have done Cholesky factorized.
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> solveSystems(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getDimension() || a.getRowSize() != a.getColumnSize()) {
throw new RuntimeException("solveSystems:: different memory size"); //$NON-NLS-1$
}
if (a.isDense() == false) {
throw new RuntimeException("solveSystems:: matrix type must be DENSE"); //$NON-NLS-1$
}
Vector<RS,RM,CS,CM> x = b.createClone();
BLAS.dtrsv("Lower", "NoTranspose", "NonUnit", a.getRowSize(), a.denseElements, a.getColumnSize(), x.getElements(), 1); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
BLAS.dtrsv("Lower", "Transpose", "NonUnit", a.getRowSize(), a.denseElements, a.getColumnSize(), x.getElements(), 1); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
return x;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getColumnSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getColumnSize() != b.getRowSize()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* Dense=Sparse*Dense
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getColumnSize() != b.getRowSize()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = b.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal(), unit);
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
for (int index = 0; index < a.nonZeroCount; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index).multiply(scalar);
if (i != j) {
RS[] bMat_j = unit.createArray(b.denseElements.length - j);
System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
RS[] retMat_i = unit.createArray(ans.denseElements.length - i);
System.arraycopy(ans.denseElements, i, retMat_i, 0, retMat_i.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_i, ans.getRowSize());
System.arraycopy(retMat_i, 0, ans.denseElements, i, retMat_i.length);
RS[] bMat_i = unit.createArray(b.denseElements.length - i);
System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_i, b.getRowSize(), retMat_j, ans.getRowSize());
System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
} else {
RS[] bMat_j = unit.createArray(b.denseElements.length - j);
System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_j, ans.getRowSize());
System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || b.isDiagonal() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = scalar.multiply(a.getDiagonalElement(j).multiply(b.diagonalElements[j]));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = scalar.multiply(a.getDiagonalElement(j)).multiply(b.diagonalElements[j]);
ans.diagonalElements[j + 1] = scalar.multiply(a.getDiagonalElement(j + 1).multiply(b.diagonalElements[j + 1]));
ans.diagonalElements[j + 2] = scalar.multiply(a.getDiagonalElement(j + 2).multiply(b.diagonalElements[j + 2]));
ans.diagonalElements[j + 3] = scalar.multiply(a.getDiagonalElement(j + 3).multiply(b.diagonalElements[j + 3]));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* TODO TODO Dense=Sparse*Dense
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getColumnSize() != b.getRowSize()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = b.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), b.getDenseOrDiagonal(), unit);
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
for (int index = 0; index < a.nonZeroCount; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index);
if (i != j) {
RS[] bMat_j = unit.createArray(b.denseElements.length - j);
System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
RS[] retMat_i = unit.createArray(ans.denseElements.length - i);
System.arraycopy(ans.denseElements, i, retMat_i, 0, retMat_i.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_i, ans.getRowSize());
System.arraycopy(retMat_i, 0, ans.denseElements, i, retMat_i.length);
RS[] bMat_i = unit.createArray(b.denseElements.length - i);
System.arraycopy(b.denseElements, i, bMat_i, 0, bMat_i.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_i, b.getRowSize(), retMat_j, ans.getRowSize());
System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
} else {
RS[] bMat_j = unit.createArray(b.denseElements.length - j);
System.arraycopy(b.denseElements, j, bMat_j, 0, bMat_j.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - j);
System.arraycopy(ans.denseElements, j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, bMat_j, b.getRowSize(), retMat_j, ans.getRowSize());
System.arraycopy(retMat_j, 0, ans.denseElements, j, retMat_j.length);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || b.isDiagonal() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.getDiagonalElement(j).multiply(b.diagonalElements[j]);
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.getDiagonalElement(j).multiply(b.diagonalElements[j]);
ans.diagonalElements[j + 1] = a.getDiagonalElement(j + 1).multiply(b.diagonalElements[j + 1]);
ans.diagonalElements[j + 2] = a.getDiagonalElement(j + 2).multiply(b.diagonalElements[j + 2]);
ans.diagonalElements[j + 3] = a.getDiagonalElement(j + 3).multiply(b.diagonalElements[j + 3]);
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* TODO MULTIPLY_NON_ATLASがFalseの場合でやるのもいいかも。 falseの場合でやると、配列のコピーの回数がここの分だけ減る。 Dense=Dense*Sparse
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getColumnSize() != b.getRowSize()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (b.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
for (int index = 0; index < b.nonZeroCount; ++index) {
int i = b.rowIndex[index];
int j = b.columnIndex[index];
RS value = b.getSparseElement(index).multiply(scalar);
if (i != j) {
RS[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_i, 1, retMat_j, 1);
System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
RS[] retMat_i = unit.createArray(ans.denseElements.length - ans.getRowSize() * i);
System.arraycopy(ans.denseElements, ans.getRowSize() * i, retMat_i, 0, retMat_i.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_i, 1);
System.arraycopy(retMat_i, 0, ans.denseElements, ans.getRowSize() * i, retMat_i.length);
} else {
RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_j, 1);
System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || a.isDiagonal() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.getDiagonalElement(j)));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.getDiagonalElement(j)));
ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.getDiagonalElement(j + 1)));
ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.getDiagonalElement(j + 2)));
ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.getDiagonalElement(j + 3)));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* TODO MULTIPLY_NON_ATLASがFalseの場合でやるのもいいかも。 falseの場合でやると、配列のコピーの回数がここの分だけ減る。 Dense=Dense*Sparse
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
if (a.getColumnSize() != b.getRowSize()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (b.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
for (int index = 0; index < b.nonZeroCount; ++index) {
int i = b.rowIndex[index];
int j = b.columnIndex[index];
RS value = b.getSparseElement(index);
if (i != j) {
RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
RS[] retMat_i = unit.createArray(ans.denseElements.length - ans.getRowSize() * i);
System.arraycopy(ans.denseElements, ans.getRowSize() * i, retMat_i, 0, retMat_i.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_i, 1);
System.arraycopy(retMat_i, 0, ans.denseElements, ans.getRowSize() * i, retMat_i.length);
RS[] aMat_i = unit.createArray(a.denseElements.length - a.getRowSize() * i);
System.arraycopy(a.denseElements, a.getRowSize() * i, aMat_i, 0, aMat_i.length);
RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_i, 1, retMat_j, 1);
System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
} else {
RS[] aMat_j = unit.createArray(a.denseElements.length - a.getRowSize() * j);
RS[] retMat_j = unit.createArray(ans.denseElements.length - ans.getRowSize() * j);
System.arraycopy(a.denseElements, a.getRowSize() * j, aMat_j, 0, aMat_j.length);
System.arraycopy(ans.denseElements, ans.getRowSize() * j, retMat_j, 0, retMat_j.length);
BLAS.daxpy(b.getColumnSize(), value, aMat_j, 1, retMat_j, 1);
System.arraycopy(retMat_j, 0, ans.denseElements, ans.getRowSize() * j, retMat_j.length);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
BLAS.dgemm(
"NoTranspose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || a.isDiagonal() == false) {
throw new RuntimeException("multiply :: different matrix type"); //$NON-NLS-1$
}
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.getDiagonalElement(j));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.getDiagonalElement(j));
ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.getDiagonalElement(j + 1));
ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.getDiagonalElement(j + 2));
ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.getDiagonalElement(j + 3));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, RS scalar) {
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
switch (ans.getDenseOrDiagonal()) {
case DENSE:
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dscal(length, scalar, ans.denseElements, 1);
return ans;
case DIAGONAL:
BLAS.dscal(ans.getColumnSize(), scalar, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, RS scalar) {
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> multiply(Vector<RS,RM,CS,CM> a, RS scalar) {
Vector<RS,RM,CS,CM> ans = a.createClone();
BLAS.dscal(ans.getDimension(), scalar, ans.getElements(), 1);
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> multiply(BlockVector<RS,RM,CS,CM> a, RS scalar) {
BlockVector<RS,RM,CS,CM> ans = a.createClone();
for (int i = 0; i < a.getBlockSize(); ++i) {
Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
if (a.getColumnSize() != b.getDimension()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
Vector<RS,RM,CS,CM> ans = b.createClone();
final RS unit = a.getElementUnit();
switch (a.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemv("NoTranspose", a.getRowSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.getElements(), 1, unit.create(0), ans.getElements(), 1); //$NON-NLS-1$
return ans;
case DIAGONAL:
int shou = ans.getDimension() / 4;
int amari = ans.getDimension() % 4;
for (int j = 0; j < amari; ++j) {
ans.setElement(j, a.diagonalElements[j].multiply(b.getElement(j)));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.setElement(j, a.diagonalElements[j].multiply(b.getElement(j)));
ans.setElement(j + 1, a.diagonalElements[j + 1].multiply(b.getElement(j + 1)));
ans.setElement(j + 2, a.diagonalElements[j + 2].multiply(b.getElement(j + 2)));
ans.setElement(j + 3, a.diagonalElements[j + 3].multiply(b.getElement(j + 3)));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b b
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> multiply(DenseMatrix<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b, RS scalar) {
if (a.getColumnSize() != b.getDimension()) {
throw new RuntimeException("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
Vector<RS,RM,CS,CM> ans = new Vector<>(a.getRowSize(), unit);
switch (a.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemv("NoTranspose", a.getRowSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.getElements(), 1, unit.create(0), ans.getElements(), 1); //$NON-NLS-1$
return ans;
case DIAGONAL:
int shou = ans.getDimension() / 4;
int amari = ans.getDimension() % 4;
for (int j = 0; j < amari; ++j) {
ans.setElement(j, scalar.multiply(a.diagonalElements[j].multiply(b.getElement(j))));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.setElement(j, scalar.multiply(a.diagonalElements[j].multiply(b.getElement(j))));
ans.setElement(j + 1, scalar.multiply(a.diagonalElements[j + 1].multiply(b.getElement(j + 1))));
ans.setElement(j + 2, scalar.multiply(a.diagonalElements[j + 2].multiply(b.getElement(j + 2))));
ans.setElement(j + 3, scalar.multiply(a.diagonalElements[j + 3].multiply(b.getElement(j + 3))));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b b
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);
for (int i = 0; i < a.getBlockSize(); ++i) {
Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
if (b.getBlockSize() != a.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockVector<RS,RM,CS,CM> ans = b.createClone();
for (int i = 0; i < b.getBlockSize(); ++i) {
Vector<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = a.createClone();
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = aMat**T * bMat
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> tran_multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemm(
"Transpose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getColumnSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = aMat**T * bMat
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> tran_multiply(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getColumnSize(), b.getColumnSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
BLAS.dgemm(
"Transpose", "NoTranspose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getColumnSize(), b.denseElements, b.getRowSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> tran_multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = tran_multiply(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> tran_multiply(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = tran_multiply(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = aMat * bMat**T
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply_tran(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
// The Point is the first argument is "NoTranspose".
BLAS.dgemm(
"NoTranspose", "Transpose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), scalar, a.denseElements, a.getRowSize(), b.denseElements, b.getColumnSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = scalar.multiply(a.diagonalElements[j].multiply(b.diagonalElements[j]));
ans.diagonalElements[j + 1] = scalar.multiply(a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]));
ans.diagonalElements[j + 2] = scalar.multiply(a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]));
ans.diagonalElements[j + 3] = scalar.multiply(a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]));
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = aMat * bMat**T
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> multiply_tran(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getColumnSize() != b.getColumnSize()) {
Tools.error("multiply :: different matrix size"); //$NON-NLS-1$
}
final RS unit = a.getElementUnit();
DenseMatrix<RS,RM,CS,CM> ans = new DenseMatrix<>(a.getRowSize(), b.getRowSize(), a.getDenseOrDiagonal(), unit);
switch (ans.getDenseOrDiagonal()) {
case DENSE:
// The Point is the first argument is "NoTranspose".
BLAS.dgemm(
"NoTranspose", "Transpose", ans.getRowSize(), ans.getColumnSize(), a.getColumnSize(), unit.create(1), a.denseElements, a.getRowSize(), b.denseElements, b.getColumnSize(), unit.create(0), ans.denseElements, //$NON-NLS-1$ //$NON-NLS-2$
ans.getRowSize());
return ans;
case DIAGONAL:
int shou = ans.getColumnSize() / 4;
int amari = ans.getColumnSize() % 4;
for (int j = 0; j < amari; ++j) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
}
int count;
int j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
ans.diagonalElements[j] = a.diagonalElements[j].multiply(b.diagonalElements[j]);
ans.diagonalElements[j + 1] = a.diagonalElements[j + 1].multiply(b.diagonalElements[j + 1]);
ans.diagonalElements[j + 2] = a.diagonalElements[j + 2].multiply(b.diagonalElements[j + 2]);
ans.diagonalElements[j + 3] = a.diagonalElements[j + 3].multiply(b.diagonalElements[j + 3]);
}
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply_tran(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("multiply:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = multiply_tran(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> multiply_tran(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
RS scalar = null;
BlockDenseMatrix<RS,RM,CS,CM> ans = multiply_tran(a, b, scalar);
return ans;
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> plus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b, RS scalar) {
if (a.getDimension() != b.getDimension()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
Vector<RS,RM,CS,CM> ans = a.createClone();
BLAS.daxpy(ans.getDimension(), scalar, b.getElements(), 1, ans.getElements(), 1);
return ans;
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> plus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
if (a.getDimension() != b.getDimension()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
Vector<RS,RM,CS,CM> ans = a.createClone();
BLAS.dxpy(ans.getDimension(), b.getElements(), 1, ans.getElements(), 1);
return ans;
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> minus(Vector<RS,RM,CS,CM> a, Vector<RS,RM,CS,CM> b) {
if (a.getDimension() != b.getDimension()) {
throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
}
Vector<RS,RM,CS,CM> ans = a.createClone();
BLAS.dxmy(ans.getDimension(), b.getElements(), 1, ans.getElements(), 1);
return ans;
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
switch (a.getDenseOrDiagonal()) {
case DENSE:
int length = a.getRowSize() * a.getColumnSize();
BLAS.daxpy(length, scalar, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
BLAS.daxpy(ans.getColumnSize(), scalar, b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
final DenseMatrix<RS,RM,CS,CM> ans;
switch (a.getDenseOrDiagonal()) {
case DENSE:
int length = a.getRowSize() * a.getColumnSize();
ans = a.createClone();
BLAS.dxpy(length, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
ans = a.createClone();
BLAS.dxpy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> minus(DenseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize() || a.getDenseOrDiagonal() != b.getDenseOrDiagonal()) {
throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
}
final DenseMatrix<RS,RM,CS,CM> ans;
switch (a.getDenseOrDiagonal()) {
case DENSE:
int length = a.getRowSize() * a.getColumnSize();
ans = a.createClone();
BLAS.dxmy(length, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
ans = a.createClone();
BLAS.dxmy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
// ret = (*scalar) * b
DenseMatrix<RS,RM,CS,CM> ans = multiply(b, scalar);
// ret += a
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int shou = a.nonZeroCount / 4;
int amari = a.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].add(value);
ans.denseElements[ji] = ans.denseElements[ji].add(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].add(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = a.rowIndex[index];
int j1 = a.columnIndex[index];
RS value1 = a.getSparseElement(index);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].add(value1);
ans.denseElements[ji] = ans.denseElements[ji].add(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].add(value1);
}
int i2 = a.rowIndex[index + 1];
int j2 = a.columnIndex[index + 1];
RS value2 = a.getSparseElement(index + 1);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].add(value2);
ans.denseElements[ji] = ans.denseElements[ji].add(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].add(value2);
}
int i3 = a.rowIndex[index + 2];
int j3 = a.columnIndex[index + 2];
RS value3 = a.getSparseElement(index + 2);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].add(value3);
ans.denseElements[ji] = ans.denseElements[ji].add(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].add(value3);
}
int i4 = a.rowIndex[index + 3];
int j4 = a.columnIndex[index + 3];
RS value4 = a.getSparseElement(index + 3);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].add(value4);
ans.denseElements[ji] = ans.denseElements[ji].add(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].add(value4);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || b.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dxpy(length, a.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || b.isDiagonal() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
BLAS.dxpy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
DenseMatrix<RS,RM,CS,CM> ans = b.createClone();
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (b.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int shou = a.nonZeroCount / 4;
int amari = a.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].add(value);
ans.denseElements[ji] = ans.denseElements[ji].add(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].add(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = a.rowIndex[index];
int j1 = a.columnIndex[index];
RS value1 = a.getSparseElement(index);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].add(value1);
ans.denseElements[ji] = ans.denseElements[ji].add(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].add(value1);
}
int i2 = a.rowIndex[index + 1];
int j2 = a.columnIndex[index + 1];
RS value2 = a.getSparseElement(index + 1);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].add(value2);
ans.denseElements[ji] = ans.denseElements[ji].add(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].add(value2);
}
int i3 = a.rowIndex[index + 2];
int j3 = a.columnIndex[index + 2];
RS value3 = a.getSparseElement(index + 2);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].add(value3);
ans.denseElements[ji] = ans.denseElements[ji].add(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].add(value3);
}
int i4 = a.rowIndex[index + 3];
int j4 = a.columnIndex[index + 3];
RS value4 = a.getSparseElement(index + 3);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].add(value4);
ans.denseElements[ji] = ans.denseElements[ji].add(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].add(value4);
}
}
return ans;
case DENSE:
if (b.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dxpy(length, a.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
if (b.isDiagonal() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
BLAS.dxpy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> minus(SparseMatrix<RS,RM,CS,CM> a, DenseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
}
DenseMatrix<RS,RM,CS,CM> ans = b.createClone();
switch (a.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (b.isDense() == false) {
throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
}
int shou = a.nonZeroCount / 4;
int amari = a.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = a.rowIndex[index];
int j = a.columnIndex[index];
RS value = a.getSparseElement(index);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = a.rowIndex[index];
int j1 = a.columnIndex[index];
RS value1 = a.getSparseElement(index);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value1);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value1);
}
int i2 = a.rowIndex[index + 1];
int j2 = a.columnIndex[index + 1];
RS value2 = a.getSparseElement(index + 1);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value2);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value2);
}
int i3 = a.rowIndex[index + 2];
int j3 = a.columnIndex[index + 2];
RS value3 = a.getSparseElement(index + 2);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value3);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value3);
}
int i4 = a.rowIndex[index + 3];
int j4 = a.columnIndex[index + 3];
RS value4 = a.getSparseElement(index + 3);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value4);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value4);
}
}
return ans;
case DENSE:
if (b.isDense() == false) {
throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
}
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dxmy(length, a.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
if (b.isDiagonal() == false) {
throw new RuntimeException("minus :: different matrix type"); //$NON-NLS-1$
}
BLAS.dxmy(ans.getColumnSize(), a.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
// ret = a
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
//ret += (*scalar) * b
switch (b.getSparseOrDenseOrDiagonal()) {
case SPARSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int shou = b.nonZeroCount / 4;
int amari = b.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = b.rowIndex[index];
int j = b.columnIndex[index];
RS value = b.getSparseElement(index).multiply(scalar);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].add(value);
ans.denseElements[ji] = ans.denseElements[ji].add(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].add(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = b.rowIndex[index];
int j1 = b.columnIndex[index];
RS value1 = b.getSparseElement(index).multiply(scalar);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].add(value1);
ans.denseElements[ji] = ans.denseElements[ji].add(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].add(value1);
}
int i2 = b.rowIndex[index + 1];
int j2 = b.columnIndex[index + 1];
RS value2 = b.getSparseElement(index + 1).multiply(scalar);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].add(value2);
ans.denseElements[ji] = ans.denseElements[ji].add(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].add(value2);
}
int i3 = b.rowIndex[index + 2];
int j3 = b.columnIndex[index + 2];
RS value3 = b.getSparseElement(index + 2).multiply(scalar);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].add(value3);
ans.denseElements[ji] = ans.denseElements[ji].add(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].add(value3);
}
int i4 = b.rowIndex[index + 3];
int j4 = b.columnIndex[index + 3];
RS value4 = b.getSparseElement(index + 3).multiply(scalar);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].add(value4);
ans.denseElements[ji] = ans.denseElements[ji].add(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].add(value4);
}
}
return ans;
case DENSE:
if (ans.isDense() == false || a.isDense() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.daxpy(length, scalar, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
if (ans.isDiagonal() == false || a.isDiagonal() == false) {
throw new RuntimeException("plus :: different matrix type"); //$NON-NLS-1$
}
BLAS.daxpy(ans.getColumnSize(), scalar, b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> plus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("plus :: different matrix size"); //$NON-NLS-1$
}
// ans = a
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
// ans += (*scalar) * b
switch (b.getSparseOrDenseOrDiagonal()) {
case SPARSE:
int shou = b.nonZeroCount / 4;
int amari = b.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = b.rowIndex[index];
int j = b.columnIndex[index];
RS value = b.getSparseElement(index);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].add(value);
ans.denseElements[ji] = ans.denseElements[ji].add(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].add(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = b.rowIndex[index];
int j1 = b.columnIndex[index];
RS value1 = b.getSparseElement(index);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].add(value1);
ans.denseElements[ji] = ans.denseElements[ji].add(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].add(value1);
}
int i2 = b.rowIndex[index + 1];
int j2 = b.columnIndex[index + 1];
RS value2 = b.getSparseElement(index + 1);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].add(value2);
ans.denseElements[ji] = ans.denseElements[ji].add(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].add(value2);
}
int i3 = b.rowIndex[index + 2];
int j3 = b.columnIndex[index + 2];
RS value3 = b.getSparseElement(index + 2);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].add(value3);
ans.denseElements[ji] = ans.denseElements[ji].add(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].add(value3);
}
int i4 = b.rowIndex[index + 3];
int j4 = b.columnIndex[index + 3];
RS value4 = b.getSparseElement(index + 3);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].add(value4);
ans.denseElements[ji] = ans.denseElements[ji].add(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].add(value4);
}
}
return ans;
case DENSE:
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dxpy(length, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
BLAS.dxpy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> minus(DenseMatrix<RS,RM,CS,CM> a, SparseMatrix<RS,RM,CS,CM> b) {
if (a.getRowSize() != b.getRowSize() || a.getColumnSize() != b.getColumnSize()) {
throw new RuntimeException("minus :: different matrix size"); //$NON-NLS-1$
}
DenseMatrix<RS,RM,CS,CM> ans = a.createClone();
switch (b.getSparseOrDenseOrDiagonal()) {
case SPARSE:
int shou = b.nonZeroCount / 4;
int amari = b.nonZeroCount % 4;
for (int index = 0; index < amari; ++index) {
int i = b.rowIndex[index];
int j = b.columnIndex[index];
RS value = b.getSparseElement(index);
if (i != j) {
int ij = i + ans.getColumnSize() * j;
int ji = j + ans.getColumnSize() * i;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value);
} else {
int ii = i + ans.getColumnSize() * i;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value);
}
}
int index;
int counter;
for (index = amari, counter = 0; counter < shou; ++counter, index += 4) {
int i1 = b.rowIndex[index];
int j1 = b.columnIndex[index];
RS value1 = b.getSparseElement(index);
if (i1 != j1) {
int ij = i1 + ans.getColumnSize() * j1;
int ji = j1 + ans.getColumnSize() * i1;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value1);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value1);
} else {
int ii = i1 + ans.getColumnSize() * i1;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value1);
}
int i2 = b.rowIndex[index + 1];
int j2 = b.columnIndex[index + 1];
RS value2 = b.getSparseElement(index + 1);
if (i2 != j2) {
int ij = i2 + ans.getColumnSize() * j2;
int ji = j2 + ans.getColumnSize() * i2;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value2);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value2);
} else {
int ii = i2 + ans.getColumnSize() * i2;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value2);
}
int i3 = b.rowIndex[index + 2];
int j3 = b.columnIndex[index + 2];
RS value3 = b.getSparseElement(index + 2);
if (i3 != j3) {
int ij = i3 + ans.getColumnSize() * j3;
int ji = j3 + ans.getColumnSize() * i3;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value3);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value3);
} else {
int ii = i3 + ans.getColumnSize() * i3;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value3);
}
int i4 = b.rowIndex[index + 3];
int j4 = b.columnIndex[index + 3];
RS value4 = b.getSparseElement(index + 3);
if (i4 != j4) {
int ij = i4 + ans.getColumnSize() * j4;
int ji = j4 + ans.getColumnSize() * i4;
ans.denseElements[ij] = ans.denseElements[ij].subtract(value4);
ans.denseElements[ji] = ans.denseElements[ji].subtract(value4);
} else {
int ii = i4 + ans.getColumnSize() * i4;
ans.denseElements[ii] = ans.denseElements[ii].subtract(value4);
}
}
return ans;
case DENSE:
int length = ans.getRowSize() * ans.getColumnSize();
BLAS.dxmy(length, b.denseElements, 1, ans.denseElements, 1);
return ans;
case DIAGONAL:
BLAS.dxmy(ans.getColumnSize(), b.diagonalElements, 1, ans.diagonalElements, 1);
return ans;
default:
throw new IllegalArgumentException();
}
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> plus(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
RS unit = a.getElementUnit();
BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);
for (int i = 0; i < a.getBlockSize(); ++i) {
Vector<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param b b
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> plus(BlockVector<RS,RM,CS,CM> a, BlockVector<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
RS unit = a.getElementUnit();
BlockVector<RS,RM,CS,CM> ans = new BlockVector<>(a.getBlockSize(), a.getBlockStructs(), unit);
for (int i = 0; i < a.getBlockSize(); ++i) {
Vector<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockSparseMatrix<RS,RM,CS,CM> a, BlockDenseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(b.getBlockSize(), b.getBlockStructs(), b.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans + a * b*scalar
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @param scalar スカラー
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i), scalar);
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a + b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> plus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("plus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = plus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a - b
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param b B
* @return 解
*/
private static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> minus(BlockDenseMatrix<RS,RM,CS,CM> a, BlockSparseMatrix<RS,RM,CS,CM> b) {
if (a.getBlockSize() != b.getBlockSize()) {
throw new RuntimeException("minus:: different nBlock size"); //$NON-NLS-1$
}
BlockDenseMatrix<RS,RM,CS,CM> ans = new BlockDenseMatrix<>(a.getBlockSize(), a.getBlockStructs(), a.getUnit());
for (int i = 0; i < a.getBlockSize(); ++i) {
DenseMatrix<RS,RM,CS,CM> blockAns = minus(a.getBlock(i), b.getBlock(i));
ans.setBlock(i, blockAns);
}
return ans;
}
/**
* ans = a '*' (*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, RS scalar) {
switch (operator) {
case '*':
Vector<RS,RM,CS,CM> ans = multiply(a, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '*' (*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockVector<RS,RM,CS,CM> let(BlockVector<RS,RM,CS,CM> a, final char operator, RS scalar) {
switch (operator) {
case '*':
BlockVector<RS,RM,CS,CM> ans = multiply(a, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '*' (*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, RS scalar) {
switch (operator) {
case '*':
BlockDenseMatrix<RS,RM,CS,CM> ans = multiply(a, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param b b
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b, RS scalar) {
Vector<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param b b
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> let(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
Vector<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' 't' 'T' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
case 't':
// ret = aMat**T * bMat
ans = tran_multiply(a, b, scalar);
return ans;
case 'T':
// ret = aMat * bMat**T
ans = multiply_tran(a, b, scalar);
return ans;
default:
throw new IllegalArgumentException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' 't' 'T' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
case 't':
// ret = aMat**T * bMat
ans = tran_multiply(a, b);
return ans;
case 'T':
// ret = aMat * bMat**T
ans = multiply_tran(a, b);
return ans;
default:
throw new IllegalArgumentException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b, RS scalar) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b, RS scalar) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> DenseMatrix<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b) {
DenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' 't' 'T' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
case 't':
// ret = aMat**T * bMat
ans = tran_multiply(a, b, scalar);
return ans;
case 'T':
// ret = aMat * bMat**T
ans = multiply_tran(a, b, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' 't' 'T' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
case 't':
// ret = aMat**T * bMat
ans = tran_multiply(a, b);
return ans;
case 'T':
// ret = aMat * bMat**T
ans = multiply_tran(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b, RS scalar) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @param scalar スカラー
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b, RS scalar) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b, scalar);
return ans;
case '-':
final RS minus_scalar = scalar.unaryMinus();
ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b, scalar);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = a '+' '-' '*' b*(*scalar)
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BlockDenseMatrix<RS,RM,CS,CM> let(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b) {
BlockDenseMatrix<RS,RM,CS,CM> ans;
switch (operator) {
case '+':
ans = plus(a, b);
return ans;
case '-':
ans = minus(a, b);
//final RS unit = a.getElementUnit();
//final RS minus_scalar = unit.create(-1);
//ans = plus(a, b, minus_scalar);
return ans;
case '*':
ans = multiply(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = aMat '*' '/' bVec
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b b
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> Vector<RS,RM,CS,CM> let(DenseMatrix<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
Vector<RS,RM,CS,CM> ans;
switch (operator) {
case '*':
ans = multiply(a, b);
return ans;
case '/':
// ret = aMat^{-1} * bVec;
// aMat is positive definite and already cholesky factorized.
ans = solveSystems(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param b b
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(Vector<RS,RM,CS,CM> a, final char operator, Vector<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(DenseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(DenseMatrix<RS,RM,CS,CM> a, final char operator, SparseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(b, a);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(SparseMatrix<RS,RM,CS,CM> a, final char operator, DenseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a a
* @param operator 演算子
* @param b b
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(BlockVector<RS,RM,CS,CM> a, final char operator, BlockVector<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(BlockSparseMatrix<RS,RM,CS,CM> a, final char operator, BlockDenseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(a, b);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
/**
* ans = inner_product(a,b) // op = '.'
*
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
* @param a A
* @param operator 演算子
* @param b B
* @return 解
*/
public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RS run(BlockDenseMatrix<RS,RM,CS,CM> a, final char operator, BlockSparseMatrix<RS,RM,CS,CM> b) {
switch (operator) {
case '.':
RS ans = getInnerProduct(b, a);
return ans;
default:
throw new RuntimeException("let:: operator error"); //$NON-NLS-1$
}
}
}