matrix.h

Go to the documentation of this file.
00001 
00002 // Matrix TCL Lite v1.13
00003 // Copyright (c) 1997-2002 Techsoft Pvt. Ltd. (See License.Txt file.)
00004 //
00005 // Matrix.h: Matrix C++ template class include file 
00006 // Web: http://www.techsoftpl.com/matrix/
00007 // Email: matrix@techsoftpl.com
00008 //
00009 
00011 // Installation:
00012 //
00013 // Copy this "matrix.h" file into include directory of your compiler.
00014 //
00015 
00017 // Note: This matrix template class defines majority of the matrix
00018 // operations as overloaded operators or methods. It is assumed that
00019 // users of this class is familiar with matrix algebra. We have not
00020 // defined any specialization of this template here, so all the instances
00021 // of matrix will be created implicitly by the compiler. The data types
00022 // tested with this class are float, double, long double, complex<float>,
00023 // complex<double> and complex<long double>. Note that this class is not 
00024 // optimized for performance.
00025 //
00026 // Since implementation of exception, namespace and template are still
00027 // not standardized among the various (mainly old) compilers, you may 
00028 // encounter compilation error with some compilers. In that case remove 
00029 // any of the above three features by defining the following macros:
00030 //
00031 //  _NO_NAMESPACE:  Define this macro to remove namespace support.
00032 //
00033 //  _NO_EXCEPTION:  Define this macro to remove exception handling
00034 //                  and use old style of error handling using function.
00035 //
00036 //  _NO_TEMPLATE:   If this macro is defined matrix class of double
00037 //                  type will be generated by default. You can also
00038 //                  generate a different type of matrix like float.
00039 //
00040 //  _SGI_BROKEN_STL: For SGI C++ v.7.2.1 compiler.
00041 //
00042 //  Since all the definitions are also included in this header file as
00043 //  inline function, some compiler may give warning "inline function
00044 //  can't be expanded". You may ignore/disable this warning using compiler
00045 //  switches. All the operators/methods defined in this class have their
00046 //  natural meaning except the followings:
00047 //
00048 //  Operator/Method                          Description
00049 //  ---------------                          -----------
00050 //   operator ()   :   This function operator can be used as a
00051 //                     two-dimensional subscript operator to get/set
00052 //                     individual matrix elements.
00053 //
00054 //   operator !    :   This operator has been used to calculate inversion
00055 //                     of matrix.
00056 //
00057 //   operator ~    :   This operator has been used to return transpose of
00058 //                     a matrix.
00059 //
00060 //   operator ^    :   It is used calculate power (by a scalar) of a matrix.
00061 //                     When using this operator in a matrix equation, care
00062 //                     must be taken by parenthesizing it because it has
00063 //                     lower precedence than addition, subtraction,
00064 //                     multiplication and division operators.
00065 //
00066 //   operator >>   :   It is used to read matrix from input stream as per
00067 //                     standard C++ stream operators.
00068 //
00069 //   operator <<   :   It is used to write matrix to output stream as per
00070 //                     standard C++ stream operators.
00071 //
00072 // Note that professional version of this package, Matrix TCL Pro 2.11
00073 // is optimized for performance and supports many more matrix operations.
00074 // It is available from our web site at <http://www.techsoftpl.com/matrix/>.
00075 //
00076 
00077 #ifndef __cplusplus
00078 #error Must use C++ for the type matrix.
00079 #endif
00080 
00081 #if !defined(__STD_MATRIX_H)
00082 #define __STD_MATRIX_H
00083 
00085 // First deal with various shortcomings and incompatibilities of
00086 // various (mainly old) versions of popular compilers available.
00087 //
00088 
00089 #if defined(__BORLANDC__)
00090 #pragma option -w-inl -w-pch
00091 #endif
00092 
00093 #if ( defined(__BORLANDC__) || _MSC_VER <= 1000 ) && !defined( __GNUG__ )
00094 #  include <stdio.h>
00095 #  include <stdlib.h>
00096 #  include <math.h>
00097 #  include <iostream.h>
00098 #  include <string.h>
00099 #else
00100 #  include <string.h>
00101 #  include <cmath>
00102 #  include <cstdio>
00103 #  include <cstdlib>
00104 #  include <string>
00105 #  include <iostream>
00106 #endif
00107 
00108 #if defined(_MSC_VER) && _MSC_VER <= 1000
00109 #  define _NO_EXCEPTION        // stdexception is not fully supported in MSVC++ 4.0
00110 typedef int bool;
00111 #  if !defined(false)
00112 #    define false  0
00113 #  endif
00114 #  if !defined(true)
00115 #    define true   1
00116 #  endif
00117 #endif
00118 
00119 #if defined(__BORLANDC__) && !defined(__WIN32__)
00120 #  define _NO_EXCEPTION        // std exception and namespace are not fully
00121 #  define _NO_NAMESPACE        // supported in 16-bit compiler
00122 #endif
00123 
00124 #if defined(_MSC_VER) && !defined(_WIN32)
00125 #  define _NO_EXCEPTION
00126 #endif
00127 
00128 #if defined(_NO_EXCEPTION)
00129 #  define _NO_THROW
00130 #  define _THROW_MATRIX_ERROR
00131 #else
00132 #  if defined(_MSC_VER)
00133 #    if _MSC_VER >= 1020
00134 #      include <stdexcept>
00135 #    else
00136 #      include <stdexcpt.h>
00137 #    endif
00138 #  elif defined(__MWERKS__)
00139 #      include <stdexcept>
00140 #  elif (__GNUC__ >= 2 || (__GNUC__ == 2 && __GNUC_MINOR__ >= 8))
00141 #     include <stdexcept>
00142 #  else
00143 #     include <stdexcep>
00144 #  endif
00145 #  define _NO_THROW               throw ()
00146 #  define _THROW_MATRIX_ERROR     throw (matrix_error)
00147 #endif
00148 
00149 #ifndef __MINMAX_DEFINED
00150 #  define max(a,b)    (((a) > (b)) ? (a) : (b))
00151 #  define min(a,b)    (((a) < (b)) ? (a) : (b))
00152 #endif
00153 
00154 #if defined(_MSC_VER)
00155 #undef _MSC_EXTENSIONS      // To include overloaded abs function definitions!
00156 #endif
00157 
00158 #if ( defined(__BORLANDC__) || _MSC_VER ) && !defined( __GNUG__ ) 
00159 inline float abs (float v) { return (float)fabs( v); } 
00160 inline double abs (double v) { return fabs( v); } 
00161 inline long double abs (long double v) { return fabsl( v); }
00162 #endif
00163 
00164 #if defined(__GNUG__) || defined(__MWERKS__) || (defined(__BORLANDC__) && (__BORLANDC__ >= 0x540))
00165 #define FRIEND_FUN_TEMPLATE <>
00166 #else
00167 #define FRIEND_FUN_TEMPLATE
00168 #endif
00169 
00170 #if defined(_MSC_VER) && _MSC_VER <= 1020   // MSVC++ 4.0/4.2 does not
00171 #  define _NO_NAMESPACE                     // support "std" namespace
00172 #endif
00173 
00174 #if !defined(_NO_NAMESPACE)
00175 #if defined( _SGI_BROKEN_STL )              // For SGI C++ v.7.2.1 compiler
00176 namespace std { }
00177 #endif
00178 using namespace std;
00179 #endif
00180 
00181 #ifndef _NO_NAMESPACE
00182 namespace math {
00183 #endif
00184 
00185 #if !defined(_NO_EXCEPTION)
00186 class matrix_error : public logic_error
00187 {
00188     public:
00189     matrix_error (const string& what_arg) : logic_error( what_arg) { }
00190 };
00191 #define REPORT_ERROR(ErrormMsg)  throw matrix_error( ErrormMsg);
00192 #else
00193 inline void _matrix_error (const char* pErrMsg)
00194 {
00195     cerr << pErrMsg << endl;
00196     exit(1);
00197 }
00198 #define REPORT_ERROR(ErrormMsg)  _matrix_error( ErrormMsg);
00199 #endif
00200 
00201 #if !defined(_NO_TEMPLATE)
00202 #  define MAT_TEMPLATE  template <class T>
00203 #  define matrixT  matrix<T>
00204 #else
00205 #  define MAT_TEMPLATE
00206 #  define matrixT  matrix
00207 #  ifdef MATRIX_TYPE
00208      typedef MATRIX_TYPE T;
00209 #  else
00210      typedef double T;
00211 #  endif
00212 #endif
00213 
00214 
00215 MAT_TEMPLATE
00216 class matrix
00217 {
00218 public:
00219    // Constructors
00220    matrix (const matrixT& m);
00221    matrix (size_t row = 6, size_t col = 6);
00222 
00223    // Destructor
00224    ~matrix ();
00225 
00226    // Assignment operators
00227    matrixT& operator = (const matrixT& m) _NO_THROW;
00228 
00229    // Value extraction method
00230    size_t RowNo () const { return _m->Row; }
00231    size_t ColNo () const { return _m->Col; }
00232 
00233    // Subscript operator
00234    T& operator () (size_t row, size_t col) _THROW_MATRIX_ERROR;
00235    T  operator () (size_t row, size_t col) const _THROW_MATRIX_ERROR;
00236 
00237    // Unary operators
00238    matrixT operator + () _NO_THROW { return *this; }
00239    matrixT operator - () _NO_THROW;
00240 
00241    // Combined assignment - calculation operators
00242    matrixT& operator += (const matrixT& m) _THROW_MATRIX_ERROR;
00243    matrixT& operator -= (const matrixT& m) _THROW_MATRIX_ERROR;
00244    matrixT& operator *= (const matrixT& m) _THROW_MATRIX_ERROR;
00245    matrixT& operator *= (const T& c) _NO_THROW;
00246    matrixT& operator /= (const T& c) _NO_THROW;
00247    matrixT& operator ^= (const size_t& pow) _THROW_MATRIX_ERROR;
00248 
00249    // Miscellaneous -methods
00250    void Null (const size_t& row, const size_t& col) _NO_THROW;
00251    void Null () _NO_THROW;
00252    void Unit (const size_t& row) _NO_THROW;
00253    void Unit () _NO_THROW;
00254    void SetSize (size_t row, size_t col) _NO_THROW;
00255 
00256    // Utility methods
00257    matrixT Solve (const matrixT& v) const _THROW_MATRIX_ERROR;
00258    matrixT Adj () _THROW_MATRIX_ERROR;
00259    matrixT Inv () _THROW_MATRIX_ERROR;
00260    T Det () const _THROW_MATRIX_ERROR;
00261    T Norm () _NO_THROW;
00262    T Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR;
00263    T Cond () _NO_THROW;
00264 
00265    // Type of matrices
00266    bool IsSquare () _NO_THROW { return (_m->Row == _m->Col); } 
00267    bool IsSingular () _NO_THROW;
00268    bool IsDiagonal () _NO_THROW;
00269    bool IsScalar () _NO_THROW;
00270    bool IsUnit () _NO_THROW;
00271    bool IsNull () _NO_THROW;
00272    bool IsSymmetric () _NO_THROW;
00273    bool IsSkewSymmetric () _NO_THROW;
00274    bool IsUpperTriangular () _NO_THROW;
00275    bool IsLowerTriangular () _NO_THROW;
00276 
00277 private:
00278     struct base_mat
00279     {
00280     T **Val;
00281     size_t Row, Col, RowSiz, ColSiz;
00282     int Refcnt;
00283 
00284     base_mat (size_t row, size_t col, T** v)
00285  {
00286         Row = row; RowSiz = row;
00287         Col = col; ColSiz = col;
00288         Refcnt = 1;
00289 
00290         Val = new T* [row];
00291         size_t rowlen = col * sizeof(T);
00292 
00293         for (size_t i=0; i < row; i++)
00294         {
00295         Val[i] = new T [col];
00296         if (v) memcpy( Val[i], v[i], rowlen);
00297         }
00298     }
00299     ~base_mat ()
00300  {
00301         for (size_t i=0; i < RowSiz; i++)
00302         delete [] Val[i];
00303         delete [] Val;
00304     }
00305     };
00306     base_mat *_m;
00307 
00308     void clone ();
00309     void realloc (size_t row, size_t col);
00310     int pivot (size_t row);
00311 };
00312 
00313 #if defined(_MSC_VER) && _MSC_VER <= 1020
00314 #  undef  _NO_THROW               // MSVC++ 4.0/4.2 does not support 
00315 #  undef  _THROW_MATRIX_ERROR     // exception specification in definition
00316 #  define _NO_THROW
00317 #  define _THROW_MATRIX_ERROR
00318 #endif
00319 
00320 // constructor
00321 MAT_TEMPLATE inline
00322 matrixT::matrix (size_t row, size_t col)
00323 {
00324   _m = new base_mat( row, col, 0);
00325 }
00326 
00327 // copy constructor
00328 MAT_TEMPLATE inline
00329 matrixT::matrix (const matrixT& m)
00330 {
00331     _m = m._m;
00332     _m->Refcnt++;
00333 }
00334 
00335 // Internal copy constructor
00336 MAT_TEMPLATE inline void
00337 matrixT::clone ()
00338 {
00339     _m->Refcnt--;
00340     _m = new base_mat( _m->Row, _m->Col, _m->Val);
00341 }
00342 
00343 // destructor
00344 MAT_TEMPLATE inline
00345 matrixT::~matrix ()
00346 {
00347    if (--_m->Refcnt == 0) delete _m;
00348 }
00349 
00350 // assignment operator
00351 MAT_TEMPLATE inline matrixT&
00352 matrixT::operator = (const matrixT& m) _NO_THROW
00353 {
00354     m._m->Refcnt++;
00355     if (--_m->Refcnt == 0) delete _m;
00356     _m = m._m;
00357     return *this;
00358 }
00359 
00360 //  reallocation method
00361 MAT_TEMPLATE inline void 
00362 matrixT::realloc (size_t row, size_t col)
00363 {
00364    if (row == _m->RowSiz && col == _m->ColSiz)
00365    {
00366       _m->Row = _m->RowSiz;
00367       _m->Col = _m->ColSiz;
00368       return;
00369    }
00370 
00371    base_mat *m1 = new base_mat( row, col, NULL);
00372    size_t colSize = min(_m->Col,col) * sizeof(T);
00373    size_t minRow = min(_m->Row,row);
00374 
00375    for (size_t i=0; i < minRow; i++)
00376       memcpy( m1->Val[i], _m->Val[i], colSize);
00377 
00378    if (--_m->Refcnt == 0) 
00379        delete _m;
00380    _m = m1;
00381 
00382    return;
00383 }
00384 
00385 // public method for resizing matrix
00386 MAT_TEMPLATE inline void
00387 matrixT::SetSize (size_t row, size_t col) _NO_THROW
00388 {
00389    size_t i,j;
00390    size_t oldRow = _m->Row;
00391    size_t oldCol = _m->Col;
00392 
00393    if (row != _m->RowSiz || col != _m->ColSiz)
00394       realloc( row, col);
00395 
00396    for (i=oldRow; i < row; i++)
00397       for (j=0; j < col; j++)
00398      _m->Val[i][j] = T(0);
00399 
00400    for (i=0; i < row; i++)                      
00401       for (j=oldCol; j < col; j++)
00402      _m->Val[i][j] = T(0);
00403 
00404    return;
00405 }
00406 
00407 // subscript operator to get/set individual elements
00408 MAT_TEMPLATE inline T&
00409 matrixT::operator () (size_t row, size_t col) _THROW_MATRIX_ERROR
00410 {
00411    if (row >= _m->Row || col >= _m->Col)
00412       REPORT_ERROR( "matrixT::operator(): Index out of range!");
00413    if (_m->Refcnt > 1) clone();
00414    return _m->Val[row][col];
00415 }
00416 
00417 // subscript operator to get/set individual elements
00418 MAT_TEMPLATE inline T
00419 matrixT::operator () (size_t row, size_t col) const _THROW_MATRIX_ERROR
00420 {
00421    if (row >= _m->Row || col >= _m->Col)
00422       REPORT_ERROR( "matrixT::operator(): Index out of range!");
00423    return _m->Val[row][col];
00424 }
00425 
00426 // input stream function
00427 MAT_TEMPLATE inline istream&
00428 operator >> (istream& istrm, matrixT& m)
00429 {
00430    for (size_t i=0; i < m.RowNo(); i++)
00431       for (size_t j=0; j < m.ColNo(); j++)
00432       {
00433          T x;
00434          istrm >> x;
00435          m(i,j) = x;
00436       }
00437    return istrm;
00438 }
00439 
00440 // output stream function
00441 MAT_TEMPLATE inline ostream&
00442 operator << (ostream& ostrm, const matrixT& m)
00443 {
00444    for (size_t i=0; i < m.RowNo(); i++)
00445    {
00446       for (size_t j=0; j < m.ColNo(); j++)
00447       {
00448          T x = m(i,j);
00449          ostrm << x << '\t';
00450       }
00451       ostrm << endl;
00452    }
00453    return ostrm;
00454 }
00455 
00456 
00457 // logical equal-to operator
00458 MAT_TEMPLATE inline bool
00459 operator == (const matrixT& m1, const matrixT& m2) _NO_THROW
00460 {
00461    if (m1.RowNo() != m2.RowNo() || m1.ColNo() != m2.ColNo())
00462       return false;
00463 
00464    for (size_t i=0; i < m1.RowNo(); i++)
00465       for (size_t j=0; j < m1.ColNo(); j++)
00466           if (m1(i,j) != m2(i,j))
00467              return false;
00468 
00469    return true;
00470 }
00471 
00472 // logical no-equal-to operator
00473 MAT_TEMPLATE inline bool
00474 operator != (const matrixT& m1, const matrixT& m2) _NO_THROW
00475 {
00476     return (m1 == m2) ? false : true;
00477 }
00478 
00479 // combined addition and assignment operator
00480 MAT_TEMPLATE inline matrixT&
00481 matrixT::operator += (const matrixT& m) _THROW_MATRIX_ERROR
00482 {
00483    if (_m->Row != m._m->Row || _m->Col != m._m->Col)
00484       REPORT_ERROR( "matrixT::operator+= : Inconsistent matrix sizes in addition!");
00485    if (_m->Refcnt > 1) clone();
00486    for (size_t i=0; i < m._m->Row; i++)
00487       for (size_t j=0; j < m._m->Col; j++)
00488      _m->Val[i][j] += m._m->Val[i][j];
00489    return *this;
00490 }
00491 
00492 // combined subtraction and assignment operator
00493 MAT_TEMPLATE inline matrixT&
00494 matrixT::operator -= (const matrixT& m) _THROW_MATRIX_ERROR
00495 {
00496    if (_m->Row != m._m->Row || _m->Col != m._m->Col)
00497       REPORT_ERROR( "matrixT::operator-= : Inconsistent matrix sizes in subtraction!");
00498    if (_m->Refcnt > 1) clone();
00499    for (size_t i=0; i < m._m->Row; i++)
00500       for (size_t j=0; j < m._m->Col; j++)
00501      _m->Val[i][j] -= m._m->Val[i][j];
00502    return *this;
00503 }
00504 
00505 // combined scalar multiplication and assignment operator
00506 MAT_TEMPLATE inline matrixT&
00507 matrixT::operator *= (const T& c) _NO_THROW
00508 {
00509     if (_m->Refcnt > 1) clone();
00510     for (size_t i=0; i < _m->Row; i++)
00511     for (size_t j=0; j < _m->Col; j++)
00512         _m->Val[i][j] *= c;
00513     return *this;
00514 }
00515 
00516 // combined matrix multiplication and assignment operator
00517 MAT_TEMPLATE inline matrixT&
00518 matrixT::operator *= (const matrixT& m) _THROW_MATRIX_ERROR
00519 {
00520    if (_m->Col != m._m->Row)
00521       REPORT_ERROR( "matrixT::operator*= : Inconsistent matrix sizes in multiplication!");
00522 
00523    matrixT temp(_m->Row,m._m->Col);
00524 
00525    for (size_t i=0; i < _m->Row; i++)
00526       for (size_t j=0; j < m._m->Col; j++)
00527       {
00528          temp._m->Val[i][j] = T(0);
00529          for (size_t k=0; k < _m->Col; k++)
00530             temp._m->Val[i][j] += _m->Val[i][k] * m._m->Val[k][j];
00531       }
00532    *this = temp;
00533 
00534    return *this;
00535 }
00536 
00537 // combined scalar division and assignment operator
00538 MAT_TEMPLATE inline matrixT&
00539 matrixT::operator /= (const T& c) _NO_THROW
00540 {
00541     if (_m->Refcnt > 1) clone();
00542     for (size_t i=0; i < _m->Row; i++)
00543     for (size_t j=0; j < _m->Col; j++)
00544         _m->Val[i][j] /= c;
00545 
00546     return *this;
00547 }
00548 
00549 // combined power and assignment operator
00550 MAT_TEMPLATE inline matrixT&
00551 matrixT::operator ^= (const size_t& pow) _THROW_MATRIX_ERROR
00552 {
00553     matrixT temp(*this);
00554 
00555     for (size_t i=2; i <= pow; i++)
00556       *this = *this * temp;
00557 
00558     return *this;
00559 }
00560 
00561 // unary negation operator
00562 MAT_TEMPLATE inline matrixT
00563 matrixT::operator - () _NO_THROW
00564 {
00565    matrixT temp(_m->Row,_m->Col);
00566 
00567    for (size_t i=0; i < _m->Row; i++)
00568       for (size_t j=0; j < _m->Col; j++)
00569      temp._m->Val[i][j] = - _m->Val[i][j];
00570 
00571    return temp;
00572 }
00573 
00574 // binary addition operator
00575 MAT_TEMPLATE inline matrixT
00576 operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
00577 {
00578    matrixT temp = m1;
00579    temp += m2;
00580    return temp;
00581 }
00582 
00583 // binary subtraction operator
00584 MAT_TEMPLATE inline matrixT
00585 operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
00586 {
00587    matrixT temp = m1;
00588    temp -= m2;
00589    return temp;
00590 }
00591 
00592 // binary scalar multiplication operator
00593 MAT_TEMPLATE inline matrixT
00594 operator * (const matrixT& m, const T& no) _NO_THROW
00595 {
00596    matrixT temp = m;
00597    temp *= no;
00598    return temp;
00599 }
00600 
00601 
00602 // binary scalar multiplication operator
00603 MAT_TEMPLATE inline matrixT
00604 operator * (const T& no, const matrixT& m) _NO_THROW
00605 {
00606    return (m * no);
00607 }
00608 
00609 // binary matrix element multiplication operator
00610 MAT_TEMPLATE inline matrixT
00611 operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
00612 {
00613     matrixT temp = m1;
00614     temp *= m2;
00615     return temp;
00616 }
00617 
00618 
00619 // binary scalar division operator
00620 MAT_TEMPLATE inline matrixT
00621 operator / (const matrixT& m, const T& no) _NO_THROW
00622 {
00623     return (m * (T(1) / no));
00624 }
00625 
00626 
00627 // binary scalar division operator
00628 MAT_TEMPLATE inline matrixT
00629 operator / (const T& no, const matrixT& m) _THROW_MATRIX_ERROR
00630 {
00631     return (!m * no);
00632 }
00633 
00634 // binary matrix division operator
00635 MAT_TEMPLATE inline matrixT
00636 operator / (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
00637 {
00638     return (m1 * !m2);
00639 }
00640 
00641 // binary power operator
00642 MAT_TEMPLATE inline matrixT
00643 operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR
00644 {
00645    matrixT temp = m;
00646    temp ^= pow;
00647    return temp;
00648 }
00649 
00650 // unary transpose operator
00651 MAT_TEMPLATE inline matrixT
00652 operator ~ (const matrixT& m) _NO_THROW
00653 {
00654    matrixT temp(m.ColNo(),m.RowNo());
00655 
00656    for (size_t i=0; i < m.RowNo(); i++)
00657       for (size_t j=0; j < m.ColNo(); j++)
00658       {
00659          T x = m(i,j);
00660           temp(j,i) = x;
00661       }
00662    return temp;
00663 }
00664 
00665 // unary inversion operator
00666 MAT_TEMPLATE inline matrixT
00667 operator ! (const matrixT m) _THROW_MATRIX_ERROR
00668 {
00669    matrixT temp = m;
00670    return temp.Inv();
00671 }
00672 
00673 // inversion function
00674 MAT_TEMPLATE inline matrixT
00675 matrixT::Inv () _THROW_MATRIX_ERROR
00676 {
00677    size_t i,j,k;
00678    T a1,a2,*rowptr;
00679 
00680    if (_m->Row != _m->Col)
00681       REPORT_ERROR( "matrixT::operator!: Inversion of a non-square matrix");
00682 
00683    matrixT temp(_m->Row,_m->Col);
00684    if (_m->Refcnt > 1) clone();
00685 
00686 
00687    temp.Unit();
00688    for (k=0; k < _m->Row; k++)
00689    {
00690       int indx = pivot(k);
00691       if (indx == -1)
00692           REPORT_ERROR( "matrixT::operator!: Inversion of a singular matrix");
00693 
00694       if (indx != 0)
00695       {
00696           rowptr = temp._m->Val[k];
00697           temp._m->Val[k] = temp._m->Val[indx];
00698           temp._m->Val[indx] = rowptr;
00699       }
00700       a1 = _m->Val[k][k];
00701       for (j=0; j < _m->Row; j++)
00702       {
00703           _m->Val[k][j] /= a1;
00704           temp._m->Val[k][j] /= a1;
00705       }
00706       for (i=0; i < _m->Row; i++)
00707        if (i != k)
00708        {
00709           a2 = _m->Val[i][k];
00710           for (j=0; j < _m->Row; j++)
00711           {
00712              _m->Val[i][j] -= a2 * _m->Val[k][j];
00713              temp._m->Val[i][j] -= a2 * temp._m->Val[k][j];
00714           }
00715        }
00716    }
00717    return temp;
00718 }
00719 
00720 // solve simultaneous equation
00721 MAT_TEMPLATE inline matrixT
00722 matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR
00723 {
00724    size_t i,j,k;
00725    T a1;
00726 
00727    if (!(_m->Row == _m->Col && _m->Col == v._m->Row))
00728       REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!");
00729 
00730    matrixT temp(_m->Row,_m->Col+v._m->Col);
00731    for (i=0; i < _m->Row; i++)
00732    {
00733       for (j=0; j < _m->Col; j++)
00734      temp._m->Val[i][j] = _m->Val[i][j];
00735       for (k=0; k < v._m->Col; k++)
00736      temp._m->Val[i][_m->Col+k] = v._m->Val[i][k];
00737    }
00738    for (k=0; k < _m->Row; k++)
00739    {
00740       int indx = temp.pivot(k);
00741       if (indx == -1)
00742      REPORT_ERROR( "matrixT::Solve(): Singular matrix!");
00743 
00744       a1 = temp._m->Val[k][k];
00745       for (j=k; j < temp._m->Col; j++)
00746      temp._m->Val[k][j] /= a1;
00747 
00748       for (i=k+1; i < _m->Row; i++)
00749       {
00750      a1 = temp._m->Val[i][k];
00751      for (j=k; j < temp._m->Col; j++)
00752        temp._m->Val[i][j] -= a1 * temp._m->Val[k][j];
00753       }
00754    }
00755    matrixT s(v._m->Row,v._m->Col);
00756    for (k=0; k < v._m->Col; k++)
00757       for (int m=int(_m->Row)-1; m >= 0; m--)
00758       {
00759      s._m->Val[m][k] = temp._m->Val[m][_m->Col+k];
00760      for (j=m+1; j < _m->Col; j++)
00761         s._m->Val[m][k] -= temp._m->Val[m][j] * s._m->Val[j][k];
00762       }
00763    return s;
00764 }
00765 
00766 // set zero to all elements of this matrix
00767 MAT_TEMPLATE inline void
00768 matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
00769 {
00770     if (row != _m->Row || col != _m->Col)
00771     realloc( row,col);
00772 
00773     if (_m->Refcnt > 1) 
00774     clone();
00775 
00776     for (size_t i=0; i < _m->Row; i++)
00777     for (size_t j=0; j < _m->Col; j++)
00778         _m->Val[i][j] = T(0);
00779     return;
00780 }
00781 
00782 // set zero to all elements of this matrix
00783 MAT_TEMPLATE inline void
00784 matrixT::Null() _NO_THROW
00785 {
00786     if (_m->Refcnt > 1) clone();   
00787     for (size_t i=0; i < _m->Row; i++)
00788     for (size_t j=0; j < _m->Col; j++)
00789         _m->Val[i][j] = T(0);
00790     return;
00791 }
00792 
00793 // set this matrix to unity
00794 MAT_TEMPLATE inline void
00795 matrixT::Unit (const size_t& row) _NO_THROW
00796 {
00797     if (row != _m->Row || row != _m->Col)
00798     realloc( row, row);
00799     
00800     if (_m->Refcnt > 1) 
00801     clone();
00802 
00803     for (size_t i=0; i < _m->Row; i++)
00804     for (size_t j=0; j < _m->Col; j++)
00805         _m->Val[i][j] = i == j ? T(1) : T(0);
00806     return;
00807 }
00808 
00809 // set this matrix to unity
00810 MAT_TEMPLATE inline void
00811 matrixT::Unit () _NO_THROW
00812 {
00813     if (_m->Refcnt > 1) clone();   
00814     size_t row = min(_m->Row,_m->Col);
00815     _m->Row = _m->Col = row;
00816 
00817     for (size_t i=0; i < _m->Row; i++)
00818     for (size_t j=0; j < _m->Col; j++)
00819         _m->Val[i][j] = i == j ? T(1) : T(0);
00820     return;
00821 }
00822 
00823 // private partial pivoting method
00824 MAT_TEMPLATE inline int
00825 matrixT::pivot (size_t row)
00826 {
00827   int k = int(row);
00828   double amax,temp;
00829 
00830   amax = -1;
00831   for (size_t i=row; i < _m->Row; i++)
00832     if ( (temp = abs( _m->Val[i][row])) > amax && temp != 0.0)
00833      {
00834        amax = temp;
00835        k = i;
00836      }
00837   if (_m->Val[k][row] == T(0))
00838      return -1;
00839   if (k != int(row))
00840   {
00841      T* rowptr = _m->Val[k];
00842      _m->Val[k] = _m->Val[row];
00843      _m->Val[row] = rowptr;
00844      return k;
00845   }
00846   return 0;
00847 }
00848 
00849 // calculate the determinant of a matrix
00850 MAT_TEMPLATE inline T
00851 matrixT::Det () const _THROW_MATRIX_ERROR
00852 {
00853    size_t i,j,k;
00854    T piv,detVal = T(1);
00855 
00856    if (_m->Row != _m->Col)
00857       REPORT_ERROR( "matrixT::Det(): Determinant a non-square matrix!");
00858    
00859    matrixT temp(*this);
00860    if (temp._m->Refcnt > 1) temp.clone();
00861 
00862    for (k=0; k < _m->Row; k++)
00863    {
00864       int indx = temp.pivot(k);
00865       if (indx == -1)
00866      return 0;
00867       if (indx != 0)
00868      detVal = - detVal;
00869       detVal = detVal * temp._m->Val[k][k];
00870       for (i=k+1; i < _m->Row; i++)
00871       {
00872      piv = temp._m->Val[i][k] / temp._m->Val[k][k];
00873      for (j=k+1; j < _m->Row; j++)
00874         temp._m->Val[i][j] -= piv * temp._m->Val[k][j];
00875       }
00876    }
00877    return detVal;
00878 }
00879 
00880 // calculate the norm of a matrix
00881 MAT_TEMPLATE inline T
00882 matrixT::Norm () _NO_THROW
00883 {
00884    T retVal = T(0);
00885 
00886    for (size_t i=0; i < _m->Row; i++)
00887       for (size_t j=0; j < _m->Col; j++)
00888      retVal += _m->Val[i][j] * _m->Val[i][j];
00889    retVal = sqrt( retVal);
00890 
00891    return retVal;
00892 }
00893 
00894 // calculate the condition number of a matrix
00895 MAT_TEMPLATE inline T
00896 matrixT::Cond () _NO_THROW
00897 {
00898    matrixT inv = ! (*this);
00899    return (Norm() * inv.Norm());
00900 }
00901 
00902 // calculate the cofactor of a matrix for a given element
00903 MAT_TEMPLATE inline T
00904 matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR
00905 {
00906    size_t i,i1,j,j1;
00907 
00908    if (_m->Row != _m->Col)
00909       REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!");
00910 
00911    if (row > _m->Row || col > _m->Col)
00912       REPORT_ERROR( "matrixT::Cofact(): Index out of range!");
00913 
00914    matrixT temp (_m->Row-1,_m->Col-1);
00915 
00916    for (i=i1=0; i < _m->Row; i++)
00917    {
00918       if (i == row)
00919     continue;
00920       for (j=j1=0; j < _m->Col; j++)
00921       {
00922      if (j == col)
00923         continue;
00924      temp._m->Val[i1][j1] = _m->Val[i][j];
00925      j1++;
00926       }
00927       i1++;
00928    }
00929    T  cof = temp.Det();
00930    if ((row+col)%2 == 1)
00931       cof = -cof;
00932 
00933    return cof;
00934 }
00935 
00936 
00937 // calculate adjoin of a matrix
00938 MAT_TEMPLATE inline matrixT
00939 matrixT::Adj () _THROW_MATRIX_ERROR
00940 {
00941    if (_m->Row != _m->Col)
00942       REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix.");
00943 
00944    matrixT temp(_m->Row,_m->Col);
00945 
00946    for (size_t i=0; i < _m->Row; i++)
00947       for (size_t j=0; j < _m->Col; j++)
00948       temp._m->Val[j][i] = Cofact(i,j);
00949    return temp;
00950 }
00951 
00952 // Determine if the matrix is singular
00953 MAT_TEMPLATE inline bool
00954 matrixT::IsSingular () _NO_THROW
00955 {
00956    if (_m->Row != _m->Col)
00957       return false;
00958    return (Det() == T(0));
00959 }
00960 
00961 // Determine if the matrix is diagonal
00962 MAT_TEMPLATE inline bool
00963 matrixT::IsDiagonal () _NO_THROW
00964 {
00965    if (_m->Row != _m->Col)
00966       return false;
00967    for (size_t i=0; i < _m->Row; i++)
00968      for (size_t j=0; j < _m->Col; j++)
00969     if (i != j && _m->Val[i][j] != T(0))
00970       return false;
00971    return true;
00972 }
00973 
00974 // Determine if the matrix is scalar
00975 MAT_TEMPLATE inline bool
00976 matrixT::IsScalar () _NO_THROW
00977 {
00978    if (!IsDiagonal())
00979      return false;
00980    T v = _m->Val[0][0];
00981    for (size_t i=1; i < _m->Row; i++)
00982      if (_m->Val[i][i] != v)
00983     return false;
00984    return true;
00985 }
00986 
00987 // Determine if the matrix is a unit matrix
00988 MAT_TEMPLATE inline bool
00989 matrixT::IsUnit () _NO_THROW
00990 {
00991    if (IsScalar() && _m->Val[0][0] == T(1))
00992      return true;
00993    return false;
00994 }
00995 
00996 // Determine if this is a null matrix
00997 MAT_TEMPLATE inline bool
00998 matrixT::IsNull () _NO_THROW
00999 {
01000    for (size_t i=0; i < _m->Row; i++)
01001       for (size_t j=0; j < _m->Col; j++)
01002      if (_m->Val[i][j] != T(0))
01003         return false;
01004    return true;
01005 }
01006 
01007 // Determine if the matrix is symmetric
01008 MAT_TEMPLATE inline bool
01009 matrixT::IsSymmetric () _NO_THROW
01010 {
01011    if (_m->Row != _m->Col)
01012       return false;
01013    for (size_t i=0; i < _m->Row; i++)
01014       for (size_t j=0; j < _m->Col; j++)
01015      if (_m->Val[i][j] != _m->Val[j][i])
01016         return false;
01017    return true;
01018 }
01019        
01020 // Determine if the matrix is skew-symmetric
01021 MAT_TEMPLATE inline bool
01022 matrixT::IsSkewSymmetric () _NO_THROW
01023 {
01024    if (_m->Row != _m->Col)
01025       return false;
01026    for (size_t i=0; i < _m->Row; i++)
01027       for (size_t j=0; j < _m->Col; j++)
01028      if (_m->Val[i][j] != -_m->Val[j][i])
01029         return false;
01030    return true;
01031 }
01032    
01033 // Determine if the matrix is upper triangular
01034 MAT_TEMPLATE inline bool
01035 matrixT::IsUpperTriangular () _NO_THROW
01036 {
01037    if (_m->Row != _m->Col)
01038       return false;
01039    for (size_t i=1; i < _m->Row; i++)
01040       for (size_t j=0; j < i-1; j++)
01041      if (_m->Val[i][j] != T(0))
01042         return false;
01043    return true;
01044 }
01045 
01046 // Determine if the matrix is lower triangular
01047 MAT_TEMPLATE inline bool
01048 matrixT::IsLowerTriangular () _NO_THROW
01049 {
01050    if (_m->Row != _m->Col)
01051       return false;
01052 
01053    for (size_t j=1; j < _m->Col; j++)
01054       for (size_t i=0; i < j-1; i++)
01055      if (_m->Val[i][j] != T(0))
01056         return false;
01057 
01058    return true;
01059 }
01060 
01061 #ifndef _NO_NAMESPACE
01062 } 
01063 #endif
01064 
01065 #endif //__STD_MATRIX_H

Generated on Thu Sep 17 23:14:16 2009 for CSL by  doxygen 1.5.8