38 numSize(0), numBorder(0), numBand(0), numCol(0) {
78 const std::vector<unsigned int> *anIndex,
79 const std::vector<double> *aVector) {
81 for (
unsigned int i = 0; i < anIndex->size(); ++i) {
82 int iIndex = (*anIndex)[i] - 1;
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
88 }
else if (jIndex < nBorder) {
89 theMixed(jIndex, iIndex - nBorder) += (*aVector)[i] * aWeight
92 unsigned int nBand = iIndex - jIndex;
93 theBand(nBand, jIndex - nBorder) += (*aVector)[i] * aWeight
112 unsigned int *anIndex,
double *aVector) {
114 for (
unsigned int i = 0; i < aSize; ++i) {
115 int iIndex = anIndex[i] - 1;
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
124 unsigned int nBand = iIndex - jIndex;
125 theBand(nBand, jIndex - nBorder) += aVector[i] * aWeight
139 const std::vector<unsigned int> anIndex)
const {
141 MatrixXd aMatrix(anIndex.size(), anIndex.size());
143 for (
unsigned int i = 0; i < anIndex.size(); ++i) {
144 int iIndex = anIndex[i] - 1;
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);
149 }
else if (jIndex < nBorder) {
150 aMatrix(i, j) = -
theMixed(jIndex, iIndex - nBorder);
152 unsigned int nBand = iIndex - jIndex;
153 aMatrix(i, j) =
theBand(nBand, jIndex - nBorder);
155 aMatrix(j, i) = aMatrix(i, j);
168 unsigned int *anIndex)
const {
170 MatrixXd aMatrix(aSize, aSize);
172 for (
unsigned int i = 0; i < aSize; ++i) {
173 int iIndex = anIndex[i] - 1;
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);
178 }
else if (jIndex < nBorder) {
179 aMatrix(i, j) = -
theMixed(jIndex, iIndex - nBorder);
181 unsigned int nBand = iIndex - jIndex;
182 aMatrix(i, j) =
theBand(nBand, jIndex - nBorder);
184 aMatrix(j, i) = aMatrix(i, j);
233 const VVector borderSolution = inverseBorder * auxVec;
237 aSolution.
putVec(borderSolution);
251 std::cout <<
"Border part " << std::endl;
253 std::cout <<
"Mixed part " << std::endl;
255 std::cout <<
"Band part " << std::endl;
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));
268 return (diagMin > 0.) ? diagMax / diagMin : 0.;
286 for (
int i = 0; i < nCol; ++i) {
287 auxVec(i) =
theBand(0, i) * 16.0;
289 for (
int i = 0; i < nCol; ++i) {
299 for (
int j = 1; j < std::min(nRow, nCol - i); ++j) {
301 for (
int k = 0; k < std::min(nRow, nCol - i) - j; ++k) {
320 VVector aSolution(aRightHandSide);
321 for (
int i = 0; i < nCol; ++i)
323 for (
int j = 1; j < std::min(nRow, nCol - i); ++j) {
324 aSolution(j + i) -=
theBand(j, i) * aSolution(i);
327 for (
int i = nCol - 1; i >= 0; i--)
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);
349 VMatrix aSolution(aRightHandSide);
350 for (
unsigned int iBorder = 0; iBorder <
numBorder; iBorder++) {
351 for (
int i = 0; i < nCol; ++i)
353 for (
int j = 1; j < std::min(nRow, nCol - i); ++j) {
354 aSolution(iBorder, j + i) -=
theBand(j, i)
355 * aSolution(iBorder, i);
358 for (
int i = nCol - 1; i >= 0; i--)
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);
364 aSolution(iBorder, i) = rxw;
378 VMatrix inverseBand(nRow, nCol);
380 for (
int i = nCol - 1; i >= 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))
387 inverseBand(i - j, j) = rxw;
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) {
408 for (
int l = 0; l < nBorder; ++l) {
409 sum += anArray(i, l) * aSymArray(l, l) * anArray(j, l);
410 for (
int k = 0; k < l; ++k) {
411 sum += anArray(i, l) * aSymArray(l, k) * anArray(j, k)
412 + anArray(i, k) * aSymArray(l, k) * anArray(j, l);
415 aBand(i - j, j) = sum;
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
virtual ~BorderedBandMatrix()
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.
VMatrix theBand
Band part.
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>
unsigned int getNumCols() const
Get number of columns.
VMatrix transpose() const
Get transposed matrix.
unsigned int getNumRows() const
Get number of rows.
void resize(const unsigned int nRows, const unsigned int nCols)
Resize Matrix.
void print() const
Print matrix.
void setZero()
Set content to 0.
Simple symmetric Matrix based on std::vector<double>
void print() const
Print matrix.
unsigned int invert()
Matrix inversion.
void resize(const unsigned int nRows)
Resize symmetric matrix.
void setZero()
Set content to 0.
Simple Vector based on std::vector<double>
void putVec(const VVector &aVector, unsigned int start=0)
Put part of vector.
VVector getVec(unsigned int len, unsigned int start=0) const
Get part of vector.
Namespace for the general broken lines package.