GeneralBrokenLines V03-01-03
using EIGEN
BorderedBandMatrix.cpp
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"
31using namespace Eigen;
32
34namespace gbl {
35
38 numSize(0), numBorder(0), numBand(0), numCol(0) {
39}
40
42}
43
45
50void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
51 unsigned int nBand) {
52 numSize = nSize;
53 numBorder = nBorder;
54 numCol = nSize - nBorder;
55 numBand = 0;
58 theBand.resize((nBand + 1), numCol);
59}
60
66}
67
69
78 const std::vector<unsigned int> *anIndex,
79 const std::vector<double> *aVector) {
80 int nBorder = numBorder;
81 for (unsigned int i = 0; i < anIndex->size(); ++i) {
82 int iIndex = (*anIndex)[i] - 1; // anIndex has to be sorted
83 for (unsigned int j = 0; j <= i; ++j) {
84 int jIndex = (*anIndex)[j] - 1;
85 if (iIndex < nBorder) {
86 theBorder(iIndex, jIndex) += (*aVector)[i] * aWeight
87 * (*aVector)[j];
88 } else if (jIndex < nBorder) {
89 theMixed(jIndex, iIndex - nBorder) += (*aVector)[i] * aWeight
90 * (*aVector)[j];
91 } else {
92 unsigned int nBand = iIndex - jIndex;
93 theBand(nBand, jIndex - nBorder) += (*aVector)[i] * aWeight
94 * (*aVector)[j];
95 numBand = std::max(numBand, nBand); // update band width
96 }
97 }
98 }
99}
100
102
111void BorderedBandMatrix::addBlockMatrix(double aWeight, unsigned int aSize,
112 unsigned int *anIndex, double *aVector) {
113 int nBorder = numBorder;
114 for (unsigned int i = 0; i < aSize; ++i) {
115 int iIndex = anIndex[i] - 1; // anIndex has to be sorted
116 for (unsigned int j = 0; j <= i; ++j) {
117 int jIndex = anIndex[j] - 1;
118 if (iIndex < nBorder) {
119 theBorder(iIndex, jIndex) += aVector[i] * aWeight * aVector[j];
120 } else if (jIndex < nBorder) {
121 theMixed(jIndex, iIndex - nBorder) += aVector[i] * aWeight
122 * aVector[j];
123 } else {
124 unsigned int nBand = iIndex - jIndex;
125 theBand(nBand, jIndex - nBorder) += aVector[i] * aWeight
126 * aVector[j];
127 numBand = std::max(numBand, nBand); // update band width
128 }
129 }
130 }
131}
132
134
139 const std::vector<unsigned int> anIndex) const {
140
141 MatrixXd aMatrix(anIndex.size(), anIndex.size());
142 int nBorder = numBorder;
143 for (unsigned int i = 0; i < anIndex.size(); ++i) {
144 int iIndex = anIndex[i] - 1; // anIndex has to be sorted
145 for (unsigned int j = 0; j <= i; ++j) {
146 int jIndex = anIndex[j] - 1;
147 if (iIndex < nBorder) {
148 aMatrix(i, j) = theBorder(iIndex, jIndex); // border part of inverse
149 } else if (jIndex < nBorder) {
150 aMatrix(i, j) = -theMixed(jIndex, iIndex - nBorder); // mixed part of inverse
151 } else {
152 unsigned int nBand = iIndex - jIndex;
153 aMatrix(i, j) = theBand(nBand, jIndex - nBorder); // band part of inverse
154 }
155 aMatrix(j, i) = aMatrix(i, j);
156 }
157 }
158 return aMatrix;
159}
160
162
167MatrixXd BorderedBandMatrix::getBlockMatrix(unsigned int aSize,
168 unsigned int *anIndex) const {
169
170 MatrixXd aMatrix(aSize, aSize);
171 int nBorder = numBorder;
172 for (unsigned int i = 0; i < aSize; ++i) {
173 int iIndex = anIndex[i] - 1; // anIndex has to be sorted
174 for (unsigned int j = 0; j <= i; ++j) {
175 int jIndex = anIndex[j] - 1;
176 if (iIndex < nBorder) {
177 aMatrix(i, j) = theBorder(iIndex, jIndex); // border part of inverse
178 } else if (jIndex < nBorder) {
179 aMatrix(i, j) = -theMixed(jIndex, iIndex - nBorder); // mixed part of inverse
180 } else {
181 unsigned int nBand = iIndex - jIndex;
182 aMatrix(i, j) = theBand(nBand, jIndex - nBorder); // band part of inverse
183 }
184 aMatrix(j, i) = aMatrix(i, j);
185 }
186 }
187 return aMatrix;
188}
189
191
218 const VVector &aRightHandSide, VVector &aSolution) {
219
220 // decompose band
222 // invert band
223 VMatrix inverseBand = invertBand();
224 if (numBorder > 0) { // need to use block matrix decomposition to solve
225 // solve for mixed part
226 const VMatrix auxMat = solveBand(theMixed); // = Xt
227 const VMatrix auxMatT = auxMat.transpose(); // = X
228 // solve for border part
229 const VVector auxVec = aRightHandSide.getVec(numBorder)
230 - auxMat * aRightHandSide.getVec(numCol, numBorder);// = b1 - Xt*b2
231 VSymMatrix inverseBorder = theBorder - theMixed * auxMatT;
232 inverseBorder.invert(); // = E
233 const VVector borderSolution = inverseBorder * auxVec; // = x1
234 // solve for band part
235 const VVector bandSolution = solveBand(
236 aRightHandSide.getVec(numCol, numBorder)); // = x
237 aSolution.putVec(borderSolution);
238 aSolution.putVec(bandSolution - auxMatT * borderSolution, numBorder);// = x2
239 // parts of inverse
240 theBorder = inverseBorder; // E
241 theMixed = inverseBorder * auxMat; // E*Xt (-mixed part of inverse) !!!
242 theBand = inverseBand + bandOfAVAT(auxMatT, inverseBorder); // band(D^-1 + X*E*Xt)
243 } else {
244 aSolution.putVec(solveBand(aRightHandSide));
245 theBand = inverseBand;
246 }
247}
248
251 std::cout << "Border part " << std::endl;
253 std::cout << "Mixed part " << std::endl;
254 theMixed.print();
255 std::cout << "Band part " << std::endl;
256 theBand.print();
257}
258
261 // get min. and max. value from diagonal matrix D
262 double diagMin = theBand(0, 0);
263 double diagMax = theBand(0, 0);
264 for (int i = 1; i < numCol; ++i) {
265 diagMin = std::min(diagMin, theBand(0, i));
266 diagMax = std::max(diagMax, theBand(0, i));
267 }
268 return (diagMin > 0.) ? diagMax / diagMin : 0.;
269}
270
271/*============================================================================
272 from Dbandmatrix.F (MillePede-II by V. Blobel, Univ. Hamburg)
273 ============================================================================*/
275
282
283 int nRow = numBand + 1;
284 int nCol = numCol;
285 VVector auxVec(nCol);
286 for (int i = 0; i < nCol; ++i) {
287 auxVec(i) = theBand(0, i) * 16.0; // save diagonal elements
288 }
289 for (int i = 0; i < nCol; ++i) {
290 if ((theBand(0, i) + auxVec(i)) != theBand(0, i)) {
291 theBand(0, i) = 1.0 / theBand(0, i);
292 if (theBand(0, i) < 0.) {
293 throw 3; // not positive definite
294 }
295 } else {
296 theBand(0, i) = 0.0;
297 throw 2; // singular
298 }
299 for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
300 double rxw = theBand(j, i) * theBand(0, i);
301 for (int k = 0; k < std::min(nRow, nCol - i) - j; ++k) {
302 theBand(k, i + j) -= theBand(k + j, i) * rxw;
303 }
304 theBand(j, i) = rxw;
305 }
306 }
307}
308
310
316VVector BorderedBandMatrix::solveBand(const VVector &aRightHandSide) const {
317
318 int nRow = theBand.getNumRows();
319 int nCol = theBand.getNumCols();
320 VVector aSolution(aRightHandSide);
321 for (int i = 0; i < nCol; ++i) // forward substitution
322 {
323 for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
324 aSolution(j + i) -= theBand(j, i) * aSolution(i);
325 }
326 }
327 for (int i = nCol - 1; i >= 0; i--) // backward substitution
328 {
329 double rxw = theBand(0, i) * aSolution(i);
330 for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
331 rxw -= theBand(j, i) * aSolution(j + i);
332 }
333 aSolution(i) = rxw;
334 }
335 return aSolution;
336}
337
339
345VMatrix BorderedBandMatrix::solveBand(const VMatrix &aRightHandSide) const {
346
347 int nRow = theBand.getNumRows();
348 int nCol = theBand.getNumCols();
349 VMatrix aSolution(aRightHandSide);
350 for (unsigned int iBorder = 0; iBorder < numBorder; iBorder++) {
351 for (int i = 0; i < nCol; ++i) // forward substitution
352 {
353 for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
354 aSolution(iBorder, j + i) -= theBand(j, i)
355 * aSolution(iBorder, i);
356 }
357 }
358 for (int i = nCol - 1; i >= 0; i--) // backward substitution
359 {
360 double rxw = theBand(0, i) * aSolution(iBorder, i);
361 for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
362 rxw -= theBand(j, i) * aSolution(iBorder, j + i);
363 }
364 aSolution(iBorder, i) = rxw;
365 }
366 }
367 return aSolution;
368}
369
371
375
376 int nRow = numBand + 1;
377 int nCol = numCol;
378 VMatrix inverseBand(nRow, nCol);
379
380 for (int i = nCol - 1; i >= 0; i--) {
381 double rxw = theBand(0, i);
382 for (int j = i; j >= std::max(0, i - nRow + 1); j--) {
383 for (int k = j + 1; k < std::min(nCol, j + nRow); ++k) {
384 rxw -= inverseBand(abs(i - k), std::min(i, k))
385 * theBand(k - j, j);
386 }
387 inverseBand(i - j, j) = rxw;
388 rxw = 0.;
389 }
390 }
391 return inverseBand;
392}
393
395
399 const VSymMatrix &aSymArray) const {
400 int nBand = numBand;
401 int nCol = numCol;
402 int nBorder = numBorder;
403 double sum;
404 VMatrix aBand((nBand + 1), nCol);
405 for (int i = 0; i < nCol; ++i) {
406 for (int j = std::max(0, i - nBand); j <= i; ++j) {
407 sum = 0.;
408 for (int l = 0; l < nBorder; ++l) { // diagonal
409 sum += anArray(i, l) * aSymArray(l, l) * anArray(j, l);
410 for (int k = 0; k < l; ++k) { // off diagonal
411 sum += anArray(i, l) * aSymArray(l, k) * anArray(j, k)
412 + anArray(i, k) * aSymArray(l, k) * anArray(j, l);
413 }
414 }
415 aBand(i - j, j) = sum;
416 }
417 }
418 return aBand;
419}
420
421}
BorderedBandMatrix definition.
double getBandCondition() const
Get condition from band (decomposition)
VMatrix theMixed
Mixed part.
void printMatrix() const
Print bordered band matrix.
void decomposeBand()
(root free) Cholesky decomposition of band part: C=LDL^T
Eigen::MatrixXd getBlockMatrix(const std::vector< unsigned int > anIndex) const
Retrieve symmetric block matrix.
VSymMatrix theBorder
Border part.
unsigned int numBand
Band width.
VVector solveBand(const VVector &aRightHandSide) const
Solve for band part.
void setZero()
Set content to 0.
unsigned int numCol
Band matrix size.
BorderedBandMatrix()
Create bordered band matrix.
VMatrix bandOfAVAT(const VMatrix &anArray, const VSymMatrix &aSymArray) const
Calculate band part of: 'anArray * aSymArray * anArray.T'.
VMatrix invertBand()
Invert band part.
unsigned int numBorder
Border size.
void solveAndInvertBorderedBand(const VVector &aRightHandSide, VVector &aSolution)
Solve linear equation system, partially calculate inverse.
unsigned int numSize
Matrix size.
void resize(unsigned int nSize, unsigned int nBorder=1, unsigned int nBand=7)
Resize bordered band matrix.
void addBlockMatrix(double aWeight, const std::vector< unsigned int > *anIndex, const std::vector< double > *aVector)
Add symmetric block matrix.
Simple Matrix based on std::vector<double>
Definition: VMatrix.h:64
unsigned int getNumCols() const
Get number of columns.
Definition: VMatrix.cpp:92
VMatrix transpose() const
Get transposed matrix.
Definition: VMatrix.cpp:70
unsigned int getNumRows() const
Get number of rows.
Definition: VMatrix.cpp:84
void resize(const unsigned int nRows, const unsigned int nCols)
Resize Matrix.
Definition: VMatrix.cpp:55
void print() const
Print matrix.
Definition: VMatrix.cpp:97
void setZero()
Set content to 0.
Definition: VMatrix.cpp:62
Simple symmetric Matrix based on std::vector<double>
Definition: VMatrix.h:88
void print() const
Print matrix.
Definition: VMatrix.cpp:199
unsigned int invert()
Matrix inversion.
Definition: VMatrix.cpp:363
void resize(const unsigned int nRows)
Resize symmetric matrix.
Definition: VMatrix.cpp:180
void setZero()
Set content to 0.
Definition: VMatrix.cpp:186
Simple Vector based on std::vector<double>
Definition: VMatrix.h:43
void putVec(const VVector &aVector, unsigned int start=0)
Put part of vector.
Definition: VMatrix.cpp:299
VVector getVec(unsigned int len, unsigned int start=0) const
Get part of vector.
Definition: VMatrix.cpp:288
Namespace for the general broken lines package.