www.pudn.com > Gesture[20040824].rar > matrix.h


/********************************************************************\ 
| Matrix.h: C++ matrix template class include file (Ver.1.0)         | 
|                                                                    | 
| Copyright (c) 1997-1998 by Somnath Kundu. (See Copying.lib file.)  | 
\********************************************************************/ 
 
/********************************************************************\ 
  Purpose: This matrix template class defines majority of the matrix 
  operations as overloaded operators or methods. Users of this class 
  have been assumed to be familiar with matrix algebra. I have not 
  defined any specialization of this template here, so all the instances 
  of matrix will be created implicitly by the compiler. The data types 
  tested with this class are float, double, long double, complex, 
  complex, complex. I think there should not be 
  any problem with very large precision floating point classes. 
 
  Since implementation of exception, namespace and template are still 
  not standardized among the various compilers, you may encounter 
  compilation error with some compiler. In that case remove any of 
  the above three features by defining the following macros: 
 
  _NO_NAMESPACE:  Define this macro to remove namespace. 
 
  _NO_EXCEPTION:  Define this macro to remove exception handling 
                  and use old style of error handling using function. 
 
  _NO_TEMPLATE:   If this macro is defined matrix class of double 
                  type will be generated by default. You can also 
                  generate a different type of matrix like float. 
 
  Since all the definitions are also included in this header file as 
  inline function, some compiler may give warning "inline function 
  can't be expanded". You may ignore/disable this warning using compiler 
  switches. All the operators/methods defined in this class have their 
  natural meaning except the followings: 
 
  Operator/Method                          Description 
  ---------------                          ----------- 
   operator ()   :   This function operator can be used as a 
                     two-dimensional subscript operator to get/set 
                     individual matrix elements. 
 
   operator !    :   This operator has been used to calculate inversion 
                     of matrix. 
 
   operator ~    :   This operator has been used to return transpose of 
                     a matrix. 
 
   operator ^    :   It is used calculate power (by a scalar) of a matrix. 
                     When using this operator in a matrix equation, care 
                     must be taken by parenthesizing it because it has 
                     lower precedence than addition, subtraction, 
                     multiplication and division operators. 
 
   operator >>   :   It is used to read matrix from input stream as per 
                     standard C++ stream operators. 
 
   operator <<   :   It is used to write matrix to output stream as per 
                     standard C++ stream operators. 
 
   If you have any suggestion/bug report, please send it to 
   somnath@kagi.com with details such as compiler used with 
   version number, OS and switches used to compile. Since you 
   have got the source code, it would be better if you solve 
   any syntactic problem for compilation yourself, because I 
   may not have access to that specific compiler/OS.  
 
\********************************************************************/ 
 
#ifndef __cplusplus 
#error Must use C++ for the type matrix. 
#endif 
 
#if !defined(__STD_MATRIX_H) 
#define __STD_MATRIX_H 
 
#if defined(__BORLANDC__) 
#pragma option -w-inl -w-pch 
#endif 
 
#if defined(__BORLANDC__) || _MSC_VER <= 1000 
#  include  
#  include  
#  include  
#  include  
#else 
#  include  
#  include  
#  include  
#endif 
 
#if defined(_MSC_VER) && _MSC_VER <= 1000 
#  define _NO_EXCEPTION   // stdexception is not fully suppoted in MSVC++ 4.0 
typedef int bool; 
#  if !defined(false) 
#    define false  0 
#  endif 
#  if !defined(true) 
#    define true   1 
#  endif 
#endif 
 
#if defined(__BORLANDC__) && !defined(__WIN32__) 
#  define _NO_EXCEPTION        // std exception and namespace are not fully 
#  define _NO_NAMESPACE        // supported in 16-bit compiler 
#endif 
 
#if defined(_MSC_VER) && !defined(_WIN32) 
#  define _NO_EXCEPTION 
#endif 
 
#if defined(_NO_EXCEPTION) 
#  define _NO_THROW 
#  define _THROW_MATRIX_ERROR 
#else 
#  if defined(_MSC_VER) 
#    if _MSC_VER >= 1020 
#      include  
#    else 
#      include  
#    endif 
#  else 
#     include  
#  endif 
#  define _NO_THROW               throw () 
#  define _THROW_MATRIX_ERROR     throw (matrix_error) 
#endif 
 
#ifndef __MINMAX_DEFINED 
#  define max(a,b)    (((a) > (b)) ? (a) : (b)) 
#  define min(a,b)    (((a) < (b)) ? (a) : (b)) 
#endif 
 
#if defined(_MSC_VER) && _MSC_VER <= 1020  // MSVC++ 4.0/4.2 does not 
#  define _NO_NAMESPACE                    // support "std" namespace 
#endif 
 
#if !defined(_NO_NAMESPACE) 
using namespace std; 
#endif 
 
#ifndef _NO_NAMESPACE 
namespace math { 
#endif 
 
#if !defined(_NO_EXCEPTION) 
class matrix_error : public logic_error 
{ 
    public: 
        matrix_error (const string& what_arg) : logic_error( what_arg) {} 
}; 
#define REPORT_ERROR(ErrormMsg)  throw matrix_error( ErrormMsg); 
#else 
inline void _matrix_error (const char* pErrMsg) 
{ 
    cout << pErrMsg << endl; 
    exit(1); 
} 
#define REPORT_ERROR(ErrormMsg)  _matrix_error( ErrormMsg); 
#endif 
 
#if !defined(_NO_TEMPLATE) 
#  define MAT_TEMPLATE  template  
#  define matrixT  matrix 
#else 
#  define MAT_TEMPLATE 
#  define matrixT  matrix 
#  ifdef MATRIX_TYPE 
     typedef MATRIX_TYPE T; 
#  else 
     typedef double T; 
#  endif 
#endif   
 
MAT_TEMPLATE 
class matrix 
{ 
private: 
   T **Val; 
   size_t Row, Col, RowSiz, ColSiz; 
 
   void realloc (size_t row, size_t col); 
   int pivot (size_t row); 
 
public: 
   // Constructors 
   matrix (const matrixT& m); 
   matrix (size_t row = 6, size_t col = 6); 
 
   // Destructor 
   ~matrix (); 
 
   // Value extraction method 
   size_t RowNo () { return Row; } 
   size_t ColNO () { return Col; } 
 
   // Subscript operator 
   T& operator () (size_t row, size_t col) _THROW_MATRIX_ERROR; 
 
   // Unary operators 
   matrixT operator + () _NO_THROW { return *this; } 
   matrixT operator - () _NO_THROW; 
 
   // Assignment operators 
   matrixT& operator = (const matrixT& m) _NO_THROW; 
 
   // Combined assignment - calculation operators 
   matrixT& operator += (const matrixT& m) _THROW_MATRIX_ERROR; 
   matrixT& operator -= (const matrixT& m) _THROW_MATRIX_ERROR; 
   matrixT& operator *= (const matrixT& m) _THROW_MATRIX_ERROR; 
   matrixT& operator *= (const T& c) _NO_THROW; 
   matrixT& operator /= (const T& c) _NO_THROW; 
   matrixT& operator ^= (const size_t& pow) _THROW_MATRIX_ERROR; 
 
   // Logical operators 
   friend bool operator == (const matrixT& m1, const matrixT& m2) _NO_THROW; 
   friend bool operator != (const matrixT& m1, const matrixT& m2) _NO_THROW; 
 
   // Calculation operators 
   friend matrixT operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR; 
   friend matrixT operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR; 
   friend matrixT operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR; 
#if !defined(_MSC_VER) || _MSC_VER > 1020 
   friend matrixT operator * (const matrixT& m, const T& no) _NO_THROW; 
   friend matrixT operator * (const T& no, const matrixT& m) _NO_THROW { return (m*no); } 
#endif 
   friend matrixT operator / (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR { return (m1 * !m2); } 
#if !defined(_MSC_VER) || _MSC_VER > 1020 
   friend matrixT operator / (const matrixT& m, const T& no) _NO_THROW { return (m*(1/no)); } 
   friend matrixT operator / (const T& no, const matrixT& m) _THROW_MATRIX_ERROR { return (!m * no); } 
#endif 
   friend matrixT operator ~ (const matrixT& m) _NO_THROW; 
   friend matrixT operator ! (matrixT m) _THROW_MATRIX_ERROR; 
   friend matrixT operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR; 
 
   // Miscellaneous -methods 
   void Null (const size_t& row, const size_t& col) _NO_THROW; 
   void Null () _NO_THROW; 
   void Unit (const size_t& row) _NO_THROW; 
   void Unit () _NO_THROW; 
   void SetSize (size_t row, size_t col) _NO_THROW; 
 
   // Utility methods 
   matrixT Solve (const matrixT& v) const _THROW_MATRIX_ERROR; 
   matrixT Adj () _THROW_MATRIX_ERROR; 
   T Det () _THROW_MATRIX_ERROR; 
   T Norm () _NO_THROW; 
   T Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR; 
   T Cond () _NO_THROW; 
 
   // Type of matrices 
   bool IsSquare () _NO_THROW { return (Row == Col); }  
   bool IsSingular () _NO_THROW; 
   bool IsDiagonal () _NO_THROW; 
   bool IsScalar () _NO_THROW; 
   bool IsUnit () _NO_THROW; 
   bool IsNull () _NO_THROW; 
   bool IsSymmetric () _NO_THROW; 
   bool IsSkewSymmetric () _NO_THROW; 
   bool IsUpperTiangular () _NO_THROW; 
   bool IsLowerTiangular () _NO_THROW; 
 
   // Io-stream operators 
   friend istream& operator >> (istream& i, matrixT& m); 
   friend ostream& operator << (ostream& o, matrixT& m); 
}; 
 
#if defined(_MSC_VER) && _MSC_VER <= 1020 
#  undef  _NO_THROW               // MSVC++ 4.0/4.2 does not support  
#  undef  _THROW_MATRIX_ERROR     // exception specification in definition 
#  define _NO_THROW 
#  define _THROW_MATRIX_ERROR 
#endif 
 
// constructor 
MAT_TEMPLATE inline 
matrixT::matrix (size_t row, size_t col) 
{ 
   RowSiz = Row = row; 
   ColSiz = Col = col; 
   Val = new T* [row]; 
   for (size_t i=0; i < row; i++) 
      Val[i] = new T [col]; 
} 
 
// copy constructor 
MAT_TEMPLATE inline 
matrixT::matrix (const matrixT& m) 
{ 
   RowSiz = Row = m.Row; 
   ColSiz = Col = m.Col; 
 
   Val = new T* [Row]; 
   size_t colsize = Col * sizeof(T); 
 
   for (size_t i=0; i < Row; i++) 
   { 
      Val[i] = new T [Col]; 
      memcpy( Val[i], m.Val[i], colsize); 
   } 
} 
 
 
// destructor 
MAT_TEMPLATE inline 
matrixT::~matrix (void) 
{ 
   for (size_t i=0; i < RowSiz; i++) 
      delete [] Val[i]; 
   delete [] Val; 
} 
 
 
//  reallocation method 
MAT_TEMPLATE inline void  
matrixT::realloc (size_t row, size_t col) 
{ 
   if (row == RowSiz && col == ColSiz) 
   { 
      Row = RowSiz; 
      Col = ColSiz; 
      return; 
   } 
   size_t i; 
   T** Val1 = new T* [row]; 
   for (i=0; i < row; i++) 
      Val1[i] = new T [col]; 
 
   size_t colSize = min(Col,col) * sizeof(T); 
   size_t minRow = min(Row,row); 
    
   for (i=0; i < minRow; i++) 
      memcpy( Val1[i], Val[i], colSize); 
 
  for (i=0; i < RowSiz; i++) 
      delete [] Val[i]; 
  delete [] Val; 
 
  RowSiz = Row = row; 
  ColSiz = Col = col; 
  Val = Val1; 
 
  return; 
} 
 
 
// public method for resizing matrix 
MAT_TEMPLATE inline void 
matrixT::SetSize (size_t row, size_t col) _NO_THROW 
{ 
   size_t i,j; 
   size_t oldRow = Row; 
   size_t oldCol = Col; 
 
   if (row != RowSiz || col != ColSiz) 
      realloc( row, col); 
 
   if (row > oldRow) 
      for (i=oldRow; i < row; i++) 
         for (j=0; j < oldCol; j++) 
            Val[i][j] = T(0); 
 
   if (col > oldCol) 
      for (i=0; i < col; i++) 
         for (j=oldCol; j < col; j++) 
            Val[i][j] = T(0); 
   return; 
} 
 
// subscript operator to get/set individual elements 
MAT_TEMPLATE inline T& 
matrixT::operator () (size_t row, size_t col) _THROW_MATRIX_ERROR 
{ 
   if (row >= Row || col >= Col) 
      REPORT_ERROR( "matrixT::operator(): Index out of range!"); 
   return Val[row][col]; 
} 
 
// input stream function 
MAT_TEMPLATE inline istream& 
operator >> (istream& istrm, matrixT& m) 
{ 
   for (size_t i=0; i < m.Row; i++) 
      for (size_t j=0; j < m.Col; j++) 
         istrm >> m.Val[i][j]; 
   return istrm; 
} 
 
// output stream function 
MAT_TEMPLATE inline ostream& 
operator << (ostream &ostrm, matrixT& m) 
{ 
   for (size_t i=0; i < m.Row; i++) 
   { 
       for (size_t j=0; j < m.Col; j++) 
          cout << m.Val[i][j] << '\t'; 
       cout << endl; 
   } 
   return ostrm; 
} 
 
// assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator = (const matrixT& m) _NO_THROW 
{ 
   if (Row != m.Row || Col != m.Col) 
      realloc( m.Row,m.Col); 
 
   size_t colbyte = m.Col * sizeof(T); 
 
   for (size_t i=0; i < m.Row; i++) 
      memcpy( Val[i], m.Val[i], colbyte); 
 
   return *this; 
} 
 
// logical equal-to operator 
MAT_TEMPLATE inline bool 
operator == (const matrixT& m1, const matrixT& m2) _NO_THROW 
{ 
   bool retVal = false; 
 
   if (m1.Row != m2.Row || m1.Col != m2.Col) 
      return retVal; 
 
   for (size_t i=0; i < m1.Row; i++) 
      for (size_t j=0; j < m1.Col; i++) 
         if (m1.Val[i][j] != m2.Val[i][j]) 
            return retVal; 
 
   return true; 
} 
 
// logical no-equal-to operator 
MAT_TEMPLATE inline bool  
operator != (const matrixT& m1, const matrixT& m2) _NO_THROW 
{ 
    return (m1 == m2) ? false : true; 
} 
 
 
// combined addition and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator += (const matrixT& m) _THROW_MATRIX_ERROR 
{ 
   if (Row != m.Row || Col != m.Col) 
      REPORT_ERROR( "matrixT::operator+= : Inconsistent matrix size in addition!"); 
   for (size_t i=0; i < m.Row; i++) 
      for (size_t j=0; j < m.Col; j++) 
         Val[i][j] += m.Val[i][j]; 
   return *this; 
} 
 
 
// combined subtraction and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator -= (const matrixT& m) _THROW_MATRIX_ERROR 
{ 
   if (Row != m.Row || Col != m.Col) 
      REPORT_ERROR( "matrixT::operator-= : Inconsistent matrix size in subtraction!"); 
 
   for (size_t i=0; i < m.Row; i++) 
      for (size_t j=0; j < m.Col; j++) 
         Val[i][j] -= m.Val[i][j]; 
   return *this; 
} 
 
// combined scalar multiplication and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator *= (const T& c) _NO_THROW 
{ 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         Val[i][j] *= c; 
   return *this; 
} 
 
 
// combined matrix multiplication and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator *= (const matrixT& m) _THROW_MATRIX_ERROR 
{ 
   if (Col != m.Row) 
      REPORT_ERROR( "matrixT::operator*= : Inconsistance matrix size in multiplication!"); 
   *this = *this * m; 
   return *this; 
} 
 
 
// combined scalar division and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator /= (const T& c) _NO_THROW 
{ 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         Val[i][j] /= c; 
 
   return *this; 
} 
 
// combined power and assignment operator 
MAT_TEMPLATE inline matrixT& 
matrixT::operator ^= (const size_t& pow) _THROW_MATRIX_ERROR 
{ 
   for (size_t i=2; i <= pow; i++) 
      *this = *this * *this; 
 
   return *this; 
} 
 
// unary negation operator 
MAT_TEMPLATE inline matrixT 
matrixT::operator - () _NO_THROW 
{ 
   matrixT temp(Row,Col); 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         temp.Val[i][j] = - Val[i][j]; 
 
   return temp; 
} 
 
// binary addition operator 
MAT_TEMPLATE inline matrixT 
operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR 
{ 
   if (m1.Row != m2.Row || m1.Col != m2.Col) 
      REPORT_ERROR( "matrixT::operator+: Inconsistent matrix size in addition!"); 
 
   matrixT temp(m1.Row,m1.Col); 
 
   for (size_t i=0; i < m1.Row; i++) 
      for (size_t j=0; j < m1.Col; j++) 
         temp.Val[i][j] = m1.Val[i][j] + m2.Val[i][j]; 
 
   return temp; 
} 
 
// binary subtraction operator 
MAT_TEMPLATE inline matrixT 
operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR 
{ 
   if (m1.Row != m2.Row || m1.Col != m2.Col) 
     REPORT_ERROR( "matrixT::operator-: Inconsistent matrix size in subtraction!"); 
 
   matrixT temp(m1.Row,m1.Col); 
 
   for (size_t i=0; i < m1.Row; i++) 
      for (size_t j=0; j < m1.Col; j++) 
         temp.Val[i][j] = m1.Val[i][j] - m2.Val[i][j]; 
 
   return temp; 
} 
 
 
// binary scalar multiplication operator 
MAT_TEMPLATE inline matrixT 
operator * (const matrixT& m, const T& no) _NO_THROW 
{ 
   matrixT temp(m.Row,m.Col); 
 
   for (size_t i=0; i < m.Row; i++) 
      for (size_t j=0; j < m.Col; j++) 
         temp.Val[i][j] = no * m.Val[i][j]; 
 
   return temp; 
} 
 
 
// binary matrix multiplication operator 
MAT_TEMPLATE inline matrixT 
operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR 
{ 
   if (m1.Col != m2.Row) 
      REPORT_ERROR( "matrixT::operator*: Inconsistent matrix size in multiplication!"); 
 
   matrixT temp(m1.Row,m2.Col); 
 
   for (size_t i=0; i < m1.Row; i++) 
      for (size_t j=0; j < m2.Col; j++) 
      { 
         temp.Val[i][j] = T(0); 
         for (size_t k=0; k < m1.Col; k++) 
            temp.Val[i][j] += m1.Val[i][k] * m2.Val[k][j]; 
      } 
   return temp; 
} 
 
// binary power operator 
MAT_TEMPLATE inline matrixT 
operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR 
{ 
   matrixT temp(m); 
 
   for (size_t i=2; i <= pow; i++) 
      temp = temp * temp; 
 
   return temp; 
} 
 
 
// unary transpose operator 
MAT_TEMPLATE inline matrixT 
operator  ~ (const matrixT& m) _NO_THROW 
{ 
   matrixT temp(m.Col,m.Row); 
 
   for (size_t i=0; i < m.Row; i++) 
      for (size_t j=0; j < m.Col; j++) 
	 temp.Val[j][i] = m.Val[i][j]; 
 
   return temp; 
} 
 
 
// unary inversion operator 
MAT_TEMPLATE inline matrixT 
operator ! (matrixT m) _THROW_MATRIX_ERROR 
{ 
   size_t i,j,k; 
   T a1,a2,*rowptr; 
 
   if (m.Row != m.Col) 
      REPORT_ERROR( "matrixT::operator!: Inversion of a non-square matrix"); 
 
   matrixT temp(m.Row,m.Col); 
 
   temp.Unit(); 
   for (k=0; k < m.Row; k++) 
   { 
      int indx = m.pivot(k); 
      if (indx == -1) 
         REPORT_ERROR( "matrixT::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 
MAT_TEMPLATE inline matrixT 
matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR 
{ 
   size_t i,j,k; 
   T a1; 
 
   if (!(Row == Col && Col == v.Row)) 
      REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!"); 
 
   matrixT 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) 
         REPORT_ERROR( "matrixT::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]; 
      } 
   } 
   matrixT 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 
MAT_TEMPLATE inline void 
matrixT::Null (const size_t& row, const size_t& col) _NO_THROW 
{ 
   if (row != Row || col != Col) 
      realloc( row,col); 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         Val[i][j] = T(0); 
   return; 
} 
 
// set zero to all elements of this matrix 
MAT_TEMPLATE inline void 
matrixT::Null() _NO_THROW 
{ 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
 	 Val[i][j] = T(0); 
   return; 
} 
 
// set this matrix to unity 
MAT_TEMPLATE inline void 
matrixT::Unit (const size_t& row) _NO_THROW 
{ 
   if (row != Row || row != Col) 
      realloc( row, row); 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         Val[i][j] = i == j ? T(1) : T(0); 
  return; 
} 
 
// set this matrix to unity 
MAT_TEMPLATE inline void 
matrixT::Unit () _NO_THROW 
{ 
   size_t row = min(Row,Col); 
   Row = Col = row; 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         Val[i][j] = i == j ? T(1) : T(0); 
   return; 
} 
 
 
// private partial pivoting method 
MAT_TEMPLATE inline int 
matrixT::pivot (size_t row) 
{ 
  int k = int(row); 
  double amax,temp; 
 
  amax = -1; 
  for (size_t i=row; i < Row; i++) 
     if ( (temp = abs( Val[i][row])) > amax && temp != 0.0) 
     { 
       amax = temp; 
       k = i; 
     } 
  if (Val[k][row] == T(0)) 
     return -1; 
  if (k != int(row)) 
  { 
     T* rowptr = Val[k]; 
     Val[k] = Val[row]; 
     Val[row] = rowptr; 
     return k; 
  } 
  return 0; 
} 
 
 
// calculate the determinant of a matrix 
MAT_TEMPLATE inline T 
matrixT::Det () _THROW_MATRIX_ERROR 
{ 
   size_t i,j,k; 
   T piv,detVal = T(1); 
 
   if (Row != Col) 
      REPORT_ERROR( "matrixT::Det(): Determinant a non-square matrix!"); 
 
   matrixT 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 
MAT_TEMPLATE inline T 
matrixT::Norm () _NO_THROW 
{ 
   T retVal = T(0); 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         retVal += Val[i][j] * Val[i][j]; 
   retVal = sqrt( retVal); 
 
   return retVal; 
} 
 
 
// calculate the condition number of a matrix 
MAT_TEMPLATE inline T 
matrixT::Cond () _NO_THROW 
{ 
   matrixT inv(Row,Col); 
 
   inv = ! (*this); 
   T retVal = Norm() * inv.Norm(); 
 
   return retVal; 
} 
 
 
// calculate the cofactor of a matrix for a given element 
MAT_TEMPLATE inline T 
matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR 
{ 
   size_t i,i1,j,j1; 
 
   if (Row != Col) 
      REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!"); 
 
   if (row > Row || col > Col) 
      REPORT_ERROR( "matrixT::Cofact(): Index out of range!"); 
 
   matrixT temp (Row-1,Col-1); 
 
   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++; 
   } 
   T  cof = temp.Det(); 
   if ((row+col)%2 == 1) 
      cof = -cof; 
 
   return cof; 
} 
 
 
// calculate adjoin of a matrix 
MAT_TEMPLATE inline matrixT 
matrixT::Adj () _THROW_MATRIX_ERROR 
{ 
   if (Row != Col) 
      REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix."); 
 
   matrixT temp(Row,Col); 
 
   for (size_t i=0; i < Row; i++) 
      for (size_t j=0; j < Col; j++) 
         temp.Val[i][j] = Cofact(i,j); 
 
   temp = ~temp; 
   return temp; 
} 
 
// Determine if the matrix is singular 
MAT_TEMPLATE inline bool 
matrixT::IsSingular () _NO_THROW 
{ 
   if (Row != Col) 
      return false; 
   return (Det() == T(0)); 
} 
 
// Determine if the matrix is diagonal 
MAT_TEMPLATE inline bool 
matrixT::IsDiagonal () _NO_THROW 
{ 
   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 
MAT_TEMPLATE inline bool 
matrixT::IsScalar () _NO_THROW 
{ 
   if (!IsDiagonal()) 
     return false; 
   T 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 
MAT_TEMPLATE inline bool 
matrixT::IsUnit () _NO_THROW 
{ 
   if (IsScalar() && Val[0][0] == T(1)) 
     return true; 
   return false; 
} 
 
// Determine if this is a null matrix 
MAT_TEMPLATE inline bool 
matrixT::IsNull () _NO_THROW 
{ 
   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 
MAT_TEMPLATE inline bool 
matrixT::IsSymmetric () _NO_THROW 
{ 
   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 skew-symmetric 
MAT_TEMPLATE inline bool 
matrixT::IsSkewSymmetric () _NO_THROW 
{ 
   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 
MAT_TEMPLATE inline bool 
matrixT::IsUpperTiangular () _NO_THROW 
{ 
   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 
MAT_TEMPLATE inline bool 
matrixT::IsLowerTiangular () _NO_THROW 
{ 
   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; 
} 
 
#ifndef _NO_NAMESPACE 
}  
#endif 
 
#endif //__STD_MATRIX_H