matrixFunctionInterface.cpp

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////////
/*! \file matrixFunctionInterface.cpp
 *  \brief Miscellaneous matrix functions
 *  \author     $Author: jayhawk_hokie $
 *  \version $Revision: 1.3 $
 *  \date    $Date: 2007/07/24 08:40:44 $
*////////////////////////////////////////////////////////////////////////////
/*
*/
#include "matrixFunctionInterface.h"

using namespace std;

using namespace O_SESSAME;

/*! \brief Takes in a square matrix and returns, by reference, the Real and Imaginary Eigenvalues and the Eigenvectors for that matrix
 *      This function does not normalize the Eigenvectors.
 *  @param _Matrix is a matrix passed to find the eigenvalues and eigenvectors
 *  @return _EigenVectors is a matrix of the eigen vectors
 *  @return _RealEigenValues is a vector of the real eigenvalues of matrix _Matrix
 *  @return _ImaginaryEigenValues is a vector of the imaginary eigenvalues of matrix _Matrix
 */
void eigenValues( Matrix _Matrix, Matrix &_EigenVectors, Vector &_RealEigenValues, Vector &_ImaginaryEigenValues )
{
        // Get matrix size.
        int rows = _Matrix[1].getIndexCount();
        int columns = _Matrix[2].getIndexCount();

        // Check for a square matrix.
        if ( rows != columns )
        {
                cerr << "WARNING: Not a square matrix. Not able to determine EigenValues in the eigenValues function." <<endl;
                return;
        }       

        // Initialize the REAL ** matrices and REAL * vectors for the eigen function.
        void *vmblock;
        REAL **matrix;
        REAL **eigenVectors;
        REAL *realEigenValues;
        REAL *imaginaryEigenValues;
        int *count;

        // Copies _Matrix into matrix for the eigen function.
        vmblock = vminit();
        matrix = (REAL **) vmalloc( vmblock, MATRIX, rows, columns );
        matrix = convertMatrix( _Matrix );
        eigenVectors = (REAL **) vmalloc( vmblock, MATRIX, rows, columns );
        realEigenValues = (REAL *) vmalloc( vmblock, VEKTOR, rows, 0 );
        imaginaryEigenValues = (REAL *) vmalloc( vmblock, VEKTOR, rows, 0 );
        count = (int *) vmalloc( vmblock, VVEKTOR, rows, sizeof( *count ) );
        
        // Takes the matrix of size (rows, rows) and finds the real and imaginary eigenvalues and the eigenvectors.
        eigen( 1, 1, 1, rows, matrix, eigenVectors, realEigenValues, imaginaryEigenValues, count );

        // Copies the Eigenvalues and Eigenvectors into a usable form (of type Matrix and Vector).
        _EigenVectors = convertMatrix(eigenVectors, rows);
        _RealEigenValues = convertVector(realEigenValues, rows);
        _ImaginaryEigenValues = convertVector(imaginaryEigenValues, rows);

        return;
}


/* /brief Calculate the greatest eigenvalue of a real
 *  square matrix and the associated eigenvector
 *  by the power method.
 *  @param _Matrix is a matrix passed to find the greatest eigenvalue and associated eigenvector
 *  @return _EigenValue the largest eigenvalue of _Matrix
 *  @return _EigenVector the eigenvector associated with _EigenValue
 *  @param _MaxIterations is the maximum iterations this function will use in finding _EigenValue                     
 */
bool greatestEigenValue( Matrix _Matrix, double &_EigenValue, Vector &_EigenVector, int _MaxIterations ) 
{
        // Maximum size of the matrix being passed in 
        int NMAX = 26;
        
        // The smallest double 
        double eps = .0000000001;

        // Required precision
        double dta = .00001;
        
        // Used to find any errors while calculating the greatest eigenvalue
        int isError;    
                
        // Get matrix size.
        int rows = _Matrix[1].getIndexCount();
        int columns = _Matrix[2].getIndexCount();

        // Check for a square matrix.
        if ( rows != columns )
        {
                cerr << "WARNING: Not a square matrix. Not able to determine EigenValues in the eigenValues function." <<endl;
                return false;
        }       
        
        // Used in the for loops
        int i, j, iterator;
        
        // phi and s are used to check the precision of the eigenvalue
        double phi, s;

        // Temporary eigenvector
        Vector eigenVectorTemp( NMAX );

        // Initialize _EigenVector to the right size
        _EigenVector = eigenVectorTemp;
 
        for( i = 1; i <= rows; i++ )
                eigenVectorTemp(i) = 1.0 / sqrt( double(i) );
        isError = -1; 
        iterator = 1;

        while ( isError == -1 && iterator <= _MaxIterations )
        {
                _EigenValue = 0.0;
                for ( i = 1; i <= rows; i++ )
                {
                        _EigenVector(i) = 0.0;
                        for ( j = 1; j <= rows; j++ )
                                _EigenVector(i) += _Matrix(i,j) * eigenVectorTemp(j);
                        if ( fabs( _EigenVector(i) ) > fabs( _EigenValue ) ) 
                                _EigenValue = _EigenVector(i);
                }

                if ( fabs( _EigenValue ) < eps ) 
                        isError = 0;
                else
                {
                        for ( i = 1; i <= rows; i++ )  
                                _EigenVector(i) /= _EigenValue;
                        phi = 0.0;
                        for ( i = 1; i <= rows; i++ ) 
                        {
                                s = fabs( _EigenVector(i) - eigenVectorTemp(i) );
                                if ( s > phi )
                                        phi = s;
                        }
                        if ( phi < dta )
                                isError = 1;
                        else
                        {
                                for ( i = 1; i <= rows; i++ ) 
                                        eigenVectorTemp(i) = _EigenVector(i);
                                iterator++;
                        }
                 }
 
        }
        
        // Check to see if any errors had occured
        switch(isError+1) {
                case 0: cerr << "WARNING: No convergence.\n";
                                return false;
                case 1: cerr << "WARNING: Method does not apply.\n"; break;
                                return false;
        }

        // If there were no errors, return true
        return true;
}


/*! \brief Takes in a square matrix and returns the norm of the matrix, which is the square root of the largest eigenvalue of A(transpose)*A 
 *  @param _Matrix is a matrix passed to find the norm
 *  @return norm2 is the norm of the matrix _Matrix
 */
double norm2( Matrix _Matrix )
{
        // Variables to me used in the eigenValue function
        Vector eigenVector;
        double largestEigenValue = 0;
        
        // Find the transpose of _Matrix
        Matrix transposeMatrix = ~_Matrix;
        
        // Variables used in this function to find the norm
        double norm2;
        
        // Find the greatest eigenvalues and its associated eigenvector
        if( greatestEigenValue( ( transposeMatrix * _Matrix ), largestEigenValue, eigenVector, 20 ) )
        {
                // Solve for the norm of _Matrix
                //norm2 = sqrt( largestEigenValue );
        }
        else
        {
                cerr << "Trying alternate method for max eigenvalue." << endl;
                Matrix _EigenVectors(3,3);
                Vector _RealEigenValues(3);
                Vector _ImaginaryEigenValues(3);
                eigenValues( transposeMatrix * _Matrix, _EigenVectors, _RealEigenValues, _ImaginaryEigenValues );

                largestEigenValue = _RealEigenValues.maxAbs( );
                
                if( largestEigenValue < 1e-4 )
                {
                        // Error if the eigenvalue function returns false
                        cerr << "WARNING: Largest eigenvalue cannot be found. Set to zero." << endl;
                }
        }
        
        // Solve for the norm of _Matrix
        norm2 = sqrt( largestEigenValue );
        
        return( norm2 );
        
}

/*! \brief Takes in a Matrix and returns a REAL ** with the same values 
 *  @param _Matrix is a Matrix 
 *  @return finalMatrix the Matrix in REAL **
 */
REAL ** convertMatrix( Matrix _Matrix )
{
        int n = _Matrix[1].getIndexCount();
        void *vmblock;
        REAL **finalMatrix;

        vmblock = vminit();
        finalMatrix = (REAL **) vmalloc( vmblock, MATRIX, n, n );

        for( int i = 1; i <= n ; i++ )
        {
                for( int j = 1; j <= n ; j++ )
                {
                        finalMatrix[i-1][j-1] = _Matrix(i,j);
                }
        }
        return finalMatrix;
}

/*! \brief Takes in a Vector and returns a REAL * with the same values 
 *  @param _Matrix is a Vector 
 *  @return finalMatrix the Matrix in REAL *
 */
REAL * convertVector( Vector _Vector )
{
        int n = _Vector.getIndexCount();
        void *vmblock;
        REAL *finalVector;

        vmblock = vminit();
        finalVector = (REAL *) vmalloc( vmblock, VEKTOR, n, 0 );

        for( int i = 1; i <= n; i++ )
        {
                finalVector[i-1] = _Vector(i);
        }
        return finalVector;
}

/*! \brief Takes in a REAL ** and returns a Matrix with the same values 
 *  @param finalMatrix the Matrix in REAL **
 *  @return _Matrix is a Matrix 
 */
Matrix convertMatrix( REAL **_Matrix, int n )
{
        Matrix finalMatrix(n,n);

        for( int i = 1; i <= n ; i++ )
        {
                for( int j = 1; j <= n ; j++ )
                {
                        finalMatrix(i,j) = _Matrix[i-1][j-1];
                }
        }
        return finalMatrix;
}

/*! \brief Takes in a REAL * and returns a Vector with the same values 
 *  @param finalMatrix the Matrix in REAL *
 *  @return _Matrix is a Vector 
 */
Vector convertVector( REAL *_Vector, int n )
{
        Vector finalVector(n);

        for( int i = 1; i <= n; i++ )
        {
                finalVector(i) = _Vector[i-1];
        }
        return finalVector;
}

/*! \brief Determines the inverse of a 3x3 matrix
 *  @param _A matrix is a 3x3 matrix to be inverted
 *  @return _Ainv is a 3x3 Matrix
 */
Matrix matrixInv3x3( Matrix _A )
{
        /* Determine cofactor Matrix */
        Matrix _coFactor(3,3);
        Matrix _temp(2,2);
                _temp(_(1,2),_(1,2) ) = _A( _(2,3),_(2,3) );
        _coFactor(1,1) = matrixDet( _temp );
                _temp( _(1,2),1 ) = _A( _(2,3),1 ); 
                _temp(_(1,2),2) = _A( _(2,3),3 ); 
        _coFactor(1,2) = matrixDet( _temp );
                _temp(_(1,2),_(1,2) ) = _A( _(2,3),_(1,2) );
        _coFactor(1,3) = matrixDet( _temp );
        
                _temp( 1,_(1,2) ) = _A( 1,_(2,3) );
                _temp( 2,_(1,2) ) = _A( 3,_(2,3) );
        _coFactor(2,1) = matrixDet( _temp );
                _temp(1,1) = _A(1,1);
                _temp(1,2) = _A(1,3);
                _temp(2,1) = _A(3,1);
                _temp(2,2) = _A(3,3);
        _coFactor(2,2) = matrixDet( _temp );
                _temp( 1,_(1,2) ) = _A( 1,_(1,2) );
                _temp( 2,_(1,2) ) = _A( 3,_(1,2) );
        _coFactor(2,3) = matrixDet( _temp );

                _temp(_(1,2),_(1,2) ) = _A( _(1,2),_(2,3) );
        _coFactor(3,1) = matrixDet( _temp );
                _temp(_(1,2),1) = _A( _(1,2),1 ); 
                _temp(_(1,2),2) = _A( _(1,2),3 ); 
        _coFactor(3,2) = matrixDet( _temp );
                _temp(_(1,2),_(1,2) ) = _A( _(1,2),_(1,2) );
        _coFactor(3,3) = matrixDet( _temp );
        
        /* Determine Determinate of A */
        double det = _A(1,1)*_coFactor(1,1) - _A(2,1)*_coFactor(2,1) + _A(3,1)*_coFactor(3,1);

        /* Apply alternating signs */
        _coFactor(1,2) = -_coFactor(1,2);
        _coFactor(2,1) = -_coFactor(2,1);
        _coFactor(2,3) = -_coFactor(2,3);
        _coFactor(3,2) = -_coFactor(3,2);

        /* Transpose */
        Matrix invPart = ~_coFactor;

        /* Inverse */
        Matrix _Ainv = invPart/det;
        
        return ( _Ainv );
} 

/*! \brief Determines the determinate of a matrix up to 3x3 
 *  @param _A matrix is a nxn matrix where n<=3
 *  @return _det is a double
 */
double matrixDet( Matrix _A )
{
        /* Check Size */
        int rows = _A[1].getIndexCount( );
        double _det;
        Matrix _Aa(2,2);
        switch (rows)
        {
                case 1:
                        _det = _A(1,1);
                        break;
                case 2:
                        _det = matrixDet2x2( _A );
                        break;
                case 3:
                        _Aa(1,1) = _A(2,1); _Aa(1,2) = _A(2,3);
                        _Aa(2,1) = _A(3,1); _Aa(2,2) = _A(3,3);
                        _det = _A(1,1)*matrixDet2x2( _A( _(2,3),_(2,3) ) ) - _A(1,2)*matrixDet2x2( _Aa ) + _A(1,3)*matrixDet2x2( _A( _(2,3), _(1,2) ) );
                        break;
                default:
                        cerr << "Error: Matrix dimension exceeds current ability of det." << endl;
                        _det = 0;
                        break;
        }
        return ( _det );
}

/*! \brief Determines the determinate of a 2x2 matrix 
 *  @param _A matrix is a 2x2 matrix
 *  @return _det is a double
 */
double matrixDet2x2( Matrix _A )
{
        double _det = _A(1,1)*_A(2,2)-_A(1,2)*_A(2,1);
        return ( _det );
}

// Do not change the comments below - they will be added automatically by CVS
/*****************************************************************************
*       $Log: matrixFunctionInterface.cpp,v $
*       Revision 1.3  2007/07/24 08:40:44  jayhawk_hokie
*       Modified matrix norm2 function. Hard coded inverse of 3x3 matrix and det.
*
*       Revision 1.2  2007/06/19 19:30:18  jayhawk_hokie
*       Added function to determine inverse for 3x3 matrix. Added function to determine det(A).
*
*       Revision 1.1  2007/04/13 16:47:52  bleslie
*       Initial submission of the matrix function interface files
*
*       Revision 1.6  2006/07/05 20:59:13  jayhawk_hokie
*       *** empty log message ***
*
*
*
*
******************************************************************************/

---