ProblemCreaterOnMemory.java
package org.mklab.sdpj.problem;
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.algorithm.Parameter;
import org.mklab.sdpj.algorithm.Solution;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.SparseDenseDiagonal;
import org.mklab.sdpj.datastructure.SparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
/**
* class <code>ProblemCreaterOnMemory</code> is defined for build a SDP problem,
* which can be computed correctly by solver <code>SDPJ</code>. With class
* <code>ProblemCreaterOnMemory</code>, it can concern to design how to deal with
* inputting data.
*
* @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 ProblemCreaterOnMemory<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>> {
/**a setting parameter of computing*/
private Parameter<RS,RM,CS,CM> parameter;
/**dimension of vector <code>x</code>*/
private int m;
/**the number of block of block-matrices array*/
private int nBlock;
/**a integer array save the block information of block-matrices array*/
private int[] blockStruct;
/**vector <code>b</code> := <code>c</code>*/
private Vector<RS,RM,CS,CM> b;
/**sparse block-matrices array <code>C</code> := <code>F0</code>*/
private BlockSparseMatrix<RS,RM,CS,CM> C;
/**sparse block-matrices array <code>A</code> := <code>Fi</code>*/
private BlockSparseMatrix<RS,RM,CS,CM>[] A;
/**the number of block matrix <code>C</code> := <code>F0</code>*/
private int[] CNonZeroCount;
/**the number of block-matrices array <code>A</code> := <code>Fi</code>*/
private int[] ANonZeroCount;
/**initial priaml-dual solution*/
private Solution<RS,RM,CS,CM> initPt;
/**current primal-dual solution*/
private Solution<RS,RM,CS,CM> currentPt;
/**block-array matrices <code>Fi</code>*/
private RS[][][][] F;
/**vector <code>x</code>*/
private RS[] vec;
/** */
private RS unit;
/**
* a class constructor, make and initialize a builded-SDP problem.
* <br>新しく生成された<code>ProblemCreaterOnMemory</code>オブジェクトを初期化します。
*
* @param F block-matrices array <code>Fi</code>
* @param vec vector <code>x</code>
* @param m dimension of vector <code>x</code>
* @param nBlock the number of block of block-matrices array
* @param blockStruct a integer array save the block information of block-matrices array
* @param parameter a setting parameter of computing
*/
public ProblemCreaterOnMemory(RS[][][][] F, RS[] vec, int m, int nBlock, int[] blockStruct, Parameter<RS,RM,CS,CM> parameter) {
this.unit = vec[0].create(1);
this.F = F;
this.setVec(vec);
this.m = m;
this.nBlock = nBlock;
this.blockStruct = blockStruct;
if (parameter != null) {
this.parameter = parameter;
} else {
this.parameter = new Parameter<>(this.unit);
this.parameter.setDefaultParameter();
}
this.CNonZeroCount = new int[nBlock];
this.ANonZeroCount = new int[nBlock * m];
this.A = new BlockSparseMatrix[m];
for (int i = 0; i < m; i++) {
this.A[i] = new BlockSparseMatrix<>();
}
this.b = new Vector<>(m, this.unit.create(0));
this.C = new BlockSparseMatrix<>(this.nBlock, this.blockStruct);
}
/**
* make a SDP problem that can be computed correctly by solver.
* <br>与えられたデータから問題を作成します。
*
* @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 F block-matrices array <code>Fi</code>
* @param vec vector <code>x</code>
* @param m dimension of vector <code>x</code>
// * @param n dimension of symmetric matrix <code>F</code> etc.
* @param nBlock the number of block of block-matrices array
* @param blockStruct a integer array save the block information of block-matrices array
* @param parameter a setting parameter of computing
* @return a SDP problem can be computed by solver
*/
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>> Problem<RS,RM,CS,CM> create(RS[][][][] F, RS[] vec, int m, int nBlock, int[] blockStruct, Parameter<RS,RM,CS,CM> parameter) {
//final RS unit = vec[0].create(1);
//Tools.setUnitNumber(unit);
ProblemCreaterOnMemory<RS,RM,CS,CM> creater = new ProblemCreaterOnMemory<>(F, vec, m, nBlock, blockStruct, parameter);
creater.b.setElements(vec);
creater.countNonzeroNumber();
creater.initializeCandA();
creater.setDataCandA(F);
creater.changeToDense();
Problem<RS,RM,CS,CM> problem = creater.getProblem();
return problem;
}
/**
* get a SDP problem that can be computed by solver.
*
* @return a builded-SDP problem
*/
public Problem<RS,RM,CS,CM> getProblem() {
return new Problem<>(this.m, this.nBlock, this.blockStruct, this.A, this.b, this.C, this.parameter, this.currentPt, this.initPt);
}
/**
* initialize matrix <code>C</code> and <code>A</code> with zero.
* <blockquote><tt>
* Note:<br>
* matrix <code>C</code> := matrix <code>F0</code><br>
* matrix <code>A</code> := matrix <code>Fi</code><br>
* </tt></blockquote>
*
* @param F block-matrices array
*/
private void setDataCandA(RS[][][][] F) {
for (int i = 0; i < this.nBlock; ++i) {
SparseMatrix<RS,RM,CS,CM> block = this.C.getBlock(i);
if (block.isSparce()) {
block.nonZeroCount = 0;
block.nonZeroEffect = 0;
}
}
for (int i = 0; i < this.m; ++i) {
for (int j = 0; j < this.nBlock; ++j) {
SparseMatrix<RS,RM,CS,CM> block = this.A[i].getBlock(j);
if (block.isSparce()) {
block.nonZeroCount = 0;
block.nonZeroEffect = 0;
}
}
}
// int counter=0;
for (int k = 0; k < this.nBlock; ++k) {
int size = this.blockStruct[k];
if (size > 0) {
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
RS tmp = F[0][k][i][j];
// counter++;
if (i <= j && !tmp.isZero()) {
SparseMatrix<RS,RM,CS,CM> target = this.C.getBlock(k);
int count = target.nonZeroCount;
target.rowIndex[count] = i;
target.columnIndex[count] = j;
// C must be opposite sign.
target.sparseElements[count] = tmp.unaryMinus();
target.nonZeroCount++;
if (i == j) {
target.nonZeroEffect++;
} else {
target.nonZeroEffect += 2;
}
}
}
}
} else {
for (int j = 0; j < -size; ++j) {
RS tmp = F[0][k][j][j];
// C must be opposite sign.
this.C.getBlock(k).diagonalElements[j] = tmp.unaryMinus();
}
}
}
for (int k = 0; k < this.m; ++k) {
// counter=0;
for (int p = 0; p < this.nBlock; ++p) {
int size = this.blockStruct[p];
if (size > 0) {
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
RS tmp = F[k + 1][p][i][j];
// counter++;
if (i <= j && !tmp.isZero()) {
SparseMatrix<RS,RM,CS,CM> target = this.A[k].getBlock(p);
int count = target.nonZeroCount;
target.rowIndex[count] = i;
target.columnIndex[count] = j;
target.sparseElements[count] = tmp;
target.nonZeroCount++;
if (i == j) {
target.nonZeroEffect++;
} else {
target.nonZeroEffect += 2;
}
}
}
}
} else {
for (int j = 0; j < -size; ++j) {
RS tmp = F[k + 1][p][j][j];
// counter++;
this.A[k].getBlock(p).diagonalElements[j] = tmp;
}
}
}
}
}
/**
* initialize matrix <code>C</code> and <code>A</code>.
* <blockquote><tt>
* Note:<br>
* matrix <code>C</code> := matrix <code>F0</code><br>
* matrix <code>A</code> := matrix <code>Fi</code><br>
* </tt></blockquote>
*/
private void initializeCandA() {
for (int i = 0; i < this.nBlock; ++i) {
int size = this.blockStruct[i];
if (size > 0) {
this.C.setBlock(i, new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, this.CNonZeroCount[i], this.unit));
} else {
this.C.setBlock(i, new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit));
}
}
for (int i = 0; i < this.m; ++i) {
this.A[i] = new BlockSparseMatrix<>(this.nBlock, this.blockStruct);
for (int j = 0; j < this.nBlock; ++j) {
int size = this.blockStruct[j];
if (size > 0) {
this.A[i].setBlock(j, new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, this.ANonZeroCount[i * this.nBlock + j], this.unit));
} else {
this.A[i].setBlock(j, new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit));
}
}
}
}
/**
* count the number of nonzero elements of matrix.
*/
private void countNonzeroNumber() {
for (int i = 0; i < this.nBlock; i++) {
this.CNonZeroCount[i] = 0;
for (int j = 0; j < this.m; j++) {
this.ANonZeroCount[j * this.nBlock + i] = 0;
}
}
// k==0
for (int l = 0; l < this.nBlock; ++l) {
int size = this.blockStruct[l];
if (size > 0) {
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
RS scalar = this.F[0][l][i][j];
if (i <= j && !scalar.isZero()) {
this.CNonZeroCount[l]++;
}
}
}
} else {
for (int j = 0; j < -size; ++j) {
this.CNonZeroCount[l]++;
}
}
}
for (int k = 0; k < this.m; ++k) {
for (int p = 0; p < this.nBlock; ++p) {
int size = this.blockStruct[p];
if (size > 0) {
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
RS scalar = this.F[k + 1][p][i][j];
if (i <= j && !scalar.isZero()) {
this.ANonZeroCount[k * this.nBlock + p]++;
}
}
}
} else {
for (int j = 0; j < -size; ++j) {
this.ANonZeroCount[k * this.nBlock + p]++;
}
}
}
}
}
/**
* change sparse matrix into dense matrix.
*/
private void changeToDense() {
// if possible , change C and A to Dense
this.C.changeToDense();
for (int i = 0; i < this.m; ++i) {
this.A[i].changeToDense();
}
}
/**
* @param vec vecを設定します。
*/
public void setVec(RS[] vec) {
this.vec = vec;
}
/**
* @return vecを返します。
*/
public RS[] getVec() {
return this.vec;
}
}