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 <cmath>
00101 #  include <cstdio>
00102 #  include <cstdlib>
00103 #  include <cstring>
00104 #  include <iostream>
00105 #endif
00106 
00107 #if defined(_MSC_VER) && _MSC_VER <= 1000
00108 #  define _NO_EXCEPTION        // stdexception is not fully supported in MSVC++ 4.0
00109 typedef int bool;
00110 #  if !defined(false)
00111 #    define false  0
00112 #  endif
00113 #  if !defined(true)
00114 #    define true   1
00115 #  endif
00116 #endif
00117 
00118 #if defined(__BORLANDC__) && !defined(__WIN32__)
00119 #  define _NO_EXCEPTION        // std exception and namespace are not fully
00120 #  define _NO_NAMESPACE        // supported in 16-bit compiler
00121 #endif
00122 
00123 #if defined(_MSC_VER) && !defined(_WIN32)
00124 #  define _NO_EXCEPTION
00125 #endif
00126 
00127 #if defined(_NO_EXCEPTION)
00128 #  define _NO_THROW
00129 #  define _THROW_MATRIX_ERROR
00130 #else
00131 #  if defined(_MSC_VER)
00132 #    if _MSC_VER >= 1020
00133 #      include <stdexcept>
00134 #    else
00135 #      include <stdexcpt.h>
00136 #    endif
00137 #  elif defined(__MWERKS__)
00138 #      include <stdexcept>
00139 #  elif (__GNUC__ >= 2 || (__GNUC__ == 2 && __GNUC_MINOR__ >= 8))
00140 #     include <stdexcept>
00141 #  else
00142 #     include <stdexcep>
00143 #  endif
00144 #  define _NO_THROW               throw ()
00145 #  define _THROW_MATRIX_ERROR     throw (matrix_error)
00146 #endif
00147 
00148 #ifndef __MINMAX_DEFINED
00149 //#  define max(a,b)    (((a) > (b)) ? (a) : (b))
00150 //#  define min(a,b)    (((a) < (b)) ? (a) : (b))
00151 #endif
00152 
00153 #if defined(_MSC_VER)
00154 #undef _MSC_EXTENSIONS      // To include overloaded abs function definitions!
00155 #endif
00156 
00157 #if ( defined(__BORLANDC__) || _MSC_VER ) && !defined( __GNUG__ )
00158 inline float abs (float v) { return (float)fabs( v); }
00159 inline double abs (double v) { return fabs( v); }
00160 inline long double abs (long double v) { return fabsl( v); }
00161 #endif
00162 
00163 #if defined(__GNUG__) || defined(__MWERKS__) || (defined(__BORLANDC__) && (__BORLANDC__ >= 0x540))
00164 #define FRIEND_FUN_TEMPLATE <>
00165 #else
00166 #define FRIEND_FUN_TEMPLATE
00167 #endif
00168 
00169 #if defined(_MSC_VER) && _MSC_VER <= 1020   // MSVC++ 4.0/4.2 does not
00170 #  define _NO_NAMESPACE                     // support "std" namespace
00171 #endif
00172 
00173 #if !defined(_NO_NAMESPACE)
00174 #if defined( _SGI_BROKEN_STL )              // For SGI C++ v.7.2.1 compiler
00175 namespace std { }
00176 #endif
00177 using namespace std;
00178 #endif
00179 
00180 #ifndef _NO_NAMESPACE
00181 namespace math {
00182 #endif
00183 
00184 #if !defined(_NO_EXCEPTION)
00185 class matrix_error : public logic_error
00186 {
00187     public:
00188         matrix_error (const string& what_arg) : logic_error( what_arg) {}
00189 };
00190 #define REPORT_ERROR(ErrormMsg)  throw matrix_error( ErrormMsg);
00191 #else
00192 inline void _matrix_error (const char* pErrMsg)
00193 {
00194     cout << pErrMsg << endl;
00195     exit(1);
00196 }
00197 #define REPORT_ERROR(ErrormMsg)  _matrix_error( ErrormMsg);
00198 #endif
00199 
00200 #if !defined(_NO_TEMPLATE)
00201 #  define MAT_TEMPLATE  template <class T>
00202 #  define matrixT  matrix<T>
00203 #else
00204 #  define MAT_TEMPLATE
00205 #  define matrixT  matrix
00206 #  ifdef MATRIX_TYPE
00207      typedef MATRIX_TYPE T;
00208 #  else
00209      typedef double T;
00210 #  endif
00211 #endif
00212 
00213 
00214 MAT_TEMPLATE
00215 class matrix
00216 {
00217 public:
00218    // Constructors
00219    matrix (const matrixT& m);
00220    matrix (size_t row = 6, size_t col = 6);
00221 
00222    // Destructor
00223    ~matrix ();
00224 
00225    // Assignment operators
00226    matrixT& operator = (const matrixT& m) _NO_THROW;
00227 
00228    // Value extraction method
00229    size_t RowNo () const { return _m->Row; }
00230    size_t ColNo () const { return _m->Col; }
00231 
00232    // Subscript operator
00233    T& operator () (size_t row, size_t col) _THROW_MATRIX_ERROR;
00234    T  operator () (size_t row, size_t col) const _THROW_MATRIX_ERROR;
00235 
00236    // Unary operators
00237    matrixT operator + () _NO_THROW { return *this; }
00238    matrixT operator - () _NO_THROW;
00239 
00240    // Combined assignment - calculation operators
00241    matrixT& operator += (const matrixT& m) _THROW_MATRIX_ERROR;
00242    matrixT& operator -= (const matrixT& m) _THROW_MATRIX_ERROR;
00243    matrixT& operator *= (const matrixT& m) _THROW_MATRIX_ERROR;
00244    matrixT& operator *= (const T& c) _NO_THROW;
00245    matrixT& operator /= (const T& c) _NO_THROW;
00246    matrixT& operator ^= (const size_t& pow) _THROW_MATRIX_ERROR;
00247 
00248    // Miscellaneous -methods
00249    void Null (const size_t& row, const size_t& col) _NO_THROW;
00250    void Null () _NO_THROW;
00251    void Unit (const size_t& row) _NO_THROW;
00252    void Unit () _NO_THROW;
00253    void SetSize (size_t row, size_t col) _NO_THROW;
00254 
00255    // Utility methods
00256    matrixT Solve (const matrixT& v) const _THROW_MATRIX_ERROR;
00257    matrixT Adj () _THROW_MATRIX_ERROR;
00258    matrixT Inv () _THROW_MATRIX_ERROR;
00259    T Det () const _THROW_MATRIX_ERROR;
00260    T Norm () _NO_THROW;
00261    T Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR;
00262    T Cond () _NO_THROW;
00263 
00264    // Type of matrices
00265    bool IsSquare () _NO_THROW { return (_m->Row == _m->Col); }
00266    bool IsSingular () _NO_THROW;
00267    bool IsDiagonal () _NO_THROW;
00268    bool IsScalar () _NO_THROW;
00269    bool IsUnit () _NO_THROW;
00270    bool IsNull () _NO_THROW;
00271    bool IsSymmetric () _NO_THROW;
00272    bool IsSkewSymmetric () _NO_THROW;
00273    bool IsUpperTriangular () _NO_THROW;
00274    bool IsLowerTriangular () _NO_THROW;
00275 
00276 private:
00277     struct base_mat
00278     {
00279         T **Val;
00280         size_t Row, Col, RowSiz, ColSiz;
00281         int Refcnt;
00282 
00283         base_mat (size_t row, size_t col, T** v)
00284         {
00285             Row = row; RowSiz = row;
00286             Col = col; ColSiz = col;
00287             Refcnt = 1;
00288 
00289             Val = new T* [row];
00290             size_t rowlen = col * sizeof(T);
00291 
00292             for (size_t i=0; i < row; i++)
00293             {
00294                 Val[i] = new T [col];
00295                 if (v) memcpy( Val[i], v[i], rowlen);
00296             }
00297         }
00298         ~base_mat ()
00299         {
00300             for (size_t i=0; i < RowSiz; i++)
00301                 delete [] Val[i];
00302             delete [] Val;
00303         }
00304     };
00305     base_mat *_m;
00306 
00307     void clone ();
00308     void realloc (size_t row, size_t col);
00309     int pivot (size_t row);
00310 };
00311 
00312 #if defined(_MSC_VER) && _MSC_VER <= 1020
00313 #  undef  _NO_THROW               // MSVC++ 4.0/4.2 does not support
00314 #  undef  _THROW_MATRIX_ERROR     // exception specification in definition
00315 #  define _NO_THROW
00316 #  define _THROW_MATRIX_ERROR
00317 #endif
00318 
00319 // constructor
00320 MAT_TEMPLATE inline
00321 matrixT::matrix (size_t row, size_t col)
00322 {
00323   _m = new base_mat( row, col, 0);
00324 }
00325 
00326 // copy constructor
00327 MAT_TEMPLATE inline
00328 matrixT::matrix (const matrixT& m)
00329 {
00330     _m = m._m;
00331     _m->Refcnt++;
00332 }
00333 
00334 // Internal copy constructor
00335 MAT_TEMPLATE inline void
00336 matrixT::clone ()
00337 {
00338     _m->Refcnt--;
00339     _m = new base_mat( _m->Row, _m->Col, _m->Val);
00340 }
00341 
00342 // destructor
00343 MAT_TEMPLATE inline
00344 matrixT::~matrix ()
00345 {
00346    if (--_m->Refcnt == 0) delete _m;
00347 }
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 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 // binary scalar division operator
00619 MAT_TEMPLATE inline matrixT
00620 operator / (const matrixT& m, const T& no) _NO_THROW
00621 {
00622     return (m * (T(1) / no));
00623 }
00624 
00625 
00626 // binary scalar division operator
00627 MAT_TEMPLATE inline matrixT
00628 operator / (const T& no, const matrixT& m) _THROW_MATRIX_ERROR
00629 {
00630     return (!m * no);
00631 }
00632 
00633 // binary matrix division operator
00634 MAT_TEMPLATE inline matrixT
00635 operator / (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
00636 {
00637     return (m1 * !m2);
00638 }
00639 
00640 // binary power operator
00641 MAT_TEMPLATE inline matrixT
00642 operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR
00643 {
00644    matrixT temp = m;
00645    temp ^= pow;
00646    return temp;
00647 }
00648 
00649 // unary transpose operator
00650 MAT_TEMPLATE inline matrixT
00651 operator ~ (const matrixT& m) _NO_THROW
00652 {
00653    matrixT temp(m.ColNo(),m.RowNo());
00654 
00655    for (size_t i=0; i < m.RowNo(); i++)
00656       for (size_t j=0; j < m.ColNo(); j++)
00657       {
00658          T x = m(i,j);
00659               temp(j,i) = x;
00660       }
00661    return temp;
00662 }
00663 
00664 // unary inversion operator
00665 MAT_TEMPLATE inline matrixT
00666 operator ! (const matrixT m) _THROW_MATRIX_ERROR
00667 {
00668    matrixT temp = m;
00669    return temp.Inv();
00670 }
00671 
00672 // inversion function
00673 MAT_TEMPLATE inline matrixT
00674 matrixT::Inv () _THROW_MATRIX_ERROR
00675 {
00676    size_t i,j,k;
00677    T a1,a2,*rowptr;
00678 
00679    if (_m->Row != _m->Col)
00680       REPORT_ERROR( "matrixT::operator!: Inversion of a non-square matrix");
00681 
00682    matrixT temp(_m->Row,_m->Col);
00683    if (_m->Refcnt > 1) clone();
00684 
00685 
00686    temp.Unit();
00687    for (k=0; k < _m->Row; k++)
00688    {
00689       int indx = pivot(k);
00690       if (indx == -1)
00691               REPORT_ERROR( "matrixT::operator!: Inversion of a singular matrix");
00692 
00693       if (indx != 0)
00694       {
00695               rowptr = temp._m->Val[k];
00696               temp._m->Val[k] = temp._m->Val[indx];
00697               temp._m->Val[indx] = rowptr;
00698       }
00699       a1 = _m->Val[k][k];
00700       for (j=0; j < _m->Row; j++)
00701       {
00702               _m->Val[k][j] /= a1;
00703               temp._m->Val[k][j] /= a1;
00704       }
00705       for (i=0; i < _m->Row; i++)
00706            if (i != k)
00707            {
00708               a2 = _m->Val[i][k];
00709               for (j=0; j < _m->Row; j++)
00710               {
00711                  _m->Val[i][j] -= a2 * _m->Val[k][j];
00712                  temp._m->Val[i][j] -= a2 * temp._m->Val[k][j];
00713               }
00714            }
00715    }
00716    return temp;
00717 }
00718 
00719 // solve simultaneous equation
00720 MAT_TEMPLATE inline matrixT
00721 matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR
00722 {
00723    size_t i,j,k;
00724    T a1;
00725 
00726    if (!(_m->Row == _m->Col && _m->Col == v._m->Row))
00727       REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!");
00728 
00729    matrixT temp(_m->Row,_m->Col+v._m->Col);
00730    for (i=0; i < _m->Row; i++)
00731    {
00732       for (j=0; j < _m->Col; j++)
00733          temp._m->Val[i][j] = _m->Val[i][j];
00734       for (k=0; k < v._m->Col; k++)
00735          temp._m->Val[i][_m->Col+k] = v._m->Val[i][k];
00736    }
00737    for (k=0; k < _m->Row; k++)
00738    {
00739       int indx = temp.pivot(k);
00740       if (indx == -1)
00741          REPORT_ERROR( "matrixT::Solve(): Singular matrix!");
00742 
00743       a1 = temp._m->Val[k][k];
00744       for (j=k; j < temp._m->Col; j++)
00745          temp._m->Val[k][j] /= a1;
00746 
00747       for (i=k+1; i < _m->Row; i++)
00748       {
00749          a1 = temp._m->Val[i][k];
00750          for (j=k; j < temp._m->Col; j++)
00751            temp._m->Val[i][j] -= a1 * temp._m->Val[k][j];
00752       }
00753    }
00754    matrixT s(v._m->Row,v._m->Col);
00755    for (k=0; k < v._m->Col; k++)
00756       for (int m=int(_m->Row)-1; m >= 0; m--)
00757       {
00758          s._m->Val[m][k] = temp._m->Val[m][_m->Col+k];
00759          for (j=m+1; j < _m->Col; j++)
00760             s._m->Val[m][k] -= temp._m->Val[m][j] * s._m->Val[j][k];
00761       }
00762    return s;
00763 }
00764 
00765 // set zero to all elements of this matrix
00766 MAT_TEMPLATE inline void
00767 matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
00768 {
00769     if (row != _m->Row || col != _m->Col)
00770         realloc( row,col);
00771 
00772     if (_m->Refcnt > 1)
00773         clone();
00774 
00775     for (size_t i=0; i < _m->Row; i++)
00776         for (size_t j=0; j < _m->Col; j++)
00777             _m->Val[i][j] = T(0);
00778     return;
00779 }
00780 
00781 // set zero to all elements of this matrix
00782 MAT_TEMPLATE inline void
00783 matrixT::Null() _NO_THROW
00784 {
00785     if (_m->Refcnt > 1) clone();
00786     for (size_t i=0; i < _m->Row; i++)
00787         for (size_t j=0; j < _m->Col; j++)
00788                 _m->Val[i][j] = T(0);
00789     return;
00790 }
00791 
00792 // set this matrix to unity
00793 MAT_TEMPLATE inline void
00794 matrixT::Unit (const size_t& row) _NO_THROW
00795 {
00796     if (row != _m->Row || row != _m->Col)
00797         realloc( row, row);
00798 
00799     if (_m->Refcnt > 1)
00800         clone();
00801 
00802     for (size_t i=0; i < _m->Row; i++)
00803         for (size_t j=0; j < _m->Col; j++)
00804             _m->Val[i][j] = i == j ? T(1) : T(0);
00805     return;
00806 }
00807 
00808 // set this matrix to unity
00809 MAT_TEMPLATE inline void
00810 matrixT::Unit () _NO_THROW
00811 {
00812     if (_m->Refcnt > 1) clone();
00813     size_t row = min(_m->Row,_m->Col);
00814     _m->Row = _m->Col = row;
00815 
00816     for (size_t i=0; i < _m->Row; i++)
00817         for (size_t j=0; j < _m->Col; j++)
00818             _m->Val[i][j] = i == j ? T(1) : T(0);
00819     return;
00820 }
00821 
00822 // private partial pivoting method
00823 MAT_TEMPLATE inline int
00824 matrixT::pivot (size_t row)
00825 {
00826   int k = int(row);
00827   double amax,temp;
00828 
00829   amax = -1;
00830   for (size_t i=row; i < _m->Row; i++)
00831     if ( (temp = abs( _m->Val[i][row])) > amax && temp != 0.0)
00832      {
00833        amax = temp;
00834        k = i;
00835      }
00836   if (_m->Val[k][row] == T(0))
00837      return -1;
00838   if (k != int(row))
00839   {
00840      T* rowptr = _m->Val[k];
00841      _m->Val[k] = _m->Val[row];
00842      _m->Val[row] = rowptr;
00843      return k;
00844   }
00845   return 0;
00846 }
00847 
00848 // calculate the determinant of a matrix
00849 MAT_TEMPLATE inline T
00850 matrixT::Det () const _THROW_MATRIX_ERROR
00851 {
00852    size_t i,j,k;
00853    T piv,detVal = T(1);
00854 
00855    if (_m->Row != _m->Col)
00856       REPORT_ERROR( "matrixT::Det(): Determinant a non-square matrix!");
00857 
00858    matrixT temp(*this);
00859    if (temp._m->Refcnt > 1) temp.clone();
00860 
00861    for (k=0; k < _m->Row; k++)
00862    {
00863       int indx = temp.pivot(k);
00864       if (indx == -1)
00865          return 0;
00866       if (indx != 0)
00867          detVal = - detVal;
00868       detVal = detVal * temp._m->Val[k][k];
00869       for (i=k+1; i < _m->Row; i++)
00870       {
00871          piv = temp._m->Val[i][k] / temp._m->Val[k][k];
00872          for (j=k+1; j < _m->Row; j++)
00873             temp._m->Val[i][j] -= piv * temp._m->Val[k][j];
00874       }
00875    }
00876    return detVal;
00877 }
00878 
00879 // calculate the norm of a matrix
00880 MAT_TEMPLATE inline T
00881 matrixT::Norm () _NO_THROW
00882 {
00883    T retVal = T(0);
00884 
00885    for (size_t i=0; i < _m->Row; i++)
00886       for (size_t j=0; j < _m->Col; j++)
00887          retVal += _m->Val[i][j] * _m->Val[i][j];
00888    retVal = sqrt( retVal);
00889 
00890    return retVal;
00891 }
00892 
00893 // calculate the condition number of a matrix
00894 MAT_TEMPLATE inline T
00895 matrixT::Cond () _NO_THROW
00896 {
00897    matrixT inv = ! (*this);
00898    return (Norm() * inv.Norm());
00899 }
00900 
00901 // calculate the cofactor of a matrix for a given element
00902 MAT_TEMPLATE inline T
00903 matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR
00904 {
00905    size_t i,i1,j,j1;
00906 
00907    if (_m->Row != _m->Col)
00908       REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!");
00909 
00910    if (row > _m->Row || col > _m->Col)
00911       REPORT_ERROR( "matrixT::Cofact(): Index out of range!");
00912 
00913    matrixT temp (_m->Row-1,_m->Col-1);
00914 
00915    for (i=i1=0; i < _m->Row; i++)
00916    {
00917       if (i == row)
00918         continue;
00919       for (j=j1=0; j < _m->Col; j++)
00920       {
00921          if (j == col)
00922             continue;
00923          temp._m->Val[i1][j1] = _m->Val[i][j];
00924          j1++;
00925       }
00926       i1++;
00927    }
00928    T  cof = temp.Det();
00929    if ((row+col)%2 == 1)
00930       cof = -cof;
00931 
00932    return cof;
00933 }
00934 
00935 
00936 // calculate adjoin of a matrix
00937 MAT_TEMPLATE inline matrixT
00938 matrixT::Adj () _THROW_MATRIX_ERROR
00939 {
00940    if (_m->Row != _m->Col)
00941       REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix.");
00942 
00943    matrixT temp(_m->Row,_m->Col);
00944 
00945    for (size_t i=0; i < _m->Row; i++)
00946       for (size_t j=0; j < _m->Col; j++)
00947           temp._m->Val[j][i] = Cofact(i,j);
00948    return temp;
00949 }
00950 
00951 // Determine if the matrix is singular
00952 MAT_TEMPLATE inline bool
00953 matrixT::IsSingular () _NO_THROW
00954 {
00955    if (_m->Row != _m->Col)
00956       return false;
00957    return (Det() == T(0));
00958 }
00959 
00960 // Determine if the matrix is diagonal
00961 MAT_TEMPLATE inline bool
00962 matrixT::IsDiagonal () _NO_THROW
00963 {
00964    if (_m->Row != _m->Col)
00965       return false;
00966    for (size_t i=0; i < _m->Row; i++)
00967      for (size_t j=0; j < _m->Col; j++)
00968         if (i != j && _m->Val[i][j] != T(0))
00969           return false;
00970    return true;
00971 }
00972 
00973 // Determine if the matrix is scalar
00974 MAT_TEMPLATE inline bool
00975 matrixT::IsScalar () _NO_THROW
00976 {
00977    if (!IsDiagonal())
00978      return false;
00979    T v = _m->Val[0][0];
00980    for (size_t i=1; i < _m->Row; i++)
00981      if (_m->Val[i][i] != v)
00982         return false;
00983    return true;
00984 }
00985 
00986 // Determine if the matrix is a unit matrix
00987 MAT_TEMPLATE inline bool
00988 matrixT::IsUnit () _NO_THROW
00989 {
00990    if (IsScalar() && _m->Val[0][0] == T(1))
00991      return true;
00992    return false;
00993 }
00994 
00995 // Determine if this is a null matrix
00996 MAT_TEMPLATE inline bool
00997 matrixT::IsNull () _NO_THROW
00998 {
00999    for (size_t i=0; i < _m->Row; i++)
01000       for (size_t j=0; j < _m->Col; j++)
01001          if (_m->Val[i][j] != T(0))
01002             return false;
01003    return true;
01004 }
01005 
01006 // Determine if the matrix is symmetric
01007 MAT_TEMPLATE inline bool
01008 matrixT::IsSymmetric () _NO_THROW
01009 {
01010    if (_m->Row != _m->Col)
01011       return false;
01012    for (size_t i=0; i < _m->Row; i++)
01013       for (size_t j=0; j < _m->Col; j++)
01014          if (_m->Val[i][j] != _m->Val[j][i])
01015             return false;
01016    return true;
01017 }
01018 
01019 // Determine if the matrix is skew-symmetric
01020 MAT_TEMPLATE inline bool
01021 matrixT::IsSkewSymmetric () _NO_THROW
01022 {
01023    if (_m->Row != _m->Col)
01024       return false;
01025    for (size_t i=0; i < _m->Row; i++)
01026       for (size_t j=0; j < _m->Col; j++)
01027          if (_m->Val[i][j] != -_m->Val[j][i])
01028             return false;
01029    return true;
01030 }
01031 
01032 // Determine if the matrix is upper triangular
01033 MAT_TEMPLATE inline bool
01034 matrixT::IsUpperTriangular () _NO_THROW
01035 {
01036    if (_m->Row != _m->Col)
01037       return false;
01038    for (size_t i=1; i < _m->Row; i++)
01039       for (size_t j=0; j < i-1; j++)
01040          if (_m->Val[i][j] != T(0))
01041             return false;
01042    return true;
01043 }
01044 
01045 // Determine if the matrix is lower triangular
01046 MAT_TEMPLATE inline bool
01047 matrixT::IsLowerTriangular () _NO_THROW
01048 {
01049    if (_m->Row != _m->Col)
01050       return false;
01051 
01052    for (size_t j=1; j < _m->Col; j++)
01053       for (size_t i=0; i < j-1; i++)
01054          if (_m->Val[i][j] != T(0))
01055             return false;
01056 
01057    return true;
01058 }
01059 
01060 #ifndef _NO_NAMESPACE
01061 }
01062 #endif
01063 
01064 #endif //__STD_MATRIX_H


ros_rt_wmp_sniffer
Author(s): Danilo Tardioli, dantard@unizar.es
autogenerated on Mon Oct 6 2014 08:27:57