MboMatrix.h

Go to the documentation of this file.

/*****************************************************************************
* Copyright (c) 2001 - 2008 Marcus Boerger.  All rights reserved.
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
* ========================================================================= */

/* ------------------------------------------------------------------------ */
/* Name:      MboMatrix.h
*
* Requires:
* - Mbo.h
*/
/* ------------------------------------------------------------------------ */

#ifndef _MBOMATRIX_H_
#define _MBOMATRIX_H_

#include <exception>

#if !defined(__MINMAX_DEFINED) && !defined(max) && !defined(min)
#  define max(a,b)    (((a) > (b)) ? (a) : (b))
#  define min(a,b)    (((a) < (b)) ? (a) : (b))
#endif

namespace mbo
{

/* ------------------------------------------------------------------------ */
class MatrixException
: public std::exception
{
    MatrixException(_IN const char *szWhat)
    : std::exception(szWhat)
    {
    }
};

/* ------------------------------------------------------------------------ */
template <class _Ty, size_t _Cols, size_t _Rows>
class Matrix
{
    public:
    // Constructors
    Matrix (const Matrix<_Ty, _Cols, _Rows> & m);

    size_t numRows () { return numRows; }
    size_t numCols () { return numCols; }

    // Subscript operator
    _Ty& operator () (size_t row, size_t col) throw(MatrixException);

    // Unary operators
    Matrix<_Ty, _Cols, _Rows> operator + ()  { return *this; }
    Matrix<_Ty, _Cols, _Rows> operator - ();

    // Assignment operators
    Matrix<_Ty, _Cols, _Rows> & operator = (const Matrix<_Ty, _Cols, _Rows> & m);

    // Combined assignment - calculation operators
    Matrix<_Ty, _Cols, _Rows> & operator += (const Matrix<_Ty, _Cols, _Rows> & m) throw(MatrixException);
    Matrix<_Ty, _Cols, _Rows> & operator -= (const Matrix<_Ty, _Cols, _Rows> & m) throw(MatrixException);
    Matrix<_Ty, _Cols, _Rows> & operator *= (const Matrix<_Ty, _Cols, _Rows> & m) throw(MatrixException);
    Matrix<_Ty, _Cols, _Rows> & operator *= (const _Ty& c);
    Matrix<_Ty, _Cols, _Rows> & operator /= (const _Ty& c);
    Matrix<_Ty, _Cols, _Rows> & operator ^= (const size_t& pow) throw(MatrixException);

    // Logical operators
    friend bool operator == (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2);
    friend bool operator != (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2);

    // Calculation operators
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator + (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2) throw(MatrixException);
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator - (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2) throw(MatrixException);
    friend template<class _Ty, size_t _Cols1Rows2, size_t _Rows1, size_t _Cols2>
        Matrix<_Ty, _Cols2, _Rows1> operator * (const Matrix<_Ty, _Cols1Rows2, _Rows1> & lhs, const Matrix<_Ty, _Cols2, _Cols1Rows2> & rhs)
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator * (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2) throw(MatrixException);
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator * (const Matrix<_Ty, _Cols, _Rows> & m, const T& no);
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator * (const T& no, const Matrix<_Ty, _Cols, _Rows> & m)  { return (m*no); }
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator / (const Matrix<_Ty, _Cols, _Rows> & m1, const Matrix<_Ty, _Cols, _Rows> & m2) throw(MatrixException) { return (m1 * !m2); }
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator / (const Matrix<_Ty, _Cols, _Rows> & m, const T& no)  { return (m*(1/no)); }
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator / (const T& no, const Matrix<_Ty, _Cols, _Rows> & m) throw(MatrixException) { return (!m * no); }
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator ~ (const Matrix<_Ty, _Cols, _Rows> & m);
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator ! (Matrix<_Ty, _Cols, _Rows> m) throw(MatrixException);
    friend template<class _Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows> operator ^ (const Matrix<_Ty, _Cols, _Rows> & m, const size_t& pow) throw(MatrixException);

    // Miscellaneous -methods
    void Null (const size_t& row, const size_t& col);
    void Null ();
    void Unit (const size_t& row);
    void Unit ();

    // Utility methods
    Matrix<_Ty, _Cols, _Rows> Solve(const Matrix<_Ty, _Cols, _Rows> & v) const throw(MatrixException);
    Matrix<_Ty, _Cols, _Rows> Adj() throw(MatrixException);
    _Ty Det() throw(MatrixException);
    _Ty Norm();
    _Ty Cofactor(size_t row, size_t col) throw(MatrixException);
    _Ty Cond();

    // Type of matrices
    bool IsSquare ()  { return (Row == Col); }
    bool IsSingular ();
    bool IsDiagonal ();
    bool IsScalar ();
    bool IsUnit ();
    bool IsNull ();
    bool IsSymmetric ();
    bool IsSkewSymmetric ();
    bool IsUpperTiangular ();
    bool IsLowerTiangular ();

    // Io-stream operators
    friend istream& operator >> (istream& i, Matrix<_Ty, _Cols, _Rows> & m);
    friend ostream& operator << (ostream& o, Matrix<_Ty, _Cols, _Rows> & m);

    private:
    _Ty   Val[_Cols][_Rows];

    int pivot(size_t row);
};

template<class _Ty, size_t _Cols>
class Vector: Matrix<_Ty, _Cols, 1>
{
};

// constructor
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows>::Matrix ()
{
}

// copy constructor
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows>::Matrix (const Matrix<_Ty, _Cols, _Rows> & m)
    : Val(m.Val)
{
}

// subscript operator to get/set individual elements
template<class _Ty, size_t _Cols, size_t _Rows>
inline _Ty& Matrix<_Ty, _Cols, _Rows>::operator () (size_t nRow, size_t nCol) throw(MatrixException)
{
    if (nRow >= _Rows || nCol >= _Cols)
    {
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::operator(): Index out of range!");
    }
    return Val[row][col];
}

// input stream function
template<class _Ty, size_t _Cols, size_t _Rows>
istream& operator >> (istream& in, Matrix<_Ty, _Cols, _Rows> & m)
{
    for (size_t i=0; i < m.Row; i++)
    {
        for (size_t j=0; j < m.Col; j++)
        {
            in >> m.Val[i][j];
        }
    }
    return in;
}

// output stream function
template<class _Ty, size_t _Cols, size_t _Rows>
ostream& operator << (ostream &out, Matrix<_Ty, _Cols, _Rows> & m)
{
    for (size_t i=0; i < m.Row; i++)
    {
        for (size_t j=0; j < m.Col; j++)
        {
            out << m.Val[i][j] << '\t';
        }
        out << endl;
    }
    return out;
}

// assignment operator
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator = (const Matrix<_Ty, _Cols, _Rows> & m)
{
    *Val = *m.Val;
    return *this;
}

// logical equal-to operator
template<class _Ty, size_t _Cols, size_t _Rows>
bool operator == (const Matrix<_Ty, _Cols, _Rows> & lhs, const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; i++)
        {
            if (lhs1.Val[i][j] != rhs2.Val[i][j])
                return false;
            }
        }
    }
    return true;
}

// logical no-equal-to operator
template<class _Ty, size_t _Cols, size_t _Rows>
bool operator != (const Matrix<_Ty, _Cols, _Rows> & lhs, const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    return !(lhs == rhs);
}


// combined addition and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator += (const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] += rhs.Val[i][j];
        }
    }
    return *this;
}


// combined subtraction and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows> 
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator -= (const Matrix<_Ty, _Cols, _Rows> & m)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] -= rhs.Val[i][j];
        }
    }
    return *this;
}

// combined scalar multiplication and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows> 
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator *= (const _Ty& c)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] *= c;
        }
    }
    return *this;
}


// combined Matrix multiplication and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator *= (const Matrix<_Ty, _Cols, _Rows> & m)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] *= rhs.Val[i][j];
        }
    }
    return *this;
}


// combined scalar division and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator /= (const T& c)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] /= c;
        }
    }
    return *this;
}

// combined power and assignment operator
template<class _Ty, size_t _Cols, size_t _Rows> 
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::operator ^= (const size_t& c) throw(MatrixException)
{
    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            Val[i][j] ^= c;
        }
    }
    return *this;
}

// unary negation operator
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows>::operator - ()
{
    Matrix<_Ty, _Cols, _Rows> ret();

    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            ret.Val[i][j] = - Val[i][j];
        }
    }
    return ret;
}

// binary addition operator
template<class _Ty, size_t _Cols, size_t _Rows> 
Matrix<_Ty, _Cols, _Rows> operator + (const Matrix<_Ty, _Cols, _Rows> & lhs, const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    Matrix<_Ty, _Cols, _Rows>  ret();

    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            ret.Val[i][j] = lhs.Val[i][j] + rhs.Val[i][j];
        }
    }
    return ret;
}

// binary subtraction operator
template<class _Ty, size_t _Cols, size_t _Rows>
Matrix<_Ty, _Cols, _Rows> operator - (const Matrix<_Ty, _Cols, _Rows> & lhs, const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    Matrix<_Ty, _Cols, _Rows> ret();

    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            ret.Val[i][j] = lhs.Val[i][j] - rhs.Val[i][j];
        }
    }
    return ret;
}


// binary scalar multiplication operator
template<class _Ty, size_t _Cols, size_t _Rows>
Matrix<_Ty, _Cols, _Rows> operator * (const Matrix<_Ty, _Cols, _Rows> & lhs, const _Ty& c)
{
    Matrix<_Ty, _Cols, _Rows> ret();

    for (size_t i=0; i < m.Row; i++)
    {
        for (size_t j=0; j < m.Col; j++)
        {
            ret.Val[i][j] = lhs.Val[i][j] * c;
        }
    }
    return ret;
}


// binary Matrix multiplication operator
template<class _Ty, size_t _Cols1Rows2, size_t _Rows1, size_t _Cols2>
Matrix<_Ty, _Cols, _Rows> operator * (const Matrix<_Ty, _Cols1Rows2, _Rows1> & lhs, const Matrix<_Ty, _Cols2, _Cols1Rows2> & rhs)
{
    Matrix<_Ty, _Cols2, _Rows1> ret();

    for (size_t i=0; i < _Rows1; i++)
    {
        for (size_t j=0; j < _Cols2; j++)
        {
            ret.Val[i][j] = _Ty(0);
            for (size_t k=0; k < _Cols1Rows2; k++)
            {
                ret.Val[i][j] += lhs.Val[i][k] * rhs[k][j];
            }
        }
    }
    return ret;
}

// binary power operator
template<class _Ty, size_t _Cols, size_t _Rows>
Matrix<_Ty, _Cols, _Rows> operator ^ (const Matrix<_Ty, _Cols, _Rows> & lhs, const size_t & rhs)
{
    assert(rhs > 0); // need matrix E to begin with?

    Matrix<_Ty, _Cols, _Rows> temp(m);

    for (size_t i=2; i <= pow; i++)
    {
        ret = ret * lhs;
    }

    return ret;
}


// unary transpose operator
template<class _Ty, size_t _Cols, size_t _Rows>
Matrix<_Ty, _Cols, _Rows> operator  ~ (const Matrix<_Ty, _Cols, _Rows> & rhs)
{
    Matrix<_Ty, _Cols, _Rows> ret();

    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            ret.Val[i][j] = rhs.Val[j][i];
        }
    }
    return ret;
}


// unary inversion operator
template<class _Ty, size_t _Cols, size_t _Rows>
Matrix<_Ty, _Cols, _Rows> operator ! (Matrix<_Ty, _Cols, _Rows> rhs)
{xxx
    size_t i,j,k;
    T a1,a2,*rowptr;

    Matrix<_Ty, _Cols, _Rows> temp(m.Row,m.Col);

    temp.Unit();
    for (k=0; k < m.Row; k++)
    {
        int indx = m.pivot(k);
        if (indx == -1)
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::operator!: Inversion of a singular Matrix");

        if (indx != 0)
        {
            rowptr = temp.Val[k];
            temp.Val[k] = temp.Val[indx];
            temp.Val[indx] = rowptr;
        }
        a1 = m.Val[k][k];
        for (j=0; j < m.Row; j++)
        {
            m.Val[k][j] /= a1;
            temp.Val[k][j] /= a1;
        }
        for (i=0; i < m.Row; i++)
        if (i != k)
        {
            a2 = m.Val[i][k];
            for (j=0; j < m.Row; j++)
            {
                m.Val[i][j] -= a2 * m.Val[k][j];
                temp.Val[i][j] -= a2 * temp.Val[k][j];
            }
        }
    }
    return temp;
}


// solve simultaneous equation
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows>::Solve(const Matrix<_Ty, _Cols, _Rows> & v) const throw(MatrixException)
{xxx
    size_t i, j, k;
    _Ty a1;

    Matrix<_Ty, _Cols, _Rows> temp(Row, Col + v.Col);
    for (i=0; i < Row; i++)
    {
        for (j=0; j < Col; j++)
        temp.Val[i][j] = Val[i][j];
        for (k=0; k < v.Col; k++)
        temp.Val[i][Col+k] = v.Val[i][k];
    }
    for (k=0; k < Row; k++)
    {
        int indx = temp.pivot(k);
        if (indx == -1)
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::Solve(): Singular Matrix!");

        a1 = temp.Val[k][k];
        for (j=k; j < temp.Col; j++)
        temp.Val[k][j] /= a1;

        for (i=k+1; i < Row; i++)
        {
            a1 = temp.Val[i][k];
            for (j=k; j < temp.Col; j++)
            temp.Val[i][j] -= a1 * temp.Val[k][j];
        }
    }
    Matrix<_Ty, _Cols, _Rows> s(v.Row,v.Col);
    for (k=0; k < v.Col; k++)
    for (int m=int(Row)-1; m >= 0; m--)
    {
        s.Val[m][k] = temp.Val[m][Col+k];
        for (j=m+1; j < Col; j++)
        s.Val[m][k] -= temp.Val[m][j] * s.Val[j][k];
    }
    return s;
}


// set zero to all elements of this Matrix
template<class _Ty, size_t _Cols, size_t _Rows> 
inline Matrix<_Ty, _Cols, _Rows> & Matrix<_Ty, _Cols, _Rows>::Null()
{
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            Val[i][j] = _Ty(0);
        }
    }
    return *this;
}

// set this Matrix to unity
template<class _Ty, size_t _Cols, size_t _Rows>
inline void Matrix<_Ty, _Cols, _Rows>::Unity()
{
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            Val[i][j] = i == j ? _Ty(1) : _Ty(0);
        }
    }
    return;
}

// private partial pivoting method
template<class _Ty, size_t _Cols, size_t _Rows> 
inline int Matrix<_Ty, _Cols, _Rows>::pivot(size_t row)
{
    size_t k = row;
    double amax = -1.0, temp;

    for (size_t i=row; i < _Rows; i++)
    {
        if ((temp = abs(Val[i][row])) > amax && temp != 0.0)
        {
            amax = temp;
            k = i;
        }
    }
    if (Val[k][row] == _Ty(0))
    {
        return -1;
    }
    if (k != row)
    {
        _Ty* rowptr = Val[k];
        Val[k] = Val[row];
        Val[row] = rowptr;
        int ret = static_cast<int>(k);
        if (ret == -1)
        {
            throw MatrixException("Matrix<_Ty, _Cols, _Rows::pivot(row): return value == -1");
        }
        return ret;
    }
    return 0;
}


// calculate the determinant of a Matrix
template<class _Ty, size_t _Cols, size_t _Rows>
inline _Ty Matrix<_Ty, _Cols, _Rows>::Det()
{
    size_t i, j, k;
    T piv, detVal = T(1);

    Matrix<_Ty, _Cols, _Rows> temp(*this);

    for (k = 0; k < Row; k++)
    {
        int indx = temp.pivot(k);
        if (indx == -1)
        {
            return 0;
        }
        if (indx != 0)
        {
            detVal = - detVal;
        }
        detVal = detVal * temp.Val[k][k];
        for (i=k+1; i < Row; i++)
        {
            piv = temp.Val[i][k] / temp.Val[k][k];
            for (j=k+1; j < Row; j++)
            {
                temp.Val[i][j] -= piv * temp.Val[k][j];
            }
        }
    }
    return detVal;
}

// calculate the norm of a Matrix
template<class _Ty, size_t _Cols, size_t _Rows>
inline _Ty Matrix<_Ty, _Cols, _Rows>::Norm()
{
    _Ty ret = _Ty(0);

    for (size_t i=0; i < _Rows; i++)
    {
        for (size_t j=0; j < _Cols; j++)
        {
            ret += Val[i][j] * Val[i][j];
        }
    }
    ret = sqrt(ret);
    return ret;
}


// calculate the condition number of a Matrix
template<class _Ty, size_t _Cols, size_t _Rows> T
inline Matrix<_Ty, _Cols, _Rows>::Cond()
{
    Matrix<_Ty, _Cols, _Rows> inv();

    inv = ! (*this);
    _Ty ret = Norm() * inv.Norm();

    return retVal;
}


// calculate the cofactor of a Matrix for a given element
template<class _Ty, size_t _Cols, size_t _Rows>
inline _Ty Matrix<_Ty, _Cols, _Rows>::Cofactor(size_t row, size_t col) throw(MatrixException)
{xxx
    size_t i, i1, j, j1;

    if (Row != Col)
    {
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::Cofactor(): Cofactor of a non-square Matrix!");
    }

    if (row > _Rows || col > _Cols)
    {
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::Cofactor(): Index out of range!");
    }

    Matrix<_Ty, _Cols-1, _Rows-1> temp();

    for (i=i1=0; i < Row; i++)
    {
        if (i == row)
        {
            continue;
        }
        for (j=j1=0; j < Col; j++)
        {
            if (j == col)
            {
                continue;
            }
            temp.Val[i1][j1] = Val[i][j];
            j1++;
        }
        i1++;
    }
    return ((row+col)%2 == 1) ? -temp.Det() : temp.Det();
}


// calculate adjoin of a Matrix
template<class _Ty, size_t _Cols, size_t _Rows>
inline Matrix<_Ty, _Cols, _Rows> Matrix<_Ty, _Cols, _Rows>::Adjoin()
{xxx
    if (Row != Col)
    {
        throw MatrixException("Matrix<_Ty, _Cols, _Rows>::Adjoin(): Adjoin of a non-square Matrix.");
    }

    Matrix<_Ty, _Cols, _Rows> ret();

    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            ret.Val[i][j] = Cofact(i,j);
        }
    }
    ret = ~ret;
    return ret;
}

// Determine if the Matrix is singular
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsSingular()
{
    if (Row != Col)
    {
        return false;
    }
    return (Det() == T(0));
}

// Determine if the Matrix is diagonal
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsDiagonal()
{
    if (Row != Col)
    {
        return false;
    }
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            if (i != j && Val[i][j] != T(0))
            {
                return false;
            }
        }
    }
    return true;
}

// Determine if the Matrix is scalar
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsScalar()
{
    if (!IsDiagonal())
    {
        return false;
    }
    _Ty v = Val[0][0];
    for (size_t i=1; i < Row; i++)
    {
        if (Val[i][i] != v)
        {
            return false;
        }
    }
    return true;
}

// Determine if the Matrix is a unit Matrix
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsUnit ()
{
    return (IsScalar() && Val[0][0] == T(1)) ? true : false;
}

// Determine if this is a null Matrix
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsNull()
{
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            if (Val[i][j] != T(0))
            {
                return false;
            }
        }
    }
    return true;
}

// Determine if the Matrix is symmetric
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsSymmetric()
{
    if (_Rows != _Cols)
    {
        return false;
    }
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            if (Val[i][j] != Val[j][i])
            {
                return false;
            }
        }
    }
    return true;
}

// Determine if the Matrix is skew-symmetric
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsSkewSymmetric()
{
    if (Row != Col)
    {
        return false;
    }
    for (size_t i=0; i < Row; i++)
    {
        for (size_t j=0; j < Col; j++)
        {
            if (Val[i][j] != -Val[j][i])
            {
                return false;
            }
        }
    }
    return true;
}

// Determine if the Matrix is upper triangular
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsUpperTiangular()
{
    if (Row != Col)
    {
        return false;
    }
    for (size_t i=1; i < Row; i++)
    {
        for (size_t j=0; j < i-1; j++)
        {
            if (Val[i][j] != T(0))
            {
                return false;
            }
        }
    }
    return true;
}

// Determine if the Matrix is lower triangular
template<class _Ty, size_t _Cols, size_t _Rows>
inline bool Matrix<_Ty, _Cols, _Rows>::IsLowerTiangular ()
{
    if (Row != Col)
    {
        return false;
    }
    for (size_t j=1; j < Col; j++)
    {
        for (size_t i=0; i < j-1; i++)
        {
            if (Val[i][j] != T(0))
            {
                return false;
            }
        }
    }
    return true;
}

} // namepace Mbo

#endif // _MBOMATRIX_H_

---