TAPsMatrixp.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsMatrixp.cpp

Matrixp class is a class for any dimension Matrixp.

Sukitti Punak   (01/19/2006)
Update          (12/16/2009)
******************************************************************************/
#include "TAPsMatrixp.hpp"
// Using Inclusion Model (i.e. definitions are included in declarations)
//                       (this name.cpp is included in name.hpp)
// Each friend is defined directly inside its declaration.

BEGIN_NAMESPACE_TAPs
//=============================================================================
// Constructors and Destructor
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::Matrixp ( int rows, int cols )
    : R( rows ), C( cols )
{
    e = new T[R*C];
    // Identity Matrix
    MakeIdentity();
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::Matrixp ( Matrixp<T> const &M )
    : R( M.R ), C( M.C )
{
    e = new T[R*C];
    for ( int i = 0; i < R*C; ++i )
        e[i] = M.e[i];
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::Matrixp ( int rows, int cols, T val )
    : R( rows ), C( cols )
{
    e = new T[R*C];
    for ( int i = 0; i < R*C; ++i )
        e[i] = val;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::Matrixp ( int rows, int cols, T const a[] )
    : R( rows ), C( cols )
{
    e = new T[R*C];
    for ( int i = 0; i < R*C; ++i )
        e[i] = a[i];
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::~Matrixp ()
{
    delete [] e;
}
//-----------------------------------------------------------------------------
//=============================================================================
// Member Access
//-----------------------------------------------------------------------------
template <typename T>
inline T & Matrixp<T>::operator[] ( int i )
{   return e[i];    }
//-----------------------------------------------------------------------------
template <typename T>
inline T const & Matrixp<T>::operator[] ( int i ) const
{   return e[i];    }
//-----------------------------------------------------------------------------
template <typename T>
inline T & Matrixp<T>::operator() ( int r, int c )
{   return e[r*C + c];  }
//-----------------------------------------------------------------------------
template <typename T>
inline T const & Matrixp<T>::operator() ( int r, int c ) const
{   return e[r*C + c];  }
//-----------------------------------------------------------------------------
//=============================================================================
// Convert to one dimension array
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::operator const T * () const
{   return e;   }
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T>::operator T * ()
{   return e;   }
//-----------------------------------------------------------------------------
//=============================================================================
// Useful Functions
//-----------------------------------------------------------------------------
template <typename T>
std::ostream & Matrixp<T>::Display ( int precision, std::ostream &o ) const
{
    if ( precision < 0 )    precision = 6;
    //--------------------------------------------------------------------
    //o << typeid(*this).name() << "( ";
    o << "Matrixp<" << typeid(T).name() << ">" 
            << "  " << R << "-by-" << C << "\n";
    o.precision( precision );
    int width = precision + 10;
    for ( int r = 0; r < R*C; r+=C ) {
        o << "| ";
        for ( int c = 0; c < C; c++ ) {
            o.width( width );   o << e[r + c];
        }
        o << " |" << endl;
    }
    o.precision( 6 );
    return o;
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::SetAllElements ( T t )
{
    for ( int i = 0; i < R*C; ++i )
        e[i] = t;
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::SetAllElements ( T const a[] )
{
    for ( int i = 0; i < R*C; ++i )
        e[i] = a[i];
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::MakeIdentity ()
{
    // Identity Matrix
    for ( int r = 0; r < R*C; r+=C )
        for ( int c = 0; c < C; ++c )
            e[r + c] = Math<T>::ZERO;
    if ( R <= C ) {
        for ( int r = 0; r < R; ++r )
            e[r*C + r] = Math<T>::ONE;
    }
    else {
        for ( int r = 0; r < C; ++r )
            e[r*C + r] = Math<T>::ONE;
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::MakeDiagonal ( T d )
{
    // Diagonal Matrix
    for ( int r = 0; r < R*C; r+=C )
        for ( int c = 0; c < C; ++c )
            e[r + c] = Math<T>::ZERO;
    for ( int r = 0; r < R; ++r )
            e[r*C + r] = d;
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::MakeDiagonal ( T const a[] )
{
    // Diagonal Matrix
    for ( int r = 0; r < R*C; r+=C )
        for ( int c = 0; c < C; ++c )
            e[r + c] = Math<T>::ZERO;
    for ( int r = 0; r < R; ++r )
            e[r*C + r] = a[r];
}
//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::MakeZero ()
{
    for ( int i = 0; i < R*C; ++i )
        e[i] = Math<T>::ZERO;
}
//-----------------------------------------------------------------------------
template <typename T>
bool Matrixp<T>::IsIdentity () const
{
    int i;
    for ( int r = 0; r < R; ++r ) {
        for ( int c = 0; c < C; ++c ) {
            i = r*C;
            if ( r != c ) {
                if ( e[i + c] != Math<T>::ZERO ) {
                    return false;
                }
            }
            else {
                if ( e[i + c] != Math<T>::ONE ) {
                    return false;
                }
            }
        }
    }
    return true;
}
//-----------------------------------------------------------------------------
template <typename T>
bool Matrixp<T>::IsSymmetric () const
{
    if ( R != C )   return false;   // Not an N-by-N matrix
    for ( int r = 0; r < R; ++r )
        for ( int c = r; c < C; ++c )
            if ( e[r*C + c] != e[c*R + r] )
                return false;
    return true;
}

//-----------------------------------------------------------------------------
template <typename T>
void Matrixp<T>::ChangeSize ( int numOfRow, int numOfCol )
{
    delete [] e;
    R = numOfRow;
    C = numOfCol;
    e = new T[R*C];
    for ( int i = 0; i < R*C; ++i )
        e[i] = T(0);
}

//=============================================================================
// Matrix Operations
//-----------------------------------------------------------------------------
// Transpose this matrix
template <typename T>
Matrixp<T> & Matrixp<T>::Transposed ()
{
    /*
    T temp; int i, j;
    if ( R >= C ) {     // If row >= col
        for ( int r = 0; r < R; ++r ) {
            for ( int c = r; c < C; ++c ) {
                i = r*C + c;
                j = c*R + r;
                temp = e[i];
                e[i] = e[j];
                e[j] = temp;
            }
        }
    }
    else {              // If col > row
        for ( int c = 0; c < C; ++c ) {
            for ( int r = c; r < R; ++r ) {
                i = r*C + c;
                j = c*R + r;
                temp = e[i];
                e[i] = e[j];
                e[j] = temp;
            }
        }
    }
    //*/

    *this = GetTranspose();
    int r = R;
    R = C;
    C = r;

    return *this;
}
//-----------------------------------------------------------------------------
// Get the transpose matrix of this matrix
template <typename T>
Matrixp<T> Matrixp<T>::GetTranspose() const
{
    Matrixp<T> O( C, R );
    int i;
    for ( int r = 0; r < R; ++r ) {
        i = r*C;
        for ( int c = 0; c < C; ++c ) {
            O[c*R + r] = e[i + c];
        }
    }
    return O;
}
/*
//-----------------------------------------------------------------------------
// Inverse the matrix
template <typename T>
inline Matrixp<T> & Matrixp<T>::Inversed ()
{
    *this = (*this).GetInverse();
    return *this;
}
//-----------------------------------------------------------------------------
// Get the matrix inverse
template <typename T>
Matrixp<T> Matrixp<T>::GetInverse () const
{
    T det = GetDeterminant();

    if ( -Math<T>::EPSILON < det && det < Math<T>::EPSILON )
        return Matrixp( Math<T>::INFINITY );

    return Matrixp<T>(   ( e[4]*e[8] - e[7]*e[5] ) / det,
                            -( e[1]*e[8] - e[7]*e[2] ) / det,
                             ( e[1]*e[5] - e[4]*e[2] ) / det,
                            -( e[3]*e[8] - e[6]*e[5] ) / det,
                             ( e[0]*e[8] - e[6]*e[2] ) / det,
                            -( e[0]*e[5] - e[3]*e[2] ) / det,
                             ( e[3]*e[7] - e[6]*e[4] ) / det,
                            -( e[0]*e[7] - e[6]*e[1] ) / det,
                             ( e[0]*e[4] - e[3]*e[1] ) / det    );
}
//-----------------------------------------------------------------------------
template <typename T>
inline T Matrixp<T>::GetDeterminant () const
{
    return      e[0] * ( e[4]*e[8] - e[5]*e[7] ) 
            -   e[3] * ( e[1]*e[8] - e[2]*e[7] )
            +   e[6] * ( e[1]*e[5] - e[2]*e[4] );
}
//*/
//-----------------------------------------------------------------------------
//=============================================================================
// Assignment Overloaded Operator
//--------------------------------------------------------------------------
// Assignment Operator
template <typename T>
inline Matrixp<T> & Matrixp<T>::operator= ( Matrixp<T> const &M )
{
    //assert( R*C == M.R*M.C );
    if ( this != &M ) {
        delete [] e;
        R = M.R;
        C = M.C;
        e = new T[R*C];
        for ( int i = 0; i < R*C; ++i )
            e[i] = M.e[i];
    }
    return *this;
}
//-----------------------------------------------------------------------------
//=============================================================================
// Unary Overloaded Operators
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator- ()
{
    Matrixp<T> O( R, C );
    for ( int i = 0; i < R*C; ++i )
        O.e[i] = -e[i];
    return O;
}
//=============================================================================
// Assign Overloaded Operators
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> & Matrixp<T>::operator+= ( Matrixp<T> const &M )
{
    assert( R == M.R && C == M.C );
    for ( int i = 0; i < R*C; ++i )
        e[i] += M.e[i];
    return *this;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> & Matrixp<T>::operator-= ( Matrixp<T> const &M )
{
    assert( R == M.R && C == M.C );
    for ( int i = 0; i < R*C; ++i )
        e[i] -= M.e[i];
    return *this;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> & Matrixp<T>::operator*= ( Matrixp<T> const &M )
{
    *this = (*this) * M;
    return *this;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> & Matrixp<T>::operator*= ( T s )
{
    for ( int i = 0; i < R*C; ++i )
        e[i] *= s;
    return *this;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> & Matrixp<T>::operator/= ( T s )
{
    for ( int i = 0; i < R*C; ++i )
        e[i] /= s;
    return *this;
}
//-----------------------------------------------------------------------------
//=============================================================================
// Binary Overloaded Operators
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator+ ( Matrixp<T> const &M ) const
{
    assert( R == M.R && C == M.C );
    Matrixp<T> O( R, C );
    for ( int i = 0; i < R*C; ++i )
        O.e[i] = e[i] + M.e[i];
    return O;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator- ( Matrixp<T> const &M ) const
{
    assert( R == M.R && C == M.C );
    Matrixp<T> O( R, C );
    for ( int i = 0; i < R*C; ++i )
        O.e[i] = e[i] - M.e[i];
    return O;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator* ( Matrixp<T> const &M ) const
{
    assert ( C == M.R );    // cols must equals rows
    Matrixp<T> O( R, M.C, Math<T>::ZERO );
    for ( int r = 0; r < R; ++r ) {
        for ( int c = 0; c < M.C; ++c ) {
            for ( int i = 0; i < C; ++i ) {
                O.e[r*M.C + c] += e[r*C + i] * M.e[i*M.C + c];
            }
        }
    }
    return O;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator* ( T s ) const
{
    Matrixp<T> O( R, C );
    for ( int i = 0; i < R*C; ++i )
        O.e[i] = e[i] * s;
    return O;
}
//-----------------------------------------------------------------------------
template <typename T>
Matrixp<T> Matrixp<T>::operator/ ( T s ) const
{
    Matrixp<T> O( R, C );
    for ( int i = 0; i < R*C; ++i )
        O.e[i] = e[i] / s;
    return O;
}
//-----------------------------------------------------------------------------
template <typename T>
inline void Matrixp<T>::MultLeft ( Matrixp<T> const &M )
{
    *this = M * (*this);
}
//-----------------------------------------------------------------------------
template <typename T>
inline void Matrixp<T>::MultRight ( Matrixp<T> const &M )
{
    *this *= M;
}
/*
//-----------------------------------------------------------------------------
// Matrixp * Vectorp
template <typename T>
Vectorp<T> Matrixp<T>::operator* ( Vectorp<T> const &V ) const
{
    Vectorp<T> O( R );
    for ( int r = 0, i = 0; r < R; ++r, i+=C ) {
        for ( int c = 0; c < C; ++c ) {
            O[r] += e[i + c] * V[c];
        }
    }
    return O;
}
//*/
//=============================================================================
// Helper Fn Members    --------------------------------------------------------
//-----------------------------------------------------------------------------
// Check that row and column numbers are in the valid ranges.
template< typename T >
bool Matrixp< T >::CheckRanges( int r, int c ) const
{
    assert( 0 <= r && r < R );
    assert( 0 <= c && c < C );
    return true;
}
//=============================================================================
END_NAMESPACE_TAPs
//345678901234567890123456789012345678901234567890123456789012345678901234567890
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8