java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/BorderedBandMatrix.java 2014-09-16 21:36:59 UTC (rev 1030)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/BorderedBandMatrix.java 2014-09-16 23:41:59 UTC (rev 1031)
@@ -3,6 +3,7 @@
import static java.lang.Math.abs;
import static java.lang.Math.min;
import static java.lang.Math.max;
+import java.util.List;
/**
*
@@ -37,8 +38,8 @@
private int numBand; ///< Band width
private int numCol; ///< Band matrix size
private VSymMatrix theBorder = new VSymMatrix(0); ///< Border part
- private VMatrix theMixed = new VMatrix(0,0); ///< Mixed part
- private VMatrix theBand = new VMatrix(0,0); ///< Band part
+ private VMatrix theMixed = new VMatrix(0, 0); ///< Mixed part
+ private VMatrix theBand = new VMatrix(0, 0); ///< Band part
/// Resize bordered band matrix.
/**
@@ -65,26 +66,21 @@
* of rows/colums to be used \param aVector [in] Vector
*/
public void addBlockMatrix(double aWeight,
- int[] anIndex,
- double[] aVector)
+ int[] anIndex,
+ double[] aVector)
{
int nBorder = numBorder;
- for (int i = 0; i < anIndex.length; ++i)
- {
+ for (int i = 0; i < anIndex.length; ++i) {
int iIndex = (anIndex)[i] - 1; // anIndex has to be sorted
- for (int j = 0; j <= i; ++j)
- {
+ for (int j = 0; j <= i; ++j) {
int jIndex = (anIndex)[j] - 1;
- if (iIndex < nBorder)
- {
+ if (iIndex < nBorder) {
theBorder.addTo(iIndex, jIndex, aVector[i] * aWeight
* aVector[j]);
- } else if (jIndex < nBorder)
- {
+ } else if (jIndex < nBorder) {
theMixed.addTo(jIndex, iIndex - nBorder, aVector[i] * aWeight
* aVector[j]);
- } else
- {
+ } else {
int nBand = iIndex - jIndex;
theBand.addTo(nBand, jIndex - nBorder, aVector[i] * aWeight
* aVector[j]);
@@ -94,34 +90,36 @@
}
}
-///// Retrieve symmetric block matrix.
-///**
-// * Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).
-// * \param anIndex [in] List of rows/colums to be used
-// */
-//TMatrixDSym BorderedBandMatrix::getBlockMatrix(
-// const std::vector<unsigned int> anIndex) const {
-//
-// TMatrixDSym aMatrix(anIndex.size());
-// int nBorder = numBorder;
-// for (unsigned int i = 0; i < anIndex.size(); ++i) {
-// int iIndex = anIndex[i] - 1; // anIndex has to be sorted
-// for (unsigned int j = 0; j <= i; ++j) {
-// int jIndex = anIndex[j] - 1;
-// if (iIndex < nBorder) {
-// aMatrix(i, j) = theBorder(iIndex, jIndex); // border part of inverse
-// } else if (jIndex < nBorder) {
-// aMatrix(i, j) = -theMixed(jIndex, iIndex - nBorder); // mixed part of inverse
-// } else {
-// unsigned int nBand = iIndex - jIndex;
-// aMatrix(i, j) = theBand(nBand, jIndex - nBorder); // band part of inverse
-// }
-// aMatrix(j, i) = aMatrix(i, j);
-// }
-// }
-// return aMatrix;
-//}
+/// Retrieve symmetric block matrix.
+ /**
+ * Get (compressed) block from bordered band matrix: aMatrix(i,j) =
+ * BBmatrix(anIndex(i),anIndex(j)). \param anIndex [in] List of rows/colums
+ * to be used
+ */
+ public SymMatrix getBlockMatrix(List<Integer> anIndex)
+ {
+
+ SymMatrix aMatrix = new SymMatrix(anIndex.size());
+ int nBorder = numBorder;
+ for (int i = 0; i < anIndex.size(); ++i) {
+ int iIndex = anIndex.get(i) - 1; // anIndex has to be sorted
+ for (int j = 0; j <= i; ++j) {
+ int jIndex = anIndex.get(j) - 1;
+ if (iIndex < nBorder) {
+ aMatrix.set(i, j, theBorder.get(iIndex, jIndex)); // border part of inverse
+ } else if (jIndex < nBorder) {
+ aMatrix.set(i, j, -theMixed.get(jIndex, iIndex - nBorder)); // mixed part of inverse
+ } else {
+ int nBand = iIndex - jIndex;
+ aMatrix.set(i, j, theBand.get(nBand, jIndex - nBorder)); // band part of inverse
+ }
+ aMatrix.set(j, i, aMatrix.get(i, j));
+ }
+ }
+ return aMatrix;
+ }
/// Solve linear equation system, partially calculate inverse.
+
/**
* Solve linear equation A*x=b system with bordered band matrix A, calculate
* bordered band part of inverse of A. Use decomposition in border and band
@@ -146,38 +144,40 @@
* \param [in] aRightHandSide Right hand side (vector) 'b' of A*x=b \param
* [out] aSolution Solution (vector) x of A*x=b
*/
-public void solveAndInvertBorderedBand(
- VVector aRightHandSide, VVector aSolution) {
+ public void solveAndInvertBorderedBand(
+ VVector aRightHandSide, VVector aSolution)
+ {
- // decompose band
- decomposeBand();
- // invert band
- VMatrix inverseBand = invertBand();
- if (numBorder > 0) { // need to use block matrix decomposition to solve
- // solve for mixed part
- VMatrix auxMat = solveBand(theMixed); // = Xt
- VMatrix auxMatT = auxMat.transpose(); // = X
- // solve for border part
- VVector auxVec = aRightHandSide.getVec(numBorder,0).minus(
- auxMat.times( aRightHandSide.getVec(numCol, numBorder) ) ); // = b1 - Xt*b2
- VSymMatrix inverseBorder = theBorder.minus(theMixed.times(auxMatT) );
- inverseBorder.invert(); // = E
- VVector borderSolution = inverseBorder.times(auxVec); // = x1
- // solve for band part
- VVector bandSolution = solveBand(
- aRightHandSide.getVec(numCol, numBorder)); // = x
- aSolution.putVec(borderSolution,0);
- aSolution.putVec(bandSolution.minus(auxMatT.times(borderSolution)), numBorder); // = x2
- // parts of inverse
- theBorder = inverseBorder; // E
- theMixed = inverseBorder.times(auxMat); // E*Xt (-mixed part of inverse) !!!
- theBand = inverseBand.plus(bandOfAVAT(auxMatT, inverseBorder)); // band(D^-1 + X*E*Xt)
- } else {
- aSolution.putVec(solveBand(aRightHandSide),0);
- theBand = inverseBand;
- }
-}
+ // decompose band
+ decomposeBand();
+ // invert band
+ VMatrix inverseBand = invertBand();
+ if (numBorder > 0) { // need to use block matrix decomposition to solve
+ // solve for mixed part
+ VMatrix auxMat = solveBand(theMixed); // = Xt
+ VMatrix auxMatT = auxMat.transpose(); // = X
+ // solve for border part
+ VVector auxVec = aRightHandSide.getVec(numBorder, 0).minus(
+ auxMat.times(aRightHandSide.getVec(numCol, numBorder))); // = b1 - Xt*b2
+ VSymMatrix inverseBorder = theBorder.minus(theMixed.times(auxMatT));
+ inverseBorder.invert(); // = E
+ VVector borderSolution = inverseBorder.times(auxVec); // = x1
+ // solve for band part
+ VVector bandSolution = solveBand(
+ aRightHandSide.getVec(numCol, numBorder)); // = x
+ aSolution.putVec(borderSolution, 0);
+ aSolution.putVec(bandSolution.minus(auxMatT.times(borderSolution)), numBorder); // = x2
+ // parts of inverse
+ theBorder = inverseBorder; // E
+ theMixed = inverseBorder.times(auxMat); // E*Xt (-mixed part of inverse) !!!
+ theBand = inverseBand.plus(bandOfAVAT(auxMatT, inverseBorder)); // band(D^-1 + X*E*Xt)
+ } else {
+ aSolution.putVec(solveBand(aRightHandSide), 0);
+ theBand = inverseBand;
+ }
+ }
/// Print bordered band matrix.
+
public void printMatrix()
{
System.out.println("Border part ");
@@ -204,29 +204,22 @@
int nRow = numBand + 1;
int nCol = numCol;
VVector auxVec = new VVector(nCol);
- for (int i = 0; i < nCol; ++i)
- {
+ for (int i = 0; i < nCol; ++i) {
auxVec.set(i, theBand.get(0, i) * 16.0); // save diagonal elements
}
- for (int i = 0; i < nCol; ++i)
- {
- if ((theBand.get(0, i) + auxVec.get(i)) != theBand.get(0, i))
- {
+ for (int i = 0; i < nCol; ++i) {
+ if ((theBand.get(0, i) + auxVec.get(i)) != theBand.get(0, i)) {
theBand.set(0, i, 1.0 / theBand.get(0, i));
- if (theBand.get(0, i) < 0.)
- {
+ if (theBand.get(0, i) < 0.) {
throw new RuntimeException("BorderedBandMatrix decomposeBand not positive definite");
}
- } else
- {
+ } else {
theBand.set(0, i, 0.0);
throw new RuntimeException("BorderedBandMatrix decomposeBand singular");
}
- for (int j = 1; j < min(nRow, nCol - i); ++j)
- {
+ for (int j = 1; j < min(nRow, nCol - i); ++j) {
double rxw = theBand.get(j, i) * theBand.get(0, i);
- for (int k = 0; k < min(nRow, nCol - i) - j; ++k)
- {
+ for (int k = 0; k < min(nRow, nCol - i) - j; ++k) {
theBand.subFrom(k, i + j, theBand.get(k + j, i) * rxw);
}
theBand.set(j, i, rxw);
@@ -245,13 +238,10 @@
int nCol = numCol;
VMatrix inverseBand = new VMatrix(nRow, nCol);
- for (int i = nCol - 1; i >= 0; i--)
- {
+ for (int i = nCol - 1; i >= 0; i--) {
double rxw = theBand.get(0, i);
- for (int j = i; j >= max(0, i - nRow + 1); j--)
- {
- for (int k = j + 1; k < min(nCol, j + nRow); ++k)
- {
+ for (int j = i; j >= max(0, i - nRow + 1); j--) {
+ for (int k = j + 1; k < min(nCol, j + nRow); ++k) {
rxw -= inverseBand.get(abs(i - k), min(i, k))
* theBand.get(k - j, j);
}
@@ -277,16 +267,14 @@
VVector aSolution = new VVector(aRightHandSide);
for (int i = 0; i < nCol; ++i) // forward substitution
{
- for (int j = 1; j < min(nRow, nCol - i); ++j)
- {
+ for (int j = 1; j < min(nRow, nCol - i); ++j) {
aSolution.subFrom(j + i, theBand.get(j, i) * aSolution.get(i));
}
}
for (int i = nCol - 1; i >= 0; i--) // backward substitution
{
double rxw = theBand.get(0, i) * aSolution.get(i);
- for (int j = 1; j < min(nRow, nCol - i); ++j)
- {
+ for (int j = 1; j < min(nRow, nCol - i); ++j) {
rxw -= theBand.get(j, i) * aSolution.get(j + i);
}
aSolution.set(i, rxw);
@@ -306,12 +294,10 @@
int nRow = theBand.getNumRows();
int nCol = theBand.getNumCols();
VMatrix aSolution = new VMatrix(aRightHandSide);
- for (int iBorder = 0; iBorder < numBorder; iBorder++)
- {
+ for (int iBorder = 0; iBorder < numBorder; iBorder++) {
for (int i = 0; i < nCol; ++i) // forward substitution
{
- for (int j = 1; j < min(nRow, nCol - i); ++j)
- {
+ for (int j = 1; j < min(nRow, nCol - i); ++j) {
aSolution.subFrom(iBorder, j + i, theBand.get(j, i)
* aSolution.get(iBorder, i));
}
@@ -319,8 +305,7 @@
for (int i = nCol - 1; i >= 0; i--) // backward substitution
{
double rxw = theBand.get(0, i) * aSolution.get(iBorder, i);
- for (int j = 1; j < min(nRow, nCol - i); ++j)
- {
+ for (int j = 1; j < min(nRow, nCol - i); ++j) {
rxw -= theBand.get(j, i) * aSolution.get(iBorder, j + i);
}
aSolution.set(iBorder, i, rxw);
@@ -329,32 +314,31 @@
return aSolution;
}
-
/// Calculate band part of: 'anArray * aSymArray * anArray.T'.
-/**
- * \return Band part of product
- */
-private VMatrix bandOfAVAT( VMatrix anArray,
- VSymMatrix aSymArray) {
- int nBand = numBand;
- int nCol = numCol;
- int nBorder = numBorder;
- double sum;
- VMatrix aBand = new VMatrix((nBand + 1), nCol);
- for (int i = 0; i < nCol; ++i) {
- for (int j = max(0, i - nBand); j <= i; ++j) {
- sum = 0.;
- for (int l = 0; l < nBorder; ++l) { // diagonal
- sum += anArray.get(i, l) * aSymArray.get(l, l) * anArray.get(j, l);
- for (int k = 0; k < l; ++k) { // off diagonal
- sum += anArray.get(i, l) * aSymArray.get(l, k) * anArray.get(j, k)
- + anArray.get(i, k) * aSymArray.get(l, k) * anArray.get(j, l);
- }
- }
- aBand.set(i - j, j, sum);
- }
- }
- return aBand;
-}
-
-}
+ /**
+ * \return Band part of product
+ */
+ private VMatrix bandOfAVAT(VMatrix anArray,
+ VSymMatrix aSymArray)
+ {
+ int nBand = numBand;
+ int nCol = numCol;
+ int nBorder = numBorder;
+ double sum;
+ VMatrix aBand = new VMatrix((nBand + 1), nCol);
+ for (int i = 0; i < nCol; ++i) {
+ for (int j = max(0, i - nBand); j <= i; ++j) {
+ sum = 0.;
+ for (int l = 0; l < nBorder; ++l) { // diagonal
+ sum += anArray.get(i, l) * aSymArray.get(l, l) * anArray.get(j, l);
+ for (int k = 0; k < l; ++k) { // off diagonal
+ sum += anArray.get(i, l) * aSymArray.get(l, k) * anArray.get(j, k)
+ + anArray.get(i, k) * aSymArray.get(l, k) * anArray.get(j, l);
+ }
+ }
+ aBand.set(i - j, j, sum);
+ }
+ }
+ return aBand;
+ }
+}
\ No newline at end of file
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/VMatrix.java 2014-09-16 21:36:59 UTC (rev 1030)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/VMatrix.java 2014-09-16 23:41:59 UTC (rev 1031)
@@ -23,33 +23,31 @@
numRows = nRows;
numCols = nCols;
theVec = new ArrayList(nRows * nCols);
- for (int i = 0; i < numRows * nCols; ++i)
- {
+ for (int i = 0; i < numRows * nCols; ++i) {
theVec.add(0.);
}
}
-
-public VMatrix(VMatrix m)
+
+ public VMatrix(VMatrix m)
{
numRows = m.numRows;
numCols = m.numCols;
theVec = new ArrayList(numRows * numCols);
- for (int i = 0; i < numRows; ++i)
- {
- for(int j=0; j<numCols; ++j)
- {
- theVec.set(numCols * i + j, m.get(i, j));
+ for (int i = 0; i < numRows * numCols; ++i) {
+ theVec.add(0.);
+ }
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < numCols; ++j) {
+ theVec.set(numCols * i + j, m.get(i, j));
}
}
- }
+ }
public VMatrix copy()
{
VMatrix aResult = new VMatrix(numRows, numCols);
- for (int i = 0; i < numRows; ++i)
- {
- for (int j = 0; j < numCols; ++j)
- {
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < numCols; ++j) {
{
aResult.set(i, j, get(i, j));
}
@@ -69,8 +67,7 @@
numRows = nRows;
numCols = nCols;
theVec = new ArrayList(numRows * nCols);
- for (int i = 0; i < numRows * nCols; ++i)
- {
+ for (int i = 0; i < numRows * nCols; ++i) {
theVec.add(0.);
}
}
@@ -82,10 +79,8 @@
public VMatrix transpose()
{
VMatrix aResult = new VMatrix(numCols, numRows);
- for (int i = 0; i < numRows; ++i)
- {
- for (int j = 0; j < numCols; ++j)
- {
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < numCols; ++j) {
//System.out.println("row: "+i+" col: "+j+" val: "+theVec.get(numCols * i + j));
aResult.set(j, i, (double) theVec.get(numCols * i + j));
}
@@ -97,14 +92,15 @@
{
theVec.set(numCols * row + col, val);
}
-
- public void addTo(int row, int col, double val)
+
+ public void addTo(int row, int col, double val)
{
- theVec.set(numCols * row + col, (double) theVec.get(numCols * row + col)+val);
+ theVec.set(numCols * row + col, (double) theVec.get(numCols * row + col) + val);
}
- public void subFrom(int row, int col, double val)
+
+ public void subFrom(int row, int col, double val)
{
- theVec.set(numCols * row + col, (double) theVec.get(numCols * row + col)-val);
+ theVec.set(numCols * row + col, (double) theVec.get(numCols * row + col) - val);
}
public double get(int row, int col)
@@ -134,11 +130,9 @@
VVector times(VVector aVector)
{
VVector aResult = new VVector(numRows);
- for (int i = 0; i < numRows; ++i)
- {
+ for (int i = 0; i < numRows; ++i) {
double sum = 0.0;
- for (int j = 0; j < numCols; ++j)
- {
+ for (int j = 0; j < numCols; ++j) {
sum += (double) theVec.get(numCols * i + j) * aVector.get(j);
}
aResult.set(i, sum);
@@ -151,13 +145,10 @@
{
VMatrix aResult = new VMatrix(numRows, aMatrix.numCols);
- for (int i = 0; i < numRows; ++i)
- {
- for (int j = 0; j < aMatrix.numCols; ++j)
- {
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < aMatrix.numCols; ++j) {
double sum = 0.0;
- for (int k = 0; k < numCols; ++k)
- {
+ for (int k = 0; k < numCols; ++k) {
sum += (double) theVec.get(numCols * i + k) * aMatrix.get(k, j);
}
aResult.set(i, j, sum);
@@ -170,10 +161,8 @@
VMatrix plus(VMatrix aMatrix)
{
VMatrix aResult = new VMatrix(numRows, numCols);
- for (int i = 0; i < numRows; ++i)
- {
- for (int j = 0; j < numCols; ++j)
- {
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < numCols; ++j) {
aResult.set(i, j, (double) theVec.get(numCols * i + j) + (double) aMatrix.get(i, j));
}
}
@@ -184,12 +173,9 @@
public void print()
{
System.out.println(" VMatrix: " + numRows + "*" + numCols);
- for (int i = 0; i < numRows; ++i)
- {
- for (int j = 0; j < numCols; ++j)
- {
- if (j % 5 == 0)
- {
+ for (int i = 0; i < numRows; ++i) {
+ for (int j = 0; j < numCols; ++j) {
+ if (j % 5 == 0) {
System.out.format("%n%4d " + "," + "%4d" + " - " + "%4d" + " : ", i, j, min(j + 4, numCols));
}
System.out.format("%13f", theVec.get(numCols * i + j));
@@ -197,5 +183,4 @@
}
System.out.print("\n\n\n");
}
-
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/VVector.java 2014-09-16 21:36:59 UTC (rev 1030)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/matrix/VVector.java 2014-09-16 23:41:59 UTC (rev 1031)
@@ -21,8 +21,7 @@
{
numRows = nRows;
theVec = new ArrayList(nRows);
- for (int i = 0; i < numRows; ++i)
- {
+ for (int i = 0; i < numRows; ++i) {
theVec.add(0.);
}
}
@@ -44,12 +43,12 @@
{
theVec.set(row, val);
}
-
+
public void addTo(int row, double val)
{
theVec.set(row, (double) theVec.get(row) + val);
}
-
+
public void subFrom(int row, double val)
{
theVec.set(row, (double) theVec.get(row) - val);
@@ -67,7 +66,7 @@
*/
public VVector getVec(int len, int start)
{
- return new VVector(theVec.subList(start, start + len - 1));
+ return new VVector(theVec.subList(start, start + len));
}
@@ -93,10 +92,8 @@
public void print()
{
System.out.println(theVec.size());
- for (int i = 0; i < numRows; ++i)
- {
- if (i % 5 == 0)
- {
+ for (int i = 0; i < numRows; ++i) {
+ if (i % 5 == 0) {
System.out.format("%n%4d " + " - " + "%4d" + " : ", i, min(i + 4, numRows));
}
System.out.format("%13f", theVec.get(i));
@@ -108,8 +105,7 @@
VVector minus(VVector aVector)
{
VVector aResult = new VVector(numRows);
- for (int i = 0; i < numRows; ++i)
- {
+ for (int i = 0; i < numRows; ++i) {
aResult.set(i, (double) theVec.get(i) - aVector.get(i));
}
return aResult;
@@ -119,8 +115,7 @@
public VVector copy()
{
VVector aResult = new VVector(numRows);
- for (int i = 0; i < numRows; ++i)
- {
+ for (int i = 0; i < numRows; ++i) {
aResult.set(i, get(i));
}
return aResult;