Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BorderedBandMatrix.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BorderedBandMatrix.cc
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 
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 }