(Symmetric) Bordered Band Matrix. More...

#include <BorderedBandMatrix.h>

Public Member Functions
void	addBlockMatrix (double aWeight, const std::vector< unsigned int > anIndex, const std::vector< double > aVector)
	Add symmetric block matrix. More...

	BorderedBandMatrix ()
	Create bordered band matrix. More...

TMatrixDSym	getBlockMatrix (const std::vector< unsigned int > anIndex) const
	Retrieve symmetric block matrix. More...

void	printMatrix () const
	Print bordered band matrix. More...

void	resize (unsigned int nSize, unsigned int nBorder=1, unsigned int nBand=5)
	Resize bordered band matrix. More...

void	solveAndInvertBorderedBand (const VVector &aRightHandSide, VVector &aSolution)
	Solve linear equation system, partially calculate inverse. More...

virtual	~BorderedBandMatrix ()

Private Member Functions
VMatrix	bandOfAVAT (const VMatrix &anArray, const VSymMatrix &aSymArray) const
	Calculate band part of: 'anArray * aSymArray * anArray.T'. More...

void	decomposeBand ()
	(root free) Cholesky decomposition of band part: C=LDL^T More...

VMatrix	invertBand ()
	Invert band part. More...

VVector	solveBand (const VVector &aRightHandSide) const
	Solve for band part. More...

VMatrix	solveBand (const VMatrix &aRightHandSide) const
	solve band part for mixed part (border rows). More...

Private Attributes
unsigned int	numBand
	Band width. More...

unsigned int	numBorder
	Border size. More...

unsigned int	numCol
	Band matrix size. More...

unsigned int	numSize
	Matrix size. More...

VMatrix	theBand
	Band part. More...

VSymMatrix	theBorder
	Border part. More...

VMatrix	theMixed
	Mixed part. More...

Detailed Description

(Symmetric) Bordered Band Matrix.

Separate storage of border, mixed and band parts (as vector<double>).

*  Example for matrix size=8 with border size and band width of two
*
*     +-                                 -+
*     |  B11 B12 M13 M14 M15 M16 M17 M18  |
*     |  B12 B22 M23 M24 M25 M26 M27 M28  |
*     |  M13 M23 C33 C34 C35  0.  0.  0.  |
*     |  M14 M24 C34 C44 C45 C46  0.  0.  |
*     |  M15 M25 C35 C45 C55 C56 C57  0.  |
*     |  M16 M26  0. C46 C56 C66 C67 C68  |
*     |  M17 M27  0.  0. C57 C67 C77 C78  |
*     |  M18 M28  0.  0.  0. C68 C78 C88  |
*     +-                                 -+
*
*  Is stored as::
*
*     +-         -+     +-                         -+
*     |  B11 B12  |     |  M13 M14 M15 M16 M17 M18  |
*     |  B12 B22  |     |  M23 M24 M25 M26 M27 M28  |
*     +-         -+     +-                         -+
*
*                       +-                         -+
*                       |  C33 C44 C55 C66 C77 C88  |
*                       |  C34 C45 C56 C67 C78  0.  |
*                       |  C35 C46 C57 C68  0.  0.  |
*                       +-                         -+
*

Definition at line 56 of file BorderedBandMatrix.h.

Constructor & Destructor Documentation

gbl::BorderedBandMatrix::BorderedBandMatrix ( )

Create bordered band matrix.

Definition at line 14 of file BorderedBandMatrix.cc.

14 : numSize(0), numBorder(0), numBand(0), numCol(0) {

15 }

gbl::BorderedBandMatrix::numCol

unsigned int numCol

Band matrix size.

Definition: BorderedBandMatrix.h:74

gbl::BorderedBandMatrix::numBorder

unsigned int numBorder

Border size.

Definition: BorderedBandMatrix.h:72

gbl::BorderedBandMatrix::numSize

unsigned int numSize

Matrix size.

Definition: BorderedBandMatrix.h:71

gbl::BorderedBandMatrix::numBand

unsigned int numBand

Band width.

Definition: BorderedBandMatrix.h:73

gbl::BorderedBandMatrix::~BorderedBandMatrix ( )

virtual

Definition at line 17 of file BorderedBandMatrix.cc.

17 {

18 }

Member Function Documentation

void gbl::BorderedBandMatrix::addBlockMatrix	(	double	aWeight,
		const std::vector< unsigned int > *	anIndex,
		const std::vector< double > *	aVector
	)

Add symmetric block matrix.

Add (extended) block matrix defined by 'aVector * aWeight * aVector.T' to bordered band matrix: BBmatrix(anIndex(i),anIndex(j)) += aVector(i) * aWeight * aVector(j).

Parameters

aWeight	[in] Weight
anIndex	[in] List of rows/colums to be used
aVector	[in] Vector

Definition at line 46 of file BorderedBandMatrix.cc.

References i, j, bookConverter::max, numBand, numBorder, theBand, theBorder, and theMixed.

Referenced by gbl::GblTrajectory::buildLinearEquationSystem().

                                           {
     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) {
                 theBorder(iIndex, jIndex) += (*aVector)[i] * aWeight
                         * (*aVector)[j];
             } else if (jIndex < nBorder) {
                 theMixed(jIndex, iIndex - nBorder) += (*aVector)[i] * aWeight
                         * (*aVector)[j];
             } else {
                 unsigned int nBand = iIndex - jIndex;
                 theBand(nBand, jIndex - nBorder) += (*aVector)[i] * aWeight
                         * (*aVector)[j];
                 numBand = std::max(numBand, nBand); // update band width
             }
         }
     }
 }

VMatrix gbl::BorderedBandMatrix::bandOfAVAT	(	const VMatrix &	anArray,
		const VSymMatrix &	aSymArray
	)		const

private

Calculate band part of: 'anArray * aSymArray * anArray.T'.

Returns: Band part of product

Definition at line 294 of file BorderedBandMatrix.cc.

References i, j, relval_2017::k, cmsLHEtoEOSManager::l, bookConverter::max, numBand, numBorder, and numCol.

Referenced by solveAndInvertBorderedBand().

                                            {
     int nBand = numBand;
     int nCol = numCol;
     int nBorder = numBorder;
     double sum;
     VMatrix aBand((nBand + 1), nCol);
     for (int i = 0; i < nCol; ++i) {
         for (int j = std::max(0, i - nBand); j <= i; ++j) {
             sum = 0.;
             for (int l = 0; l < nBorder; ++l) { // diagonal
                 sum += anArray(i, l) * aSymArray(l, l) * anArray(j, l);
                 for (int k = 0; k < l; ++k) { // off diagonal
                     sum += anArray(i, l) * aSymArray(l, k) * anArray(j, k)
                             + anArray(i, k) * aSymArray(l, k) * anArray(j, l);
                 }
             }
             aBand(i - j, j) = sum;
         }
     }
     return aBand;
 }

void gbl::BorderedBandMatrix::decomposeBand ( )

private

(root free) Cholesky decomposition of band part: C=LDL^T

Decompose band matrix into diagonal matrix D and lower triangular band matrix L (diagonal=1). Overwrite band matrix with D and off-diagonal part of L.

Exceptions

2	: matrix is singular.
3	: matrix is not positive definite.

Definition at line 177 of file BorderedBandMatrix.cc.

References i, j, relval_2017::k, min(), numBand, numCol, and theBand.

Referenced by solveAndInvertBorderedBand().

                                        {
 
     int nRow = numBand + 1;
     int nCol = numCol;
     VVector auxVec(nCol);
     for (int i = 0; i < nCol; ++i) {
         auxVec(i) = theBand(0, i) * 16.0; // save diagonal elements
     }
     for (int i = 0; i < nCol; ++i) {
         if ((theBand(0, i) + auxVec(i)) != theBand(0, i)) {
             theBand(0, i) = 1.0 / theBand(0, i);
             if (theBand(0, i) < 0.) {
                 throw 3; // not positive definite
             }
         } else {
             theBand(0, i) = 0.0;
             throw 2; // singular
         }
         for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
             double rxw = theBand(j, i) * theBand(0, i);
             for (int k = 0; k < std::min(nRow, nCol - i) - j; ++k) {
                 theBand(k, i + j) -= theBand(k + j, i) * rxw;
             }
             theBand(j, i) = rxw;
         }
     }
 }

TMatrixDSym gbl::BorderedBandMatrix::getBlockMatrix ( const std::vector< unsigned int > anIndex ) const

Retrieve symmetric block matrix.

Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).

Parameters

anIndex [in] List of rows/colums to be used

Definition at line 75 of file BorderedBandMatrix.cc.

References i, j, numBorder, theBand, theBorder, and theMixed.

Referenced by gbl::GblTrajectory::getResAndErr(), and gbl::GblTrajectory::getResults().

                                                      {
 
     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;
 }

VMatrix gbl::BorderedBandMatrix::invertBand ( )

private

Invert band part.

Returns: Inverted band

Definition at line 270 of file BorderedBandMatrix.cc.

References funct::abs(), i, j, relval_2017::k, bookConverter::max, min(), numBand, numCol, and theBand.

Referenced by solveAndInvertBorderedBand().

                                        {
 
     int nRow = numBand + 1;
     int nCol = numCol;
     VMatrix inverseBand(nRow, nCol);
 
     for (int i = nCol - 1; i >= 0; i--) {
         double rxw = theBand(0, i);
         for (int j = i; j >= std::max(0, i - nRow + 1); j--) {
             for (int k = j + 1; k < std::min(nCol, j + nRow); ++k) {
                 rxw -= inverseBand(abs(i - k), std::min(i, k))
                         * theBand(k - j, j);
             }
             inverseBand(i - j, j) = rxw;
             rxw = 0.;
         }
     }
     return inverseBand;
 }

void gbl::BorderedBandMatrix::printMatrix ( ) const

Print bordered band matrix.

Definition at line 158 of file BorderedBandMatrix.cc.

References gather_cfg::cout, gbl::VMatrix::print(), gbl::VSymMatrix::print(), theBand, theBorder, and theMixed.

Referenced by gbl::GblTrajectory::printTrajectory().

                                            {
     std::cout << "Border part " << std::endl;
     theBorder.print();
     std::cout << "Mixed  part " << std::endl;
     theMixed.print();
     std::cout << "Band   part " << std::endl;
     theBand.print();
 }

void gbl::BorderedBandMatrix::resize	(	unsigned int	nSize,
		unsigned int	nBorder = `1`,
		unsigned int	nBand = `5`
	)

Resize bordered band matrix.

Parameters

nSize	[in] Size of matrix
nBorder	[in] Size of border (=1 for q/p + additional local parameters)
nBand	[in] Band width (usually = 5, for simplified jacobians = 4)

Definition at line 26 of file BorderedBandMatrix.cc.

References numBand, numBorder, numCol, numSize, gbl::VMatrix::resize(), gbl::VSymMatrix::resize(), theBand, theBorder, and theMixed.

Referenced by Vispa.Gui.TextDialog.TextDialog::__init__(), Vispa.Plugins.ConfigEditor.ToolDialog.ToolDialog::__init__(), Vispa.Main.MainWindow.MainWindow::_loadIni(), gbl::GblTrajectory::buildLinearEquationSystem(), and Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                             {
     numSize = nSize;
     numBorder = nBorder;
     numCol = nSize - nBorder;
     numBand = 0;
     theBorder.resize(numBorder);
     theMixed.resize(numBorder, numCol);
     theBand.resize((nBand + 1), numCol);
 }

void gbl::BorderedBandMatrix::solveAndInvertBorderedBand	(	const VVector &	aRightHandSide,
		VVector &	aSolution
	)

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 part for block matrix algebra:

| A  Ct |   | x1 |   | b1 |        , A  is the border part
|       | * |    | = |    |        , Ct is the mixed part
| C  D  |   | x2 |   | b2 |        , D  is the band part

Explicit inversion of D is avoided by using solution X of D*X=C (X=D^-1*C, obtained from Cholesky decomposition and forward/backward substitution)

| x1 |   | E*b1 - E*Xt*b2 |        , E^-1 = A-Ct*D^-1*C = A-Ct*X
|    | = |                |
| x2 |   |  x   - X*x1    |        , x is solution of D*x=b2 (x=D^-1*b2)

Inverse matrix is:

|  E   -E*Xt          |
|                     |            , only band part of (D^-1 + X*E*Xt)
| -X*E  D^-1 + X*E*Xt |              is calculated

Parameters

[in]	aRightHandSide	Right hand side (vector) 'b' of A*x=b
[out]	aSolution	Solution (vector) x of A*x=b

Definition at line 125 of file BorderedBandMatrix.cc.

References bandOfAVAT(), decomposeBand(), gbl::VVector::getVec(), gbl::VSymMatrix::invert(), invertBand(), numBorder, numCol, gbl::VVector::putVec(), solveBand(), theBand, theBorder, theMixed, and gbl::VMatrix::transpose().

Referenced by gbl::GblTrajectory::fit().

                                                            {
 
     // decompose band
     decomposeBand();
     // invert band
     VMatrix inverseBand = invertBand();
     if (numBorder > 0) { // need to use block matrix decomposition to solve
         // solve for mixed part
         const VMatrix auxMat = solveBand(theMixed); // = Xt
         const VMatrix auxMatT = auxMat.transpose(); // = X
         // solve for border part
         const VVector auxVec = aRightHandSide.getVec(numBorder)
                 - auxMat * aRightHandSide.getVec(numCol, numBorder); // = b1 - Xt*b2
         VSymMatrix inverseBorder = theBorder - theMixed * auxMatT;
         inverseBorder.invert(); // = E
         const VVector borderSolution = inverseBorder * auxVec; // = x1
         // solve for band part
         const VVector bandSolution = solveBand(
                 aRightHandSide.getVec(numCol, numBorder)); // = x
         aSolution.putVec(borderSolution);
         aSolution.putVec(bandSolution - auxMatT * borderSolution, numBorder); // = x2
         // parts of inverse
         theBorder = inverseBorder; // E
         theMixed = inverseBorder * auxMat; // E*Xt (-mixed part of inverse) !!!
         theBand = inverseBand + bandOfAVAT(auxMatT, inverseBorder); // band(D^-1 + X*E*Xt)
     } else {
         aSolution.putVec(solveBand(aRightHandSide));
         theBand = inverseBand;
     }
 }

VVector gbl::BorderedBandMatrix::solveBand ( const VVector & aRightHandSide ) const

private

Solve for band part.

Solve C*x=b for band part using decomposition C=LDL^T and forward (L*z=b) and backward substitution (L^T*x=D^-1*z).

Parameters

[in] aRightHandSide Right hand side (vector) 'b' of C*x=b

Returns: Solution (vector) 'x' of C*x=b

Definition at line 212 of file BorderedBandMatrix.cc.

References gbl::VMatrix::getNumCols(), gbl::VMatrix::getNumRows(), i, j, min(), and theBand.

Referenced by solveAndInvertBorderedBand().

                                                                          {
 
     int nRow = theBand.getNumRows();
     int nCol = theBand.getNumCols();
     VVector aSolution(aRightHandSide);
     for (int i = 0; i < nCol; ++i) // forward substitution
             {
         for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
             aSolution(j + i) -= theBand(j, i) * aSolution(i);
         }
     }
     for (int i = nCol - 1; i >= 0; i--) // backward substitution
             {
         double rxw = theBand(0, i) * aSolution(i);
         for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
             rxw -= theBand(j, i) * aSolution(j + i);
         }
         aSolution(i) = rxw;
     }
     return aSolution;
 }

VMatrix gbl::BorderedBandMatrix::solveBand ( const VMatrix & aRightHandSide ) const

private

solve band part for mixed part (border rows).

Solve C*X=B for mixed part using decomposition C=LDL^T and forward and backward substitution.

Parameters

[in] aRightHandSide Right hand side (matrix) 'B' of C*X=B

Returns: Solution (matrix) 'X' of C*X=B

Definition at line 241 of file BorderedBandMatrix.cc.

References gbl::VMatrix::getNumCols(), gbl::VMatrix::getNumRows(), i, j, min(), numBorder, and theBand.

                                                                          {
 
     int nRow = theBand.getNumRows();
     int nCol = theBand.getNumCols();
     VMatrix aSolution(aRightHandSide);
     for (unsigned int iBorder = 0; iBorder < numBorder; iBorder++) {
         for (int i = 0; i < nCol; ++i) // forward substitution
                 {
             for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
                 aSolution(iBorder, j + i) -= theBand(j, i)
                         * aSolution(iBorder, i);
             }
         }
         for (int i = nCol - 1; i >= 0; i--) // backward substitution
                 {
             double rxw = theBand(0, i) * aSolution(iBorder, i);
             for (int j = 1; j < std::min(nRow, nCol - i); ++j) {
                 rxw -= theBand(j, i) * aSolution(iBorder, j + i);
             }
             aSolution(iBorder, i) = rxw;
         }
     }
     return aSolution;
 }