gbl::VSymMatrix Class Reference

Simple symmetric Matrix based on std::vector<double> More...

#include <VMatrix.h>

Public Member Functions
	VSymMatrix (const unsigned int nRows=0)
virtual	~VSymMatrix ()
void	resize (const unsigned int nRows)
	Resize symmetric matrix.
unsigned int	invert ()
	Matrix inversion.
double &	operator() (unsigned int i, unsigned int j)
	access element (i,j) assuming i>=j
double	operator() (unsigned int i, unsigned int j) const
	access element (i,j) assuming i>=j
unsigned int	getNumRows () const
	Get number of rows (= number of colums).
void	print () const
	Print matrix.
VSymMatrix	operator- (const VMatrix &aMatrix) const
	Subtraction SymMatrix-(sym)Matrix.
VVector	operator* (const VVector &aVector) const
	Multiplication SymMatrix*Vector.
VMatrix	operator* (const VMatrix &aMatrix) const
	Multiplication SymMatrix*Matrix.

Private Attributes
unsigned int	numRows
	Number of rows.
std::vector< double >	theVec
	Data (symmetric storage)

Detailed Description

Simple symmetric Matrix based on std::vector<double>

Definition at line 65 of file VMatrix.h.

Constructor & Destructor Documentation

VSymMatrix()

gbl::VSymMatrix::VSymMatrix

(

const unsigned int

nRows = 0

)

Definition at line 142 of file VMatrix.cc.

                                               :
        numRows(nRows), theVec((nRows * nRows + nRows) / 2) {
}

~VSymMatrix()

gbl::VSymMatrix::~VSymMatrix ( )

virtual

Definition at line 146 of file VMatrix.cc.

                        {
}

Member Function Documentation

getNumRows()

unsigned int gbl::VSymMatrix::getNumRows

(

)

const

Get number of rows (= number of colums).

Returns

Number of rows.

Definition at line 162 of file VMatrix.cc.

                                          {
    return numRows;
}

invert()

unsigned int gbl::VSymMatrix::invert

(

)

Matrix inversion.

Invert symmetric N-by-N matrix V in symmetric storage mode V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, . . . replaced by inverse matrix

Method of solution is by elimination selecting the pivot on the diagonal each stage. The rank of the matrix is returned in NRANK. For NRANK ne N, all remaining rows and cols of the resulting matrix V are set to zero.

Exceptions

1	: matrix is singular.

Returns

Rank of matrix.

Definition at line 326 of file VMatrix.cc.

                                {
 
    const double eps = 1.0E-10;
    unsigned int aSize = numRows;
    std::vector<int> next(aSize);
    std::vector<double> diag(aSize);
    int nSize = aSize;
 
    int first = 1;
    for (int i = 1; i <= nSize; ++i) {
        next[i - 1] = i + 1; // set "next" pointer
        diag[i - 1] = fabs(theVec[(i * i + i) / 2 - 1]); // save abs of diagonal elements
    }
    next[aSize - 1] = -1; // end flag
 
    unsigned int nrank = 0;
    for (int i = 1; i <= nSize; ++i) { // start of loop
        int k = 0;
        double vkk = 0.0;
 
        int jCandidate = first;
        int previous = 0;
        int last = previous;
        // look for pivot
        while (jCandidate > 0) {
            int jj = (jCandidate * jCandidate + jCandidate) / 2 - 1;
            if (fabs(theVec[jj]) > std::max(fabs(vkk), eps * diag[jCandidate - 1])) {
                vkk = theVec[jj];
                k = jCandidate;
                last = previous;
            }
            previous = jCandidate;
            jCandidate = next[jCandidate - 1];
        }
        // pivot found
        if (k > 0) {
            int kk = (k * k + k) / 2 - 1;
            if (last <= 0) {
                first = next[k - 1];
            } else {
                next[last - 1] = next[k - 1];
            }
            next[k - 1] = 0; // index is used, reset
            nrank++; // increase rank and ...
 
            vkk = 1.0 / vkk;
            theVec[kk] = -vkk;
            int jk = kk - k;
            int jl = -1;
            for (int j = 1; j <= nSize; ++j) { // elimination
                if (j == k) {
                    jk = kk;
                    jl += j;
                } else {
                    if (j < k) {
                        ++jk;
                    } else {
                        jk += j - 1;
                    }
 
                    double vjk = theVec[jk];
                    theVec[jk] = vkk * vjk;
                    int lk = kk - k;
                    if (j >= k) {
                        for (int l = 1; l <= k - 1; ++l) {
                            ++jl;
                            ++lk;
                            theVec[jl] -= theVec[lk] * vjk;
                        }
                        ++jl;
                        lk = kk;
                        for (int l = k + 1; l <= j; ++l) {
                            ++jl;
                            lk += l - 1;
                            theVec[jl] -= theVec[lk] * vjk;
                        }
                    } else {
                        for (int l = 1; l <= j; ++l) {
                            ++jl;
                            ++lk;
                            theVec[jl] -= theVec[lk] * vjk;
                        }
                    }
                }
            }
        } else {
            for (int n = 1; n <= nSize; ++n) {
                if (next[n - 1] >= 0) {
                    int nn = (n * n - n) / 2 - 1;
                    for (int j = 1; j <= nSize; ++j) {
                        if (next[j - 1] >= 0) {
                            theVec[nn + j] = 0.0; // clear matrix row/col
                        }
                    }
                }
            }
            throw 1; // singular
        }
    }
    for (int ij = 0; ij < (nSize * nSize + nSize) / 2; ++ij) {
        theVec[ij] = -theVec[ij]; // finally reverse sign of all matrix elements
    }
    return nrank;
}

operator()() [1/2]

double & gbl::VSymMatrix::operator() ( unsigned int  i, unsigned int  j  )

inline

access element (i,j) assuming i>=j

Definition at line 104 of file VMatrix.h.

                                                                          {
    return theVec[(iRow * iRow + iRow) / 2 + iCol]; // assuming iCol <= iRow
}

operator()() [2/2]

double gbl::VSymMatrix::operator() ( unsigned int  i, unsigned int  j  ) const

inline

access element (i,j) assuming i>=j

Definition at line 109 of file VMatrix.h.

                                 {
    return theVec[(iRow * iRow + iRow) / 2 + iCol]; // assuming iCol <= iRow
}

operator*() [1/2]

VMatrix gbl::VSymMatrix::operator*

(

const VMatrix &

aMatrix

)

const

Multiplication SymMatrix*Matrix.

Definition at line 207 of file VMatrix.cc.

                                                          {
    unsigned int nCol = aMatrix.getNumCols();
    VMatrix aResult(numRows, nCol);
    for (unsigned int l = 0; l < nCol; ++l) {
        for (unsigned int i = 0; i < numRows; ++i) {
            aResult(i, l) = theVec[(i * i + i) / 2 + i] * aMatrix(i, l);
            for (unsigned int j = 0; j < i; ++j) {
                aResult(j, l) += theVec[(i * i + i) / 2 + j] * aMatrix(i, l);
                aResult(i, l) += theVec[(i * i + i) / 2 + j] * aMatrix(j, l);
            }
        }
    }
    return aResult;
}

operator*() [2/2]

VVector gbl::VSymMatrix::operator*

(

const VVector &

aVector

)

const

Multiplication SymMatrix*Vector.

Definition at line 194 of file VMatrix.cc.

                                                          {
    VVector aResult(numRows);
    for (unsigned int i = 0; i < numRows; ++i) {
        aResult(i) = theVec[(i * i + i) / 2 + i] * aVector(i);
        for (unsigned int j = 0; j < i; ++j) {
            aResult(j) += theVec[(i * i + i) / 2 + j] * aVector(i);
            aResult(i) += theVec[(i * i + i) / 2 + j] * aVector(j);
        }
    }
    return aResult;
}

operator-()

VSymMatrix gbl::VSymMatrix::operator-

(

const VMatrix &

aMatrix

)

const

Subtraction SymMatrix-(sym)Matrix.

Definition at line 183 of file VMatrix.cc.

                                                             {
    VSymMatrix aResult(numRows);
    for (unsigned int i = 0; i < numRows; ++i) {
        for (unsigned int j = 0; j <= i; ++j) {
            aResult(i, j) = theVec[(i * i + i) / 2 + j] - aMatrix(i, j);
        }
    }
    return aResult;
}

print()

void gbl::VSymMatrix::print

(

)

const

Print matrix.

Definition at line 167 of file VMatrix.cc.

                             {
    std::cout << " VSymMatrix: " << numRows << "*" << numRows << std::endl;
    for (unsigned int i = 0; i < numRows; ++i) {
        for (unsigned int j = 0; j <= i; ++j) {
            if (j % 5 == 0) {
                std::cout << std::endl << std::setw(4) << i << ","
                        << std::setw(4) << j << "-" << std::setw(4)
                        << std::min(j + 4, i) << ":";
            }
            std::cout << std::setw(13) << theVec[(i * i + i) / 2 + j];
        }
    }
    std::cout << std::endl << std::endl;
}

resize()

void gbl::VSymMatrix::resize

(

const unsigned int

nRows

)

Resize symmetric matrix.

Parameters

[in]

nRows

Number of rows.

Definition at line 153 of file VMatrix.cc.

                                                {
    numRows = nRows;
    theVec.resize((nRows * nRows + nRows) / 2);
}

Member Data Documentation

numRows

unsigned int gbl::VSymMatrix::numRows

private

Number of rows.

Definition at line 79 of file VMatrix.h.

theVec

std::vector<double> gbl::VSymMatrix::theVec

private

Data (symmetric storage)

Definition at line 80 of file VMatrix.h.

---

The documentation for this class was generated from the following files:

• genfit/GBL/include/VMatrix.h
• genfit/GBL/src/VMatrix.cc