Newton.java
package org.mklab.sdpj.algorithm;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.PrintStream;
import org.mklab.cga.round.FPURoundModeSelector;
import org.mklab.gap.Guarantor;
import org.mklab.gap.MPFloatIntervalMatrix;
import org.mklab.gap.linear.MPFloatIntlabMethod;
import org.mklab.mpf.ap.MPFloatMatrix;
import org.mklab.mpf.ap.MPFloatRoundModeSelector;
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.nfc.util.RoundModeManager;
import org.mklab.sdpj.algebra.Algebra;
import org.mklab.sdpj.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.BooleanObject;
import org.mklab.sdpj.datastructure.DenseDiagonal;
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.tool.Tools;
/**
* @author koga
* @version $Revision$, 2009/04/24
* @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 Newton<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {
/** */
private DenseMatrix<RS,RM,CS,CM> bMat;
/** */
private DenseMatrix<RS,RM,CS,CM> inverseOfbMat;
/** */
private Vector<RS,RM,CS,CM> gVec;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> DxMat;
/** */
private Vector<RS,RM,CS,CM> DyVec;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> DzMat;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> rMat;
/** */
private int[] upNonZeroCount;
/** */
private FormulaType[] useFormula;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> invzMat;
// temporary matrix
/** */
private BlockDenseMatrix<RS,RM,CS,CM> fMat;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> gMat;
/** */
private BlockDenseMatrix<RS,RM,CS,CM> DxMatDzMat;
/**
* コンストラクタ
*
* @param m ,
* @param blockSize ,
* @param blockStruct .
* @param unit 単位
*/
public Newton(int m, int blockSize, int[] blockStruct, RS unit) {
this.bMat = new DenseMatrix<>(m, m, DenseDiagonal.DENSE, unit);
this.gVec = new Vector<>(m, unit.create(0));
this.DxMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.DyVec = new Vector<>(m, unit.create(0));
this.DzMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.rMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.invzMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.fMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.gMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.DxMatDzMat = new BlockDenseMatrix<>(blockSize, blockStruct, unit);
this.upNonZeroCount = new int[m * blockSize];
this.useFormula = new FormulaType[m * blockSize];
}
/**
* get instance bMat
*
* @return bMat
*/
public DenseMatrix<RS,RM,CS,CM> getBMat() {
return this.bMat;
}
/**
* dxMatを返します。
*
* @return dxMat
*/
public BlockDenseMatrix<RS,RM,CS,CM> getDxMat() {
return this.DxMat;
}
/**
* dyVecを返します。
*
* @return dyVec
*/
public Vector<RS,RM,CS,CM> getDyVec() {
return this.DyVec;
}
/**
* dzMatを返します。
*
* @return dzMat
*/
public BlockDenseMatrix<RS,RM,CS,CM> getDzMat() {
return this.DzMat;
}
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 633
*
* @param m ,
* @param A ,
* @param kappa .
*/
public void computeFormula(int m, BlockSparseMatrix<RS,RM,CS,CM>[] A, RS kappa) {
final RS unit = kappa.createUnit();
if (this.upNonZeroCount == null || this.useFormula == null) {
Tools.error("rNewton:: failed initialization"); //$NON-NLS-1$
}
// Count sum of number of elements
// that each number of elements are less than own.
int blockSize = A[0].getBlockSize();
for (int i = 0; i < blockSize; ++i) {
if (A[0].getTargetBlockStruct(i) < 0) {
// in Diagonal case, we don't have to calculate.
continue;
}
for (int j = 0; j < m; ++j) {
int up = A[j].getBlock(i).nonZeroEffect;
for (int k = 0; k < m; ++k) {
if (A[k].getBlock(i).nonZeroEffect < A[j].getBlock(i).nonZeroEffect) {
up += A[k].getBlock(i).nonZeroEffect;
} else if (A[k].getBlock(i).nonZeroEffect == A[j].getBlock(i).nonZeroEffect && k < j) {
up += A[k].getBlock(i).nonZeroEffect;
}
}
this.upNonZeroCount[j * blockSize + i] = up;
}
}
// Determine which formula
for (int i = 0; i < blockSize; ++i) {
if (A[0].getTargetBlockStruct(i) < 0) {
// in Diagonal case, we don't have to calculate.
continue;
}
for (int j = 0; j < m; ++j) {
RS n = unit.create(A[j].getBlock(i).getRowSize());
RS up = unit.create(this.upNonZeroCount[j * blockSize + i]);
RS nonzero = unit.create(A[j].getBlock(i).nonZeroEffect);
RS f1 = (kappa.multiply(n).multiply(nonzero)).add(n.power(3)).add(kappa.multiply(up));
RS f2 = (kappa.multiply(n.multiply(nonzero))).add(kappa.multiply((n.add(unit.create(1))).multiply(up)));
RS f3 = kappa.multiply((kappa.multiply(unit.create(2).multiply(nonzero)).add(1))).multiply(up);
if (A[j].getBlock(i).isDense()) {
// if DENSE, we use only F1 or F2, that is we don't use F3
if (f1.isLessThan(f2)) {
this.useFormula[j * blockSize + i] = FormulaType.F1;
} else {
this.useFormula[j * blockSize + i] = FormulaType.F2;
}
} else {
// this case is SPARSE
if (f1.isLessThan(f2) && f1.isLessThan(f3)) {
this.useFormula[j * blockSize + i] = FormulaType.F1;
} else if (f2.isLessThan(f3)) {
this.useFormula[j * blockSize + i] = FormulaType.F2;
} else {
this.useFormula[j * blockSize + i] = FormulaType.F3;
}
}
}
}
return;
}
/**
* BLAS.ddot (int,NumericalScalar,int,NumericalScalar[],int)とは別に(int,N,int,N,int )に対応できるようにする。
*
* @param G G
* @param Aj Aj
* @return F1
*/
private RS calF1(DenseMatrix<RS,RM,CS,CM> G, SparseMatrix<RS,RM,CS,CM> Aj) {
return Algebra.run(Aj, '.', G);
}
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 752
*
* @param F F
* @param G G
* @param X X
* @param Aj Aj
* @param hasF2Gcal hasF2Gcal
* @return F2
*/
private RS calF2(DenseMatrix<RS,RM,CS,CM> F, DenseMatrix<RS,RM,CS,CM> G, DenseMatrix<RS,RM,CS,CM> X, SparseMatrix<RS,RM,CS,CM> Aj, BooleanObject hasF2Gcal) {
final RS unit = Aj.getElementUnit();
RS ret = unit.create(0);
int n = Aj.getRowSize();
int index = 0;
switch (Aj.getSparseOrDenseOrDiagonal()) {
case SPARSE:
for (index = 0; index < Aj.nonZeroCount; ++index) {
int alpha = Aj.rowIndex[index];
int beta = Aj.columnIndex[index];
RS value1 = Aj.getSparseElement(index);
RS[] x_alpha = unit.createArray(X.denseElements.length - alpha);
RS[] f_nbeta = unit.createArray(F.denseElements.length - n * beta);
System.arraycopy(X.denseElements, alpha, x_alpha, 0, x_alpha.length);
System.arraycopy(F.denseElements, n * beta, f_nbeta, 0, f_nbeta.length);
RS value2 = BLAS.ddot(n, x_alpha, n, f_nbeta, 1);
ret = ret.add(value1.multiply(value2));
if (alpha != beta) {
RS[] x_beta = unit.createArray(X.denseElements.length - beta);
RS[] f_nalpha = unit.createArray(F.denseElements.length - n * alpha);
System.arraycopy(X.denseElements, beta, x_beta, 0, x_beta.length);
System.arraycopy(F.denseElements, n * alpha, f_nalpha, 0, f_nalpha.length);
value2 = BLAS.<RS, RM, CS, CM> ddot(n, x_beta, n, f_nalpha, 1);
ret = ret.add(value1.multiply(value2));
}
}
break;
case DENSE:
// G is temporary matrix
if (hasF2Gcal.getBoolean() == false) {
DenseMatrix<RS,RM,CS,CM> ans = Algebra.let(X, '*', F);
G.copyFrom(ans);
hasF2Gcal.setBoolean(true);
}
ret = Algebra.run(Aj, '.', G);
break;
case DIAGONAL:
// DIAGONAL is something wrong
Tools.message("F2::DIAGONAL"); //$NON-NLS-1$
Tools.message("calF2:: miss condition"); //$NON-NLS-1$
break;
default:
throw new UnsupportedOperationException();
}
return ret;
}
//
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 803
*
* @param X X
* @param invZ invZ
* @param Ai Ai
* @param Aj Aj
* @return F3
*/
private RS calF3(DenseMatrix<RS,RM,CS,CM> X, DenseMatrix<RS,RM,CS,CM> invZ, SparseMatrix<RS,RM,CS,CM> Ai, SparseMatrix<RS,RM,CS,CM> Aj) {
final RS unit = Aj.getElementUnit();
// Ai and Aj are SPARSE
RS ret = unit.createZero();
for (int i = 0; i < Aj.nonZeroCount; ++i) {
int alpha = Aj.rowIndex[i];
int beta = Aj.columnIndex[i];
RS value1 = Aj.getSparseElement(i);
RS sum = unit.createZero();
for (int j = 0; j < Ai.nonZeroCount; ++j) {
int gamma = Ai.rowIndex[j];
int delta = Ai.columnIndex[j];
RS value2 = Ai.getSparseElement(j);
RS plu = value2.multiply(invZ.denseElements[delta + invZ.getColumnSize() * beta]).multiply(X.denseElements[alpha + X.getColumnSize() * gamma]);
sum = sum.add(plu);
if (gamma != delta) {
RS plu2 = value2.multiply(invZ.denseElements[gamma + invZ.getColumnSize() * beta].multiply(X.denseElements[alpha + X.getColumnSize() * delta]));
sum = sum.add(plu2);
}
}
ret = ret.add(value1.multiply(sum));
if (alpha == beta) {
continue;
}
sum = unit.createZero();
for (int j = 0; j < Ai.nonZeroCount; ++j) {
int gamma = Ai.rowIndex[j];
int delta = Ai.columnIndex[j];
RS value2 = Ai.getSparseElement(j);
RS plu = value2.multiply(invZ.denseElements[delta + invZ.getColumnSize() * alpha]).multiply(X.denseElements[beta + X.getColumnSize() * gamma]);
sum = sum.add(plu);
if (gamma != delta) {
RS plu2 = value2.multiply(invZ.denseElements[gamma + invZ.getColumnSize() * alpha]).multiply(X.denseElements[beta + X.getColumnSize() * delta]);
sum = sum.add(plu2);
}
}
ret = ret.add(value1.multiply(sum));
}
return ret;
}
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 1007
*
* @param m m
* @param A A
* @param xMat xMat
* @param invzMat invzMat
*/
private void compute_bMat(int m, BlockSparseMatrix<RS,RM,CS,CM>[] A, BlockDenseMatrix<RS,RM,CS,CM> xMat, BlockDenseMatrix<RS,RM,CS,CM> invzMat) {
int blockSize = A[0].getBlockSize();
this.bMat.setZero();
for (int k = 0; k < blockSize; ++k) {
if (A[0].getTargetBlockStruct(k) < 0) {
// case Diagonal
for (int i = 0; i < m; ++i) {
this.fMat.setBlock(k, Algebra.let(A[i].getBlock(k), '*', invzMat.getBlock(k)));
this.gMat.setBlock(k, Algebra.let(xMat.getBlock(k), '*', this.fMat.getBlock(k)));
for (int j = 0; j <= i; ++j) {
RS value = Algebra.run(this.gMat.getBlock(k), '.', A[j].getBlock(k));
if (i != j) {
this.bMat.denseElements[i + this.bMat.getColumnSize() * j] = this.bMat.denseElements[i + this.bMat.getColumnSize() * j].add(value);
this.bMat.denseElements[j + this.bMat.getColumnSize() * i] = this.bMat.denseElements[j + this.bMat.getColumnSize() * i].add(value);
} else {
this.bMat.denseElements[i + this.bMat.getColumnSize() * i] = this.bMat.denseElements[i + this.bMat.getColumnSize() * i].add(value);
}
}
}
} else {
// case not Diagonal
for (int i = 0; i < m; ++i) {
final int A_i_ele_l_NonZeroEffect = A[i].getBlock(k).nonZeroEffect;
if (A_i_ele_l_NonZeroEffect == 0) {
continue;
}
FormulaType formula = this.useFormula[i * blockSize + k];
// ---------------------------------------------------
// formula = F3; // this is force change
// ---------------------------------------------------
BooleanObject hasF2Gcal = new BooleanObject(false);
if (formula == FormulaType.F1) {
this.fMat.setBlock(k, Algebra.let(A[i].getBlock(k), '*', invzMat.getBlock(k)));
this.gMat.setBlock(k, Algebra.let(xMat.getBlock(k), '*', this.fMat.getBlock(k)));
} else if (formula == FormulaType.F2) {
this.fMat.setBlock(k, Algebra.let(A[i].getBlock(k), '*', invzMat.getBlock(k)));
hasF2Gcal.setBoolean(false);
}
for (int j = 0; j < m; ++j) {
final int A_j_ele_l_NonZeroEffect = A[j].getBlock(k).nonZeroEffect;
if (A_j_ele_l_NonZeroEffect == 0) {
continue;
}
// Select the formula A[i] or the formula A[j].
// Use formula that has more NonZeroEffects than others.
// We must calculate i==j.
if (A_i_ele_l_NonZeroEffect < A_j_ele_l_NonZeroEffect || ((A_i_ele_l_NonZeroEffect == A_j_ele_l_NonZeroEffect) && i < j)) {
continue;
}
RS value;
switch (formula) {
case F1:
value = calF1(this.gMat.getBlock(k), A[j].getBlock(k));
break;
case F2:
value = calF2(this.fMat.getBlock(k), this.gMat.getBlock(k), xMat.getBlock(k), A[j].getBlock(k), hasF2Gcal);
break;
case F3:
value = calF3(xMat.getBlock(k), invzMat.getBlock(k), A[i].getBlock(k), A[j].getBlock(k));
break;
default:
throw new IllegalArgumentException();
}
if (i != j) {
this.bMat.denseElements[i + this.bMat.getColumnSize() * j] = this.bMat.denseElements[i + this.bMat.getColumnSize() * j].add(value);
this.bMat.denseElements[j + this.bMat.getColumnSize() * i] = this.bMat.denseElements[j + this.bMat.getColumnSize() * i].add(value);
} else {
this.bMat.denseElements[i + this.bMat.getColumnSize() * i] = this.bMat.denseElements[i + this.bMat.getColumnSize() * i].add(value);
}
}
}
}
}
}
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 1123
*
* @param direction direction
* @param mu mu
* @param beta beta
* @param currentPt currentPt
*/
private void compute_rMat(WHICH_DIRECTION direction, AverageComplementarity<RS, RM, CS, CM> mu, DirectionParameter<RS, RM, CS, CM> beta, Solution<RS,RM,CS,CM> currentPt) {
// CORRECTOR :: R = -XZ -dXdZ + mu I
// not CORRECTOR :: R = -XZ + mu I
final RS unit = beta.getValue().createUnit();
this.rMat = Algebra.let(currentPt.xMatzMat, '*', unit.create(-1));
if (direction == WHICH_DIRECTION.CORRECTOR) {
this.DxMatDzMat = Algebra.let(this.DxMat, '*', this.DzMat);
this.rMat = Algebra.let(this.rMat, '+', this.DxMatDzMat, unit.create(-1));
}
RS target_mu = beta.getValue().multiply(mu.getValue());
for (int i = 0; i < this.rMat.getBlockSize(); ++i) {
int size = this.rMat.getBlockStruct(i);
if (size < 0) {
// case Diagonal
size = -size;
RS[] target = this.rMat.getBlock(i).diagonalElements;
int shou = size / 4;
int amari = size % 4;
for (int j = 0; j < amari; ++j) {
target[j] = target[j].add(target_mu);
}
int count, j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
target[j] = target[j].add(target_mu);
target[j + 1] = target[j + 1].add(target_mu);
target[j + 2] = target[j + 2].add(target_mu);
target[j + 3] = target[j + 3].add(target_mu);
}
} else { // case non Diagonal
RS[] target = this.rMat.getBlock(i).denseElements;
int shou = size / 4;
int amari = size % 4;
for (int j = 0; j < amari; ++j) {
target[j + j * size] = target[j + j * size].add(target_mu);
}
int count, j;
for (j = amari, count = 0; count < shou; ++count, j += 4) {
target[j + j * size] = target[j + j * size].add(target_mu);
target[(j + 1) + (j + 1) * size] = target[(j + 1) + (j + 1) * size].add(target_mu);
target[(j + 2) + (j + 2) * size] = target[(j + 2) + (j + 2) * size].add(target_mu);
target[(j + 3) + (j + 3) * size] = target[(j + 3) + (j + 3) * size].add(target_mu);
}
}
}
}
/**
* Refer to rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 1188
*
* @param direction ,
* @param m ,
* @param A ,
* @param mu ,
* @param beta ,
* @param phase ,
* @param currentPt ,
* @param currentRes ,
* @return boolean
*/
public boolean Mehrotra(WHICH_DIRECTION direction, int m, BlockSparseMatrix<RS,RM,CS,CM>[] A, AverageComplementarity<RS, RM, CS, CM> mu, DirectionParameter<RS, RM, CS, CM> beta, Phase<RS,RM,CS,CM> phase,
Solution<RS,RM,CS,CM> currentPt, Residuals<RS,RM,CS,CM> currentRes) {
final boolean guarantee = false;
final RS unit = A[0].getElementUnit();
if (direction == WHICH_DIRECTION.PREDICTOR) {
// invzMat = Z^{-1}
if (!guarantee) {
for (int i = 0; i < currentPt.invCholeskyZ.getBlockSize(); ++i) {
if (currentPt.invCholeskyZ.getBlock(i).isDense()) {
this.invzMat.setBlock(i, currentPt.invCholeskyZ.getBlock(i).createClone());
BLAS.dtrmm("Left", "Lower", "Transpose", "NonUnitDiag", currentPt.invCholeskyZ.getBlock(i).getRowSize(), currentPt.invCholeskyZ.getBlock(i).getColumnSize(), unit.create(1), //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
currentPt.invCholeskyZ.getBlock(i).denseElements, currentPt.invCholeskyZ.getBlock(i).getRowSize(), this.invzMat.getBlock(i).denseElements, this.invzMat.getBlock(i).getRowSize());
} else {
for (int j = 0; j < this.invzMat.getBlock(i).getRowSize(); ++j) {
this.invzMat.getBlock(i).diagonalElements[j] = currentPt.invCholeskyZ.getBlock(i).diagonalElements[j].multiply(currentPt.invCholeskyZ.getBlock(i).diagonalElements[j]);
}
}
}
} else {
getInverseMatrixOfXWithVerifier(currentPt);
}
}
compute_rMat(direction, mu, beta, currentPt);
if (phase.getValue() == PhaseType.pFEAS || phase.getValue() == PhaseType.noINFO) {
// currentPt is infeasible, that is the residual
// dualMat is not 0. fMat, gMat is temporary, gMat = R-XD.
this.fMat = Algebra.let(currentPt.getXmat(), '*', currentRes.getDualMat());
this.gMat = Algebra.let(this.rMat, '+', this.fMat, unit.create(-1));
} else {
// dualMat == 0
this.gMat = this.rMat.createClone();
}
this.fMat = Algebra.let(this.gMat, '*', this.invzMat, unit.create(-1));
// fMat = (-1)*(R-XD)*Z^{-1}
for (int i = 0; i < m; ++i) {
this.gVec.setElement(i, Algebra.run(this.fMat, '.', A[i]));
}
this.gVec = Algebra.let(this.gVec, '+', currentRes.getPrimalVec());
if (direction == WHICH_DIRECTION.PREDICTOR) {
compute_bMat(m, A, currentPt.getXmat(), this.invzMat);
if (!guarantee) {
boolean ret = Algebra.choleskyFactorWithAdjust(this.bMat);
if (ret == false) {
return false;
}
} else {
getInverseMatrixOfBWithVerifier();
}
}
// bMat is already cholesky factorized.
if (!guarantee) {
this.DyVec = Algebra.let(this.bMat, '/', this.gVec);
} else {
this.DyVec = Algebra.let(this.inverseOfbMat, '*', this.gVec);
}
if (phase.getValue() == PhaseType.pFEAS || phase.getValue() == PhaseType.noINFO) {
this.DzMat = currentRes.getDualMat().createClone();
} else {
this.DzMat.setZero();
}
for (int i = 0; i < m; ++i) {
this.DzMat = Algebra.let(this.DzMat, '-', A[i], this.DyVec.getElement(i));
}
// fMat,gMat are temporary.
this.fMat = Algebra.let(currentPt.getXmat(), '*', this.DzMat);
this.gMat = Algebra.let(this.rMat, '+', this.fMat, unit.create(-1));
this.DxMat = Algebra.let(this.gMat, '*', this.invzMat);
Algebra.getSymmetrize(this.DxMat);
return true;
}
/**
* Compute a inverse matrix of B based on reliable computing
*/
private void getInverseMatrixOfBWithVerifier() {
RoundModeManager manager = RoundModeManager.getManager();
manager.add(new FPURoundModeSelector());
manager.add(new MPFloatRoundModeSelector());
RS[][] matrixA = this.bMat.getArrayOfElements();
RM a = matrixA[0][0].createGrid(matrixA);
//NumericalMatrix<RS> a = new BaseNumericalMatrix<>(matrixA);
RS unit = this.bMat.getElementUnit();
int rowSize = this.bMat.getRowSize();
int columnSize = this.bMat.getColumnSize();
RM b = unit.createUnitGrid(rowSize, columnSize);
MPFloatIntlabMethod verifier = new MPFloatIntlabMethod((MPFloatMatrix)a,(MPFloatMatrix)b);
Guarantor<MPFloatIntlabMethod> guarantor = new Guarantor<>(verifier);
guarantor.solveWithEffectiveDigits(10);
MPFloatIntervalMatrix x = guarantor.getGuranteedVerifier().getSolution();
//IntervalMatrix<RS> x = ((LinearEquationVerifiedSolution<RS>)guarantor.getSolution()).getX();
this.inverseOfbMat = new DenseMatrix<>(rowSize, columnSize, this.bMat.getDenseOrDiagonal(), unit);
this.inverseOfbMat.setElements((RM)x.getMiddle());
}
/**
* Compute a inverse matrix of X based on reliable computing
*
* @param currentPt information of current point
*/
private void getInverseMatrixOfXWithVerifier(Solution<RS,RM,CS,CM> currentPt) {
RoundModeManager manager = RoundModeManager.getManager();
manager.add(new FPURoundModeSelector());
manager.add(new MPFloatRoundModeSelector());
RS[][] matrixA = currentPt.getZmat().getArrayOfElements();
RM a = matrixA[0][0].createGrid(matrixA);
//NumericalMatrix<RS> a = new BaseNumericalMatrix<>(matrixA);
RS unit = this.bMat.getElementUnit();
int rowSize = currentPt.getZmat().getRowSize();
int columnSize = currentPt.getZmat().getColumnSize();
RM b = unit.createUnitGrid(rowSize, columnSize);
MPFloatIntlabMethod verifier = new MPFloatIntlabMethod((MPFloatMatrix)a,(MPFloatMatrix)b);
Guarantor<MPFloatIntlabMethod> guarantor = new Guarantor<>(verifier);
guarantor.solveWithEffectiveDigits(10);
MPFloatIntlabMethod guaranteedVerifier = guarantor.getGuranteedVerifier();
MPFloatIntervalMatrix x = guaranteedVerifier.getSolution();
//IntervalMatrix<RS> x = ((LinearEquationVerifiedSolution<RS>)guarantor.getSolution()).getX();
this.invzMat.setElements((RM)x.getMiddle());
}
/**
* @param fpout write a file
*/
public void display(File fpout) {
if (fpout == null) {
return;
}
try (PrintStream ps = new PrintStream(new FileOutputStream(fpout))) {
Tools.message(this.toString(), ps);
ps.close();
} catch (FileNotFoundException e) {
throw new RuntimeException(e);
}
}
/**
*
*/
public void display() {
Tools.message(this.toString());
}
/**
* @see java.lang.Object#toString()
*/
@Override
public String toString() {
StringBuffer buff = new StringBuffer();
buff.append("rNewton.DxMat = \n"); //$NON-NLS-1$
buff.append(this.DxMat.toString());
buff.append("rNewton.DyVec = \n"); //$NON-NLS-1$
buff.append(this.DyVec.toString());
buff.append("rNewton.DzMat = \n"); //$NON-NLS-1$
buff.append(this.DzMat.toString());
return buff.toString();
}
}