MatrixOps.h File Reference

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	DUMMY   0
#define	EPSILON   1e-15
#define	FIRST   1
#define	SECOND   2

Typedefs
typedef vector< vector< double > >	Matrix
typedef vector< double >	Vector

Functions
double	doPartialPivot (Matrix *, unsigned int, unsigned int, vector< unsigned int > *)
Matrix *	matAlloc (unsigned int)
double	matColNorm (Matrix *)
Matrix *	matEyeAdd (const Matrix *, double)
void	matEyeAdd (Matrix *, double, unsigned int)
void	matInv (Matrix *, vector< unsigned int > *, Matrix *)
Matrix *	matMatAdd (const Matrix *, const Matrix *, double, double)
void	matMatAdd (Matrix *, Matrix *, double, double, unsigned int)
Matrix *	matMatMul (Matrix *, Matrix *)
void	matMatMul (Matrix *, Matrix *, unsigned int)
void	matPermMul (Matrix *, vector< unsigned int > *)
void	matPrint (Matrix *)
Matrix *	matScalShift (const Matrix *, double, double)
void	matScalShift (Matrix *, double, double, unsigned int)
double	matTrace (Matrix *)
Matrix *	matTrans (Matrix *)
Vector *	matVecMul (Matrix *, Vector *)
void	triMatMul (Matrix *, Matrix *)
Vector *	vecAlloc (unsigned int)
Vector *	vecMatMul (const Vector *, Matrix *)
void	vecPrint (Vector *)
Vector *	vecScalShift (const Vector *, double, double)
void	vecScalShift (Vector *, double, double, unsigned int)
Vector *	vecVecScalAdd (const Vector *, const Vector *, double, double)
void	vecVecScalAdd (Vector *, Vector *, double, double, unsigned int)

Macro Definition Documentation

#define DUMMY   0

Definition at line 33 of file MatrixOps.h.

Referenced by MarkovSolver::computeMatrixExponential(), and MarkovSolver::computePadeApproximant().

#define EPSILON   1e-15

Definition at line 28 of file MatrixOps.h.

Referenced by Dsolve::calcJnDiff(), ReadSwc::cleanZeroLength(), coordSystem(), Ksolve::getEstimatedDt(), NeuroMesh::insertSingleDummy(), integ(), CylBase::matchCubeMeshEntries(), CylMesh::matchCubeMeshEntries(), CylMesh::matchCylMeshEntries(), matInv(), and CylBase::nearest().

#define FIRST   1

Definition at line 30 of file MatrixOps.h.

Referenced by MarkovSolver::computeMatrixExponential(), MarkovSolver::computePadeApproximant(), matMatAdd(), and matMatMul().

#define SECOND   2

Definition at line 31 of file MatrixOps.h.

Referenced by MarkovSolver::computePadeApproximant(), matMatAdd(), and matMatMul().

Typedef Documentation

typedef vector< vector< double > > Matrix

Definition at line 22 of file MatrixOps.h.

typedef vector< double > Vector

Definition at line 23 of file MatrixOps.h.

Function Documentation

double doPartialPivot	(	Matrix *	,
		unsigned	int,
		unsigned	int,
		vector< unsigned int > *
	)

Definition at line 314 of file MatrixOps.cpp.

References doubleEq().

Referenced by matInv().

{
    unsigned int n = A->size(), pivotRow = row;
    double pivot = (*A)[row][col];

    for( unsigned i = row; i < n; ++i )
    {
        if( fabs( (*A)[i][col] ) > pivot )
        {
            pivotRow = i;
            pivot = (*A)[i][col];
        }
    }

    //If pivot is non-zero, do the row swap.
    if ( !doubleEq(pivot,0.0) && pivotRow != row )
    {
        Matrix::iterator pivotRowItr = (A->begin() + pivotRow);
        Matrix::iterator currRowItr = (A->begin() + row);
        swap( *pivotRowItr, *currRowItr );

        //The row numbers interchanged are stored as a 2-digit number and pushed
        //into this vector. This information is later used in creating the
        //permutation matrices.
        swaps->push_back( 10 * pivotRow + row );
        return pivot;
    }
    else if ( !doubleEq(pivot, 0.0) && pivotRow == row )
        return (*A)[row][col];          //Choice of pivot is unchanged.
    else
        return 0;                                       //Matrix is singular!
}

Here is the call graph for this function:

Here is the caller graph for this function:

Matrix* matAlloc

(

unsigned

int

)

Definition at line 480 of file MatrixOps.cpp.

Referenced by MarkovSolver::computePadeApproximant(), MarkovSolverBase::init(), matEyeAdd(), matInv(), matMatAdd(), matMatMul(), matScalShift(), and matTrans().

{
    Matrix* A = new Matrix;

    A->resize( n );
    for ( unsigned int i = 0; i < n; ++i )
        (*A)[i].resize( n, 0.0 );

    return A;
}

Here is the caller graph for this function:

double matColNorm

(

Matrix *

)

Definition at line 282 of file MatrixOps.cpp.

Referenced by MarkovSolver::computeMatrixExponential().

{
    double norm = 0, colSum = 0;
    unsigned int n = A->size();

    for( unsigned int i = 0; i < n; ++i )
    {
        colSum = 0;
        for( unsigned int j = 0; j < n; ++j )
            colSum += fabs( (*A)[j][i] );

        if ( colSum > norm )
            norm = colSum;
    }

    return norm;
}

Here is the caller graph for this function:

Matrix* matEyeAdd	(	const Matrix *	,
		double
	)

Definition at line 147 of file MatrixOps.cpp.

References matAlloc().

Referenced by MarkovSolver::computeMatrixExponential(), and MarkovSolver::computePadeApproximant().

{
    unsigned int n = A->size();
    Matrix* B = matAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
        {
            if ( i == j )
                (*B)[i][j] = (*A)[i][j] + k;
            else
                (*B)[i][j] = (*A)[i][j];
        }
    }

    return B;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void matEyeAdd	(	Matrix *	,
		double	,
		unsigned	int
	)

Definition at line 166 of file MatrixOps.cpp.

{
    unsigned int n = A->size();
    dummy = 0;

    for( unsigned int i = 0; i < n; ++i )
        (*A)[i][i] += k;
}

void matInv	(	Matrix *	,
		vector< unsigned int > *	,
		Matrix *
	)

Definition at line 348 of file MatrixOps.cpp.

References doPartialPivot(), doubleEq(), EPSILON, matAlloc(), matPermMul(), and triMatMul().

Referenced by MarkovSolver::computePadeApproximant().

{
    Matrix *L, *invL;
    unsigned int n = A->size(), i, j, diagPos;
    double pivot, rowMultiplier1, rowMultiplier2;

    L = matAlloc( n );
    invL = matAlloc( n );
    *invA = *A;

    //The upper triangular portion is stored and inverted in invA

    //Creating a copy of the input matrix, as well as initializing the
    //lower triangular matrix L.
    for (i = 0; i < n; ++i)
        (*L)[i][i] = 1;

    //Calculation of LU decomposition.
    for ( unsigned int k = 0; k < n; ++k )
        pivot = doPartialPivot( invA, k, k, swaps );

    diagPos = 0;
    pivot = 0;
    i = 1;
    j = 0;
    while( diagPos < n - 1 )
    {
        rowMultiplier1 = (*invA)[diagPos][j];       //Pivot element.
        rowMultiplier2 = (*invA)[i][j];

        (*invA)[i][j] = 0;
        for( unsigned int k = j + 1; k < n; ++k )
            (*invA)[i][k] = ( (*invA)[i][k] * rowMultiplier1 -
                                       (*invA)[diagPos][k] *rowMultiplier2 ) / rowMultiplier1;


        (*L)[i][j] = rowMultiplier2 / rowMultiplier1;

        if ( i != n - 1 )
            ++i;
        else
        {
            ++j;
            ++diagPos;
            //In case of zero pivot, divide by a very small number whose value we
            //know.
            if ( doubleEq( (*invA)[diagPos][diagPos], 0.0 ) )
            {
                cerr << "Warning : Singularity detected. Proceeding with computation"
                                "anyway.\n";
                (*invA)[diagPos][diagPos] = EPSILON;
            }

            i = diagPos + 1;
        }
    }
    //End of computation of L and U (which is stored in invA).

    //Obtaining the inverse of invA and L, which is obtained by solving the
    //simple systems Ux = I and Lx= I.

    double sum = 0;

    //We serially solve for the equations Ux = e_n, Ux=e_{n-1} ..., Ux = e1.
    //where, e_k is the k'th elementary basis vector.
    for( int k = n - 1; k >= 0; --k )
    {
        for ( int l = k; l >= 0; --l )
        {
            sum = 0;

            if ( l != k )
            {
                for ( int m = k; m > l; --m )
                    sum += (*invA)[l][m] * (*invA)[m][k];
            }

            if ( l == k )
                (*invA)[l][k] = (1 - sum)/(*invA)[l][l];
            else
                (*invA)[l][k] = -sum/(*invA)[l][l];
        }
    }
    //Similarly as above, we find the inverse of the lower triangular matrix by
    //forward-substitution.

    *invL = *L;
    double invTemp;

    //Always true when the diagonal elements are 1.
    for( unsigned int k = 0; k < n - 1; ++k )
        (*invL)[k+1][k] = -(*invL)[k+1][k];

    for( unsigned int k = 0; k <= n - 1; ++k )
    {
        for( unsigned int l = k + 2; l <= n - 1; ++l )
        {
            invTemp = 0;
            for ( unsigned int m = k + 1; m <= n - 1; ++m )
                invTemp -= (*invL)[m][k] * (*L)[l][m];

            (*invL)[l][k] = invTemp;
        }
    }

    //End of computation of invL and invU. Note that they have been computed in
    //place, which means the original copies of L and U are now gone.

    //Constructing the inverse of the permutation matrix P.
    //P is calculated only if there was any pivoting carried out.
    //At this stage, L^(-1) has already been calculated.
    //P^(-1) = P^T.
    //Since we have computed PA = LU,
    //the final inverse is given by U^(-1)*L^(-1)*P^(-1).
    //If P was not calculated i.e. there were no exchanges, then the
    //inverse is just U^(-1) * L^(-1).
    triMatMul( invA, invL );
    if ( !swaps->empty() )
        matPermMul( invA, swaps );

    delete invL;
    delete L;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Matrix* matMatAdd	(	const Matrix *	,
		const Matrix *	,
		double	,
		double
	)

Definition at line 111 of file MatrixOps.cpp.

References matAlloc().

Referenced by MarkovSolver::computePadeApproximant().

{
    unsigned int n = A->size();
    Matrix *C = matAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*C)[i][j] = alpha * (*A)[i][j] + beta * (*B)[i][j];
    }

    return C;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void matMatAdd	(	Matrix *	,
		Matrix *	,
		double	,
		double	,
		unsigned	int
	)

Definition at line 125 of file MatrixOps.cpp.

References FIRST, and SECOND.

{
    Matrix *C;
    unsigned int n = A->size();

    if ( resIndex == FIRST ) {
        C = A;
    } else if ( resIndex == SECOND ) {
        C = B;
    } else {
        cerr << "matMatAdd : Invalid index supplied to store result.\n";
        return;
    }

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*C)[i][j] = alpha * (*A)[i][j] + beta * (*B)[i][j];
    }
}

Matrix* matMatMul	(	Matrix *	,
		Matrix *
	)

Definition at line 36 of file MatrixOps.cpp.

References matAlloc().

Referenced by MarkovSolver::computeMatrixExponential(), MarkovSolver::computePadeApproximant(), and matMatMul().

{
    unsigned int n = A->size();
    Matrix *C = matAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
        {
            for( unsigned int k = 0; k < n; ++k )
                (*C)[i][j] += (*A)[i][k] * (*B)[k][j];
        }
    }

    return C;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void matMatMul	(	Matrix *	,
		Matrix *	,
		unsigned	int
	)

Definition at line 53 of file MatrixOps.cpp.

References FIRST, matMatMul(), and SECOND.

{
    Matrix *C;

    C = matMatMul( A, B );

    if ( resIndex == FIRST )
        *A = *C;
    else if ( resIndex == SECOND )
        *B = *C;

    delete C;
}

Here is the call graph for this function:

void matPermMul	(	Matrix *	,
		vector< unsigned int > *
	)

Definition at line 88 of file MatrixOps.cpp.

Referenced by matInv().

{
    unsigned int i,j,n,index;
    double temp;
    n = A->size();

    while (!swaps->empty())
    {
        index = swaps->back();
        swaps->pop_back();
        i = index % 10;
        j = (index / 10 ) % 10;

        //Swapping the columns.
        for( unsigned int l = 0; l < n; ++l )
        {
            temp = (*A)[l][i];
            (*A)[l][i] = (*A)[l][j];
            (*A)[l][j] = temp;
        }
    }
}

Here is the caller graph for this function:

void matPrint

(

Matrix *

)

Definition at line 19 of file MatrixOps.cpp.

{
    for( unsigned int i = 0; i < A->size(); ++i )
    {
        for( unsigned int j = 0; j < A->size(); ++j )
            cout << (*A)[i][j] << " ";
        cout << endl;
    }
}

Matrix* matScalShift	(	const Matrix *	,
		double	,
		double
	)

Definition at line 175 of file MatrixOps.cpp.

References matAlloc().

Referenced by MarkovSolver::computeMatrixExponential(), and MarkovSolver::computePadeApproximant().

{
    unsigned int n = A->size();
    Matrix *C = matAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*C)[i][j] = (*A)[i][j] * mul + add;
    }

    return C;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void matScalShift	(	Matrix *	,
		double	,
		double	,
		unsigned	int
	)

Definition at line 189 of file MatrixOps.cpp.

{
    unsigned int n = A->size();
    dummy = 0;

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*A)[i][j] = (*A)[i][j] * mul + add;
    }
}

double matTrace

(

Matrix *

)

Definition at line 271 of file MatrixOps.cpp.

References moose::trace.

Referenced by MarkovSolver::computeMatrixExponential().

{
    unsigned int n = A->size();
    double trace = 0;

    for ( unsigned int i = 0; i < n; ++i )
        trace += (*A)[i][i];

    return trace;
}

Here is the caller graph for this function:

Matrix* matTrans

(

Matrix *

)

Definition at line 300 of file MatrixOps.cpp.

References matAlloc().

{
    unsigned int n = A->size();
    Matrix* At = matAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*At)[i][j] = (*A)[j][i];
    }

    return At;
}

Here is the call graph for this function:

Vector* matVecMul	(	Matrix *	,
		Vector *
	)

Definition at line 235 of file MatrixOps.cpp.

References vecAlloc().

{
    unsigned int n = A->size();
    Vector* w = vecAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*w)[i] += (*A)[i][j] * (*v)[j];
    }

    return w;
}

Here is the call graph for this function:

void triMatMul	(	Matrix *	,
		Matrix *
	)

Definition at line 67 of file MatrixOps.cpp.

Referenced by matInv().

{
    unsigned int n = A->size();
    unsigned int k;
    double temp;

    for ( unsigned int i = 0; i < n; ++i )
    {
        for ( unsigned int j = 0; j < n; ++j )
        {
            k = j >= i ? j : i;
            temp = (*A)[i][j];

            for ( ; k < n; ++k )
                (*A)[i][j] += (*A)[i][k] * (*B)[k][j];

            (*A)[i][j] -= temp;
        }
    }
}

Here is the caller graph for this function:

Vector* vecAlloc

(

unsigned

int

)

Definition at line 491 of file MatrixOps.cpp.

Referenced by matVecMul(), vecMatMul(), vecScalShift(), and vecVecScalAdd().

{
    Vector* vec = new Vector;

    vec->resize( n, 0.0 );

    return vec;
}

Here is the caller graph for this function:

Vector* vecMatMul	(	const Vector *	,
		Matrix *
	)

Definition at line 201 of file MatrixOps.cpp.

References vecAlloc().

Referenced by MarkovSolverBase::bilinearInterpolate(), and MarkovSolverBase::linearInterpolate().

{
    unsigned int n = A->size();
    Vector* w = vecAlloc( n );

    for( unsigned int i = 0; i < n; ++i )
    {
        for( unsigned int j = 0; j < n; ++j )
            (*w)[i] += (*v)[j] * (*A)[j][i];
    }

    return w;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void vecPrint

(

Vector *

)

Definition at line 29 of file MatrixOps.cpp.

{
    for( unsigned int j = 0; j < v->size(); ++j )
        cout << (*v)[j] << " ";
    cout << endl;
}

Vector* vecScalShift	(	const Vector *	,
		double	,
		double
	)

Definition at line 215 of file MatrixOps.cpp.

References vecAlloc().

{
    unsigned int n = v->size();
    Vector* w = vecAlloc( n );

    for ( unsigned int i = 0; i < n; ++i )
        (*w)[i] = scal * (*v)[i] + shift;

    return w;
}

Here is the call graph for this function:

void vecScalShift	(	Vector *	,
		double	,
		double	,
		unsigned	int
	)

Definition at line 226 of file MatrixOps.cpp.

{
    unsigned int n = v->size();
    dummy = 0;

    for ( unsigned int i = 0; i < n; ++i )
        (*v)[i] += scal * (*v)[i] + shift;
}

Vector* vecVecScalAdd	(	const Vector *	,
		const Vector *	,
		double	,
		double
	)

Definition at line 249 of file MatrixOps.cpp.

References vecAlloc().

Referenced by MarkovSolverBase::bilinearInterpolate(), and MarkovSolverBase::linearInterpolate().

{
    unsigned int n = v1->size();
    Vector *w = vecAlloc( n );

    for ( unsigned int i = 0; i < n; ++i )
        (*w)[i] = (*v1)[i] * alpha + (*v2)[i] * beta;

    return w;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void vecVecScalAdd	(	Vector *	,
		Vector *	,
		double	,
		double	,
		unsigned	int
	)

Definition at line 261 of file MatrixOps.cpp.

{
    unsigned int n = v1->size();
    dummy = 0;

    for ( unsigned int i = 0; i < n; ++i )
        (*v1)[i] = (*v1)[i] * alpha + (*v2)[i] * beta;
}