BorderedBandMatrix.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BorderedBandMatrix.cc

/*
 * BorderedBandMatrix.cpp
 *
 *  Created on: Aug 14, 2011
 *      Author: kleinwrt
 */

#include "BorderedBandMatrix.h"

namespace gbl {

BorderedBandMatrix::BorderedBandMatrix() : numSize(0), numBorder(0), numBand(0), numCol(0) {
}

BorderedBandMatrix::~BorderedBandMatrix() {
}


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


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


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

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


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

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();
}

/*============================================================================
 from Dbandmatrix.F (MillePede-II by V. Blobel, Univ. Hamburg)
 ============================================================================*/

void BorderedBandMatrix::decomposeBand() {

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


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

        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 BorderedBandMatrix::solveBand(const VMatrix &aRightHandSide) const {

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


VMatrix BorderedBandMatrix::invertBand() {

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


VMatrix BorderedBandMatrix::bandOfAVAT(const VMatrix &anArray,
                const VSymMatrix &aSymArray) const {
        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;
}

}