matrix.h

一个方便进行矩阵运算的C＋＋函数包

/********************************************************************\
| 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<float>,
  complex<double>, complex<long double>. 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
   [e-mail] 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 <stdlib.h>
#  include <math.h>
#  include <iostream.h>
#  include <string.h>
#else
#  include <cstdlib>
#  include <cmath>
#  include <iostream>
#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 <stdexcept>
#    else
#      include <stdexcpt.h>
#    endif
#  else
#     include <stdexcep>
#  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 <class T>
#  define matrixT  matrix<T>
#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;