Commit 7bfab4e9 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

cpp: cleanup and fixes for first Eigen based release

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@128 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 5a07f1e8
# - Try to find Eigen3 lib
#
# This module supports requiring a minimum version, e.g. you can do
# find_package(Eigen3 3.1.2)
# to require version 3.1.2 or newer of Eigen3.
#
# Once done this will define
#
# EIGEN3_FOUND - system has eigen lib with correct version
# EIGEN3_INCLUDE_DIR - the eigen include directory
# EIGEN3_VERSION - eigen version
# Copyright (c) 2006, 2007 Montel Laurent, <montel@kde.org>
# Copyright (c) 2008, 2009 Gael Guennebaud, <g.gael@free.fr>
# Copyright (c) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
# Redistribution and use is allowed according to the terms of the 2-clause BSD license.
if(NOT Eigen3_FIND_VERSION)
if(NOT Eigen3_FIND_VERSION_MAJOR)
set(Eigen3_FIND_VERSION_MAJOR 2)
endif(NOT Eigen3_FIND_VERSION_MAJOR)
if(NOT Eigen3_FIND_VERSION_MINOR)
set(Eigen3_FIND_VERSION_MINOR 91)
endif(NOT Eigen3_FIND_VERSION_MINOR)
if(NOT Eigen3_FIND_VERSION_PATCH)
set(Eigen3_FIND_VERSION_PATCH 0)
endif(NOT Eigen3_FIND_VERSION_PATCH)
set(Eigen3_FIND_VERSION "${Eigen3_FIND_VERSION_MAJOR}.${Eigen3_FIND_VERSION_MINOR}.${Eigen3_FIND_VERSION_PATCH}")
endif(NOT Eigen3_FIND_VERSION)
macro(_eigen3_check_version)
file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header)
string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}")
set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}")
string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen3_major_version_match "${_eigen3_version_header}")
set(EIGEN3_MAJOR_VERSION "${CMAKE_MATCH_1}")
string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen3_minor_version_match "${_eigen3_version_header}")
set(EIGEN3_MINOR_VERSION "${CMAKE_MATCH_1}")
set(EIGEN3_VERSION ${EIGEN3_WORLD_VERSION}.${EIGEN3_MAJOR_VERSION}.${EIGEN3_MINOR_VERSION})
if(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION})
set(EIGEN3_VERSION_OK FALSE)
else(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION})
set(EIGEN3_VERSION_OK TRUE)
endif(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION})
if(NOT EIGEN3_VERSION_OK)
message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in ${EIGEN3_INCLUDE_DIR}, "
"but at least version ${Eigen3_FIND_VERSION} is required")
endif(NOT EIGEN3_VERSION_OK)
endmacro(_eigen3_check_version)
if (EIGEN3_INCLUDE_DIR)
# in cache already
_eigen3_check_version()
set(EIGEN3_FOUND ${EIGEN3_VERSION_OK})
else (EIGEN3_INCLUDE_DIR)
find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library
PATHS
${CMAKE_INSTALL_PREFIX}/include
${KDE4_INCLUDE_DIR}
${EIGEN3_DIR}/include
PATH_SUFFIXES eigen3 eigen
)
if(EIGEN3_INCLUDE_DIR)
_eigen3_check_version()
endif(EIGEN3_INCLUDE_DIR)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(Eigen3 DEFAULT_MSG EIGEN3_INCLUDE_DIR EIGEN3_VERSION_OK)
mark_as_advanced(EIGEN3_INCLUDE_DIR)
endif(EIGEN3_INCLUDE_DIR)
......@@ -86,8 +86,10 @@ public:
const std::vector<double>* aVector);
void addBlockMatrix(double aWeight, unsigned int nSimple,
unsigned int* anIndex, double* aVector);
Eigen::MatrixXd getBlockMatrix(const std::vector<unsigned int> anIndex) const;
Eigen::MatrixXd getBlockMatrix(unsigned int aSize, unsigned int* anIndex) const;
Eigen::MatrixXd getBlockMatrix(
const std::vector<unsigned int> anIndex) const;
Eigen::MatrixXd getBlockMatrix(unsigned int aSize,
unsigned int* anIndex) const;
void printMatrix() const;
private:
......
......@@ -36,12 +36,13 @@
#include "VMatrix.h"
#include "Eigen/Core"
typedef Eigen::Matrix<double, 5, 5> Matrix5d;
typedef Eigen::Matrix<double, 2, 7> Matrix27d;
//! Namespace for the general broken lines package
namespace gbl {
typedef Eigen::Matrix<double, 5, 5> Matrix5d;
typedef Eigen::Matrix<double, 2, 7> Matrix27d;
enum dataBlockType {
None, InternalMeasurement, InternalKink, ExternalSeed, ExternalMeasurement
};
......@@ -58,8 +59,8 @@ public:
virtual ~GblData();
void addDerivatives(unsigned int iRow,
const std::vector<unsigned int> &labDer, const Matrix5d &matDer,
unsigned int iOff, const Eigen::MatrixXd &derLocal, unsigned int nLocal,
const Eigen::MatrixXd &derTrans);
unsigned int iOff, const Eigen::MatrixXd &derLocal,
unsigned int nLocal, const Eigen::MatrixXd &derTrans);
void addDerivatives(unsigned int iRow,
const std::vector<unsigned int> &labDer, const Matrix27d &matDer,
unsigned int nLocal, const Eigen::MatrixXd &derTrans);
......
......@@ -87,14 +87,16 @@ public:
const TMatrixD &aDerivatives);
#endif
// input via Eigen
void addMeasurement(const Eigen::MatrixXd &aProjection, const Eigen::VectorXd &aResiduals,
void addMeasurement(const Eigen::MatrixXd &aProjection,
const Eigen::VectorXd &aResiduals,
const Eigen::MatrixXd &aPrecision, double minPrecision = 0.);
void addMeasurement(const Eigen::VectorXd &aResiduals, const Eigen::MatrixXd &aPrecision,
double minPrecision = 0.);
void addScatterer(const Eigen::Vector2d &aResiduals, const Eigen::MatrixXd &aPrecision);
void addMeasurement(const Eigen::VectorXd &aResiduals,
const Eigen::MatrixXd &aPrecision, double minPrecision = 0.);
void addScatterer(const Eigen::Vector2d &aResiduals,
const Eigen::MatrixXd &aPrecision);
void addLocals(const Eigen::MatrixXd &aDerivatives);
void addGlobals(const std::vector<int> &aLabels,
const Eigen::MatrixXd &aDerivatives);
const Eigen::MatrixXd &aDerivatives);
//
unsigned int hasMeasurement() const;
double getMeasPrecMin() const;
......@@ -102,8 +104,8 @@ public:
Vector5d &aPrecision) const;
void getMeasTransformation(Eigen::MatrixXd &aTransformation) const;
bool hasScatterer() const;
void getScatterer(Eigen::Matrix2d &aTransformation, Eigen::Vector2d &aResiduals,
Eigen::Vector2d &aPrecision) const;
void getScatterer(Eigen::Matrix2d &aTransformation,
Eigen::Vector2d &aResiduals, Eigen::Vector2d &aPrecision) const;
void getScatTransformation(Eigen::Matrix2d &aTransformation) const;
unsigned int getNumLocals() const;
const Eigen::MatrixXd& getLocalDerivatives() const;
......@@ -115,8 +117,8 @@ public:
unsigned int getLabel() const;
int getOffset() const;
const Matrix5d& getP2pJacobian() const;
void getDerivatives(int aDirection, Eigen::Matrix2d &matW, Eigen::Matrix2d &matWJ,
Eigen::Vector2d &vecWd) const;
void getDerivatives(int aDirection, Eigen::Matrix2d &matW,
Eigen::Matrix2d &matWJ, Eigen::Vector2d &vecWd) const;
void printPoint(unsigned int level = 0) const;
private:
......
......@@ -97,7 +97,8 @@ public:
TVectorD &aDownWeights);
#endif
unsigned int getLabels(std::vector<unsigned int> &aLabelList) const;
unsigned int getLabels(std::vector<std::vector<unsigned int> > &aLabelList) const;
unsigned int getLabels(
std::vector<std::vector<unsigned int> > &aLabelList) const;
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight,
std::string optionList = "", unsigned int aLabel = 0);
void milleOut(MilleBinary &aMille);
......
......@@ -71,9 +71,9 @@ public:
bool doublePrec = false, unsigned int aSize = 2000);
virtual ~MilleBinary();
void addData(double aMeas, double aErr, unsigned int numLocal,
unsigned int* indLocal, double* derLocal,
const std::vector<int> &labGlobal,
const std::vector<double> &derGlobal);
unsigned int* indLocal, double* derLocal,
const std::vector<int> &labGlobal,
const std::vector<double> &derGlobal);
void writeRecord();
private:
......
......@@ -34,369 +34,369 @@ using namespace Eigen;
namespace gbl {
/// Create bordered band matrix.
BorderedBandMatrix::BorderedBandMatrix() :
numSize(0), numBorder(0), numBand(0), numCol(0) {
}
BorderedBandMatrix::BorderedBandMatrix() :
numSize(0), numBorder(0), numBand(0), numCol(0) {
}
BorderedBandMatrix::~BorderedBandMatrix() {
}
BorderedBandMatrix::~BorderedBandMatrix() {
}
/// Resize bordered band matrix.
/**
* \param nSize [in] Size of matrix
* \param nBorder [in] Size of border (=1 for q/p + additional local parameters)
* \param nBand [in] Band width (usually = 5, for simplified jacobians = 4)
*/
void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
unsigned int nBand) {
numSize = nSize;
numBorder = nBorder;
numCol = nSize - nBorder;
numBand = 0;
theBorder.resize(numBorder);
theMixed.resize(numBorder, numCol);
theBand.resize((nBand + 1), numCol);
}
/**
* \param nSize [in] Size of matrix
* \param nBorder [in] Size of border (=1 for q/p + additional local parameters)
* \param nBand [in] Band width (usually = 5, for simplified jacobians = 4)
*/
void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
unsigned int nBand) {
numSize = nSize;
numBorder = nBorder;
numCol = nSize - nBorder;
numBand = 0;
theBorder.resize(numBorder);
theMixed.resize(numBorder, numCol);
theBand.resize((nBand + 1), numCol);
}
/// 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).
* \param aWeight [in] Weight
* \param anIndex [in] List of rows/colums to be used
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight,
const std::vector<unsigned int>* anIndex,
const std::vector<double>* aVector) {
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
}
/**
* Add (extended) block matrix defined by 'aVector * aWeight * aVector.T'
* to bordered band matrix:
* BBmatrix(anIndex(i),anIndex(j)) += aVector(i) * aWeight * aVector(j).
* \param aWeight [in] Weight
* \param anIndex [in] List of rows/colums to be used
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight,
const std::vector<unsigned int>* anIndex,
const std::vector<double>* aVector) {
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
}
}
}
}
/// 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).
* \param aWeight [in] Weight
* \param aSize [in] Size of block matrix
* \param anIndex [in] List of rows/colums to be used
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight, unsigned int aSize,
unsigned int* anIndex, double* aVector) {
int nBorder = numBorder;
for (unsigned int i = 0; i < aSize; ++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
}
/**
* Add (extended) block matrix defined by 'aVector * aWeight * aVector.T'
* to bordered band matrix:
* BBmatrix(anIndex(i),anIndex(j)) += aVector(i) * aWeight * aVector(j).
* \param aWeight [in] Weight
* \param aSize [in] Size of block matrix
* \param anIndex [in] List of rows/colums to be used
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight, unsigned int aSize,
unsigned int* anIndex, double* aVector) {
int nBorder = numBorder;
for (unsigned int i = 0; i < aSize; ++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
}
}
}
}
/// Retrieve symmetric block matrix.
/**
* Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).
* \param anIndex [in] List of rows/colums to be used
*/
MatrixXd BorderedBandMatrix::getBlockMatrix(
const std::vector<unsigned int> anIndex) const {
/**
* Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).
* \param anIndex [in] List of rows/colums to be used
*/
MatrixXd BorderedBandMatrix::getBlockMatrix(
const std::vector<unsigned int> anIndex) const {
MatrixXd aMatrix(anIndex.size(), 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);
MatrixXd aMatrix(anIndex.size(), 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;
}
return aMatrix;
}
/// Retrieve symmetric block matrix.
/**
* Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).
* \param aSize [in] Matrix size
* \param anIndex [in] Array of rows/colums to be used
*/
MatrixXd BorderedBandMatrix::getBlockMatrix(unsigned int aSize,
unsigned int* anIndex) const {
/**
* Get (compressed) block from bordered band matrix: aMatrix(i,j) = BBmatrix(anIndex(i),anIndex(j)).
* \param aSize [in] Matrix size
* \param anIndex [in] Array of rows/colums to be used
*/
MatrixXd BorderedBandMatrix::getBlockMatrix(unsigned int aSize,
unsigned int* anIndex) const {
MatrixXd aMatrix(aSize, aSize);
int nBorder = numBorder;
for (unsigned int i = 0; i < aSize; ++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);
MatrixXd aMatrix(aSize, aSize);
int nBorder = numBorder;
for (unsigned int i = 0; i < aSize; ++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;
}
return aMatrix;
}
/// 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
*
*
* \param [in] aRightHandSide Right hand side (vector) 'b' of A*x=b
* \param [out] aSolution Solution (vector) x of A*x=b
*/
void BorderedBandMatrix::solveAndInvertBorderedBand(
const VVector &aRightHandSide, VVector &aSolution) {
/**
* 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
*
*
* \param [in] aRightHandSide Right hand side (vector) 'b' of A*x=b
* \param [out] aSolution Solution (vector) x of A*x=b
*/
void BorderedBandMatrix::solveAndInvertBorderedBand(
const VVector &aRightHandSide, VVector &aSolution) {
// 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;
}
// 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;
}
}
/// Print bordered band matrix.
void BorderedBandMatrix::printMatrix() const {
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 BorderedBandMatrix::printMatrix() const {
std::cout << "Border part " << std::endl;
theBorder.print();
std::cout << "Mixed part " << std::endl;
theMixed.print();
std::cout << "Band part " << std::endl;
theBand.print();