Belle II Software  release-05-01-25
BorderedBandMatrix.cc
Go to the documentation of this file.
1 /*
2  * BorderedBandMatrix.cpp
3  *
4  * Created on: Aug 14, 2011
5  * Author: kleinwrt
6  */
7 
30 #include "BorderedBandMatrix.h"
31 
33 namespace gbl {
34 
36 BorderedBandMatrix::BorderedBandMatrix() : numSize(0), numBorder(0), numBand(0), numCol(0) {
37 }
38 
39 BorderedBandMatrix::~BorderedBandMatrix() {
40 }
41 
43 
48 void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
49  unsigned int nBand) {
50  numSize = nSize;
51  numBorder = nBorder;
52  numCol = nSize - nBorder;
53  numBand = 0;
56  theBand.resize((nBand + 1), numCol);
57 }
58 
60 
69  const std::vector<unsigned int>* anIndex,
70  const std::vector<double>* aVector) {
71  int nBorder = numBorder;
72  for (unsigned int i = 0; i < anIndex->size(); ++i) {
73  int iIndex = (*anIndex)[i] - 1; // anIndex has to be sorted
74  for (unsigned int j = 0; j <= i; ++j) {
75  int jIndex = (*anIndex)[j] - 1;
76  if (iIndex < nBorder) {
77  theBorder(iIndex, jIndex) += (*aVector)[i] * aWeight
78  * (*aVector)[j];
79  } else if (jIndex < nBorder) {
80  theMixed(jIndex, iIndex - nBorder) += (*aVector)[i] * aWeight
81  * (*aVector)[j];
82  } else {
83  unsigned int nBand = iIndex - jIndex;
84  theBand(nBand, jIndex - nBorder) += (*aVector)[i] * aWeight
85  * (*aVector)[j];
86  numBand = std::max(numBand, nBand); // update band width
87  }
88  }
89  }
90 }
91 
93 
98  const std::vector<unsigned int> anIndex) const {
99 
100  TMatrixDSym aMatrix(anIndex.size());
101  int nBorder = numBorder;
102  for (unsigned int i = 0; i < anIndex.size(); ++i) {
103  int iIndex = anIndex[i] - 1; // anIndex has to be sorted
104  for (unsigned int j = 0; j <= i; ++j) {
105  int jIndex = anIndex[j] - 1;
106  if (iIndex < nBorder) {
107  aMatrix(i, j) = theBorder(iIndex, jIndex); // border part of inverse
108  } else if (jIndex < nBorder) {
109  aMatrix(i, j) = -theMixed(jIndex, iIndex - nBorder); // mixed part of inverse
110  } else {
111  unsigned int nBand = iIndex - jIndex;
112  aMatrix(i, j) = theBand(nBand, jIndex - nBorder); // band part of inverse
113  }
114  aMatrix(j, i) = aMatrix(i, j);
115  }
116  }
117  return aMatrix;
118 }
119 
121 
148  const VVector &aRightHandSide, VVector &aSolution) {
149 
150  // decompose band
151  decomposeBand();
152  // invert band
153  VMatrix inverseBand = invertBand();
154  if (numBorder > 0) { // need to use block matrix decomposition to solve
155  // solve for mixed part
156  const VMatrix auxMat = solveBand(theMixed); // = Xt
157  const VMatrix auxMatT = auxMat.transpose(); // = X
158  // solve for border part
159  const VVector auxVec = aRightHandSide.getVec(numBorder)
160  - auxMat * aRightHandSide.getVec(numCol, numBorder); // = b1 - Xt*b2
161  VSymMatrix inverseBorder = theBorder - theMixed * auxMatT;
162  inverseBorder.invert(); // = E
163  const VVector borderSolution = inverseBorder * auxVec; // = x1
164  // solve for band part
165  const VVector bandSolution = solveBand(
166  aRightHandSide.getVec(numCol, numBorder)); // = x
167  aSolution.putVec(borderSolution);
168  aSolution.putVec(bandSolution - auxMatT * borderSolution, numBorder); // = x2
169  // parts of inverse
170  theBorder = inverseBorder; // E
171  theMixed = inverseBorder * auxMat; // E*Xt (-mixed part of inverse) !!!
172  theBand = inverseBand + bandOfAVAT(auxMatT, inverseBorder); // band(D^-1 + X*E*Xt)
173  } else {
174  aSolution.putVec(solveBand(aRightHandSide));
175  theBand = inverseBand;
176  }
177 }
178 
181  std::cout << "Border part " << std::endl;
182  theBorder.print();
183  std::cout << "Mixed part " << std::endl;
184  theMixed.print();
185  std::cout << "Band part " << std::endl;
186  theBand.print();
187 }
188 
189 /*============================================================================
190  from Dbandmatrix.F (MillePede-II by V. Blobel, Univ. Hamburg)
191  ============================================================================*/
193 
200 
201  int nRow = numBand + 1;
202  int nCol = numCol;
203  VVector auxVec(nCol);
204  for (int i = 0; i < nCol; ++i) {
205  auxVec(i) = theBand(0, i) * 16.0; // save diagonal elements
206  }
207  for (int i = 0; i < nCol; ++i) {
208  if ((theBand(0, i) + auxVec(i)) != theBand(0, i)) {
209  theBand(0, i) = 1.0 / theBand(0, i);
210  if (theBand(0, i) < 0.) {
211  throw 3; // not positive definite
212  }
213  } else {
214  theBand(0, i) = 0.0;
215  throw 2; // singular
216  }
217  for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
218  double rxw = theBand(j, i) * theBand(0, i);
219  for (int k = 0; k < std::min(nRow, nCol - i) - j; ++k) {
220  theBand(k, i + j) -= theBand(k + j, i) * rxw;
221  }
222  theBand(j, i) = rxw;
223  }
224  }
225 }
226 
228 
234 VVector BorderedBandMatrix::solveBand(const VVector &aRightHandSide) const {
235 
236  int nRow = theBand.getNumRows();
237  int nCol = theBand.getNumCols();
238  VVector aSolution(aRightHandSide);
239  for (int i = 0; i < nCol; ++i) // forward substitution
240  {
241  for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
242  aSolution(j + i) -= theBand(j, i) * aSolution(i);
243  }
244  }
245  for (int i = nCol - 1; i >= 0; i--) // backward substitution
246  {
247  double rxw = theBand(0, i) * aSolution(i);
248  for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
249  rxw -= theBand(j, i) * aSolution(j + i);
250  }
251  aSolution(i) = rxw;
252  }
253  return aSolution;
254 }
255 
257 
263 VMatrix BorderedBandMatrix::solveBand(const VMatrix &aRightHandSide) const {
264 
265  int nRow = theBand.getNumRows();
266  int nCol = theBand.getNumCols();
267  VMatrix aSolution(aRightHandSide);
268  for (unsigned int iBorder = 0; iBorder < numBorder; iBorder++) {
269  for (int i = 0; i < nCol; ++i) // forward substitution
270  {
271  for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
272  aSolution(iBorder, j + i) -= theBand(j, i)
273  * aSolution(iBorder, i);
274  }
275  }
276  for (int i = nCol - 1; i >= 0; i--) // backward substitution
277  {
278  double rxw = theBand(0, i) * aSolution(iBorder, i);
279  for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
280  rxw -= theBand(j, i) * aSolution(iBorder, j + i);
281  }
282  aSolution(iBorder, i) = rxw;
283  }
284  }
285  return aSolution;
286 }
287 
289 
293 
294  int nRow = numBand + 1;
295  int nCol = numCol;
296  VMatrix inverseBand(nRow, nCol);
297 
298  for (int i = nCol - 1; i >= 0; i--) {
299  double rxw = theBand(0, i);
300  for (int j = i; j >= std::max(0, i - nRow + 1); j--) {
301  for (int k = j + 1; k < std::min(nCol, j + nRow); ++k) {
302  rxw -= inverseBand(abs(i - k), std::min(i, k))
303  * theBand(k - j, j);
304  }
305  inverseBand(i - j, j) = rxw;
306  rxw = 0.;
307  }
308  }
309  return inverseBand;
310 }
311 
313 
317  const VSymMatrix &aSymArray) const {
318  int nBand = numBand;
319  int nCol = numCol;
320  int nBorder = numBorder;
321  double sum;
322  VMatrix aBand((nBand + 1), nCol);
323  for (int i = 0; i < nCol; ++i) {
324  for (int j = std::max(0, i - nBand); j <= i; ++j) {
325  sum = 0.;
326  for (int l = 0; l < nBorder; ++l) { // diagonal
327  sum += anArray(i, l) * aSymArray(l, l) * anArray(j, l);
328  for (int k = 0; k < l; ++k) { // off diagonal
329  sum += anArray(i, l) * aSymArray(l, k) * anArray(j, k)
330  + anArray(i, k) * aSymArray(l, k) * anArray(j, l);
331  }
332  }
333  aBand(i - j, j) = sum;
334  }
335  }
336  return aBand;
337 }
338 
339 }
gbl::VMatrix::getNumCols
unsigned int getNumCols() const
Get number of columns.
Definition: VMatrix.cc:87
gbl::VVector::putVec
void putVec(const VVector &aVector, unsigned int start=0)
Put part of vector.
Definition: VMatrix.cc:284
gbl::BorderedBandMatrix::theBorder
VSymMatrix theBorder
Border part.
Definition: BorderedBandMatrix.h:98
gbl::VMatrix::resize
void resize(const unsigned int nRows, const unsigned int nCols)
Resize Matrix.
Definition: VMatrix.cc:55
gbl::BorderedBandMatrix::numSize
unsigned int numSize
Matrix size.
Definition: BorderedBandMatrix.h:94
gbl::VSymMatrix::resize
void resize(const unsigned int nRows)
Resize symmetric matrix.
Definition: VMatrix.cc:175
gbl::BorderedBandMatrix::BorderedBandMatrix
BorderedBandMatrix()
Create bordered band matrix.
Definition: BorderedBandMatrix.cc:36
gbl::VVector
Simple Vector based on std::vector<double>
Definition: VMatrix.h:43
gbl::BorderedBandMatrix::solveBand
VVector solveBand(const VVector &aRightHandSide) const
Solve for band part.
Definition: BorderedBandMatrix.cc:234
gbl
Namespace for the general broken lines package.
Definition: BorderedBandMatrix.h:44
gbl::VMatrix::transpose
VMatrix transpose() const
Get transposed matrix.
Definition: VMatrix.cc:65
gbl::BorderedBandMatrix::bandOfAVAT
VMatrix bandOfAVAT(const VMatrix &anArray, const VSymMatrix &aSymArray) const
Calculate band part of: 'anArray * aSymArray * anArray.T'.
Definition: BorderedBandMatrix.cc:316
gbl::BorderedBandMatrix::theBand
VMatrix theBand
Band part.
Definition: BorderedBandMatrix.h:100
gbl::BorderedBandMatrix::printMatrix
void printMatrix() const
Print bordered band matrix.
Definition: BorderedBandMatrix.cc:180
BorderedBandMatrix.h
gbl::VMatrix::getNumRows
unsigned int getNumRows() const
Get number of rows.
Definition: VMatrix.cc:79
gbl::VMatrix::print
void print() const
Print matrix.
Definition: VMatrix.cc:92
gbl::VMatrix
Simple Matrix based on std::vector<double>
Definition: VMatrix.h:63
gbl::BorderedBandMatrix::numBorder
unsigned int numBorder
Border size.
Definition: BorderedBandMatrix.h:95
gbl::BorderedBandMatrix::addBlockMatrix
void addBlockMatrix(double aWeight, const std::vector< unsigned int > *anIndex, const std::vector< double > *aVector)
Add symmetric block matrix.
Definition: BorderedBandMatrix.cc:68
gbl::BorderedBandMatrix::solveAndInvertBorderedBand
void solveAndInvertBorderedBand(const VVector &aRightHandSide, VVector &aSolution)
Solve linear equation system, partially calculate inverse.
Definition: BorderedBandMatrix.cc:147
gbl::BorderedBandMatrix::decomposeBand
void decomposeBand()
(root free) Cholesky decomposition of band part: C=LDL^T
Definition: BorderedBandMatrix.cc:199
gbl::VSymMatrix::invert
unsigned int invert()
Matrix inversion.
Definition: VMatrix.cc:348
gbl::BorderedBandMatrix::numCol
unsigned int numCol
Band matrix size.
Definition: BorderedBandMatrix.h:97
gbl::BorderedBandMatrix::resize
void resize(unsigned int nSize, unsigned int nBorder=1, unsigned int nBand=5)
Resize bordered band matrix.
Definition: BorderedBandMatrix.cc:48
gbl::VSymMatrix::print
void print() const
Print matrix.
Definition: VMatrix.cc:189
gbl::BorderedBandMatrix::getBlockMatrix
TMatrixDSym getBlockMatrix(const std::vector< unsigned int > anIndex) const
Retrieve symmetric block matrix.
Definition: BorderedBandMatrix.cc:97
gbl::VSymMatrix
Simple symmetric Matrix based on std::vector<double>
Definition: VMatrix.h:86
gbl::VVector::getVec
VVector getVec(unsigned int len, unsigned int start=0) const
Get part of vector.
Definition: VMatrix.cc:273
gbl::BorderedBandMatrix::numBand
unsigned int numBand
Band width.
Definition: BorderedBandMatrix.h:96
gbl::BorderedBandMatrix::invertBand
VMatrix invertBand()
Invert band part.
Definition: BorderedBandMatrix.cc:292
gbl::BorderedBandMatrix::theMixed
VMatrix theMixed
Mixed part.
Definition: BorderedBandMatrix.h:99