MatrixOps.cpp File Reference

#include <vector>
#include <math.h>
#include "doubleEq.h"
#include <iostream>
#include "MatrixOps.h"

Include dependency graph for MatrixOps.cpp:

Go to the source code of this file.

Functions
double	doPartialPivot (Matrix *A, unsigned int row, unsigned int col, vector< unsigned int > *swaps)
Matrix *	matAlloc (unsigned int n)
double	matColNorm (Matrix *A)
Matrix *	matEyeAdd (const Matrix *A, double k)
void	matEyeAdd (Matrix *A, double k, unsigned int dummy)
void	matInv (Matrix *A, vector< unsigned int > *swaps, Matrix *invA)
Matrix *	matMatAdd (const Matrix *A, const Matrix *B, double alpha, double beta)
void	matMatAdd (Matrix *A, Matrix *B, double alpha, double beta, unsigned int resIndex)
Matrix *	matMatMul (Matrix *A, Matrix *B)
void	matMatMul (Matrix *A, Matrix *B, unsigned int resIndex)
void	matPermMul (Matrix *A, vector< unsigned int > *swaps)
void	matPrint (Matrix *A)
Matrix *	matScalShift (const Matrix *A, double mul, double add)
void	matScalShift (Matrix *A, double mul, double add, unsigned int dummy)
double	matTrace (Matrix *A)
Matrix *	matTrans (Matrix *A)
Vector *	matVecMul (Matrix *A, Vector *v)
void	triMatMul (Matrix *A, Matrix *B)
Vector *	vecAlloc (unsigned int n)
Vector *	vecMatMul (const Vector *v, Matrix *A)
void	vecPrint (Vector *v)
Vector *	vecScalShift (const Vector *v, double scal, double shift)
void	vecScalShift (Vector *v, double scal, double shift, unsigned int dummy)
Vector *	vecVecScalAdd (const Vector *v1, const Vector *v2, double alpha, double beta)
void	vecVecScalAdd (Vector *v1, Vector *v2, double alpha, double beta, unsigned int dummy)

Function Documentation

double doPartialPivot	(	Matrix *	A,
		unsigned int	row,
		unsigned int	col,
		vector< unsigned int > *	swaps
	)

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

n

)

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 *

A

)

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 *	A,
		double	k
	)

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 *	A,
		double	k,
		unsigned int	dummy
	)

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 *	A,
		vector< unsigned int > *	swaps,
		Matrix *	invA
	)

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 *	A,
		const Matrix *	B,
		double	alpha,
		double	beta
	)

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 *	A,
		Matrix *	B,
		double	alpha,
		double	beta,
		unsigned int	resIndex
	)

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 *	A,
		Matrix *	B
	)

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 *	A,
		Matrix *	B,
		unsigned int	resIndex
	)

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 *	A,
		vector< unsigned int > *	swaps
	)

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 *

A

)

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 *	A,
		double	mul,
		double	add
	)

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 *	A,
		double	mul,
		double	add,
		unsigned int	dummy
	)

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 *

A

)

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 *

A

)

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 *	A,
		Vector *	v
	)

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 *	A,
		Matrix *	B
	)

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

n

)

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 *	v,
		Matrix *	A
	)

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 *

v

)

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 *	v,
		double	scal,
		double	shift
	)

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 *	v,
		double	scal,
		double	shift,
		unsigned int	dummy
	)

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 *	v1,
		const Vector *	v2,
		double	alpha,
		double	beta
	)

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 *	v1,
		Vector *	v2,
		double	alpha,
		double	beta,
		unsigned int	dummy
	)

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