00001 #ifndef MATRIX_H
00002 #define MATRIX_H
00003
00004 #include <cmath>
00005 #include <cstdio>
00006 #include <cstdlib>
00007 #include <iostream>
00008 #include <fstream>
00009 #include <cstring>
00010
00011
00012
00013 template <class T>
00014 class Matrix
00015 {
00016 public:
00017
00018
00019 Matrix (size_t row = 0, size_t col = 0)
00020 {
00021 SetSize(row,col);
00022 }
00023
00024
00025 virtual ~Matrix()
00026 {
00027 FreeMatrix();
00028 }
00029
00030
00031 Matrix (const Matrix<T>& m)
00032 {
00033 SetSize(m.GetRows(),m.GetCols());
00034 for (size_t i=0;i<rows;i++)
00035 {
00036 for (size_t j=0;j<cols;j++)
00037 {
00038 matrix[i][j] = m[i][j];
00039 }
00040 }
00041 }
00042
00043
00044 void SetSize (size_t row, size_t col)
00045 {
00046 rows = row;
00047 cols = col;
00048 if (rows > 0 && cols > 0)
00049 {
00050 AllocateMatrix();
00051 }
00052 else
00053 {
00054 matrix = NULL;
00055 }
00056 }
00057
00058
00059 void ReSize (size_t row, size_t col)
00060 {
00061 if (row == rows && col == cols)
00062 {
00063 Zeros();
00064 return;
00065 }
00066 FreeMatrix();
00067 SetSize(row,col);
00068 }
00069
00070
00071 Matrix<T>& operator = (const Matrix<T>& m)
00072 {
00073 ReSize(m.GetRows(),m.GetCols());
00074 for (size_t i=0;i<rows;i++)
00075 {
00076 for (size_t j=0;j<cols;j++)
00077 {
00078 matrix[i][j] = m[i][j];
00079 }
00080 }
00081 return *this;
00082 }
00083
00084
00085 size_t GetRows () const
00086 {
00087 return rows;
00088 }
00089
00090
00091 size_t GetCols () const
00092 {
00093 return cols;
00094 }
00095
00096
00097 T*& operator[](int row)
00098 {
00099 return matrix[row];
00100 }
00101
00102
00103 T* operator[](int row) const
00104 {
00105 return matrix[row];
00106 }
00107
00108 void set(size_t i,size_t j,T value)
00109 {
00110 matrix[i][j] = value;
00111 }
00112
00113 T get(size_t i,size_t j)
00114 {
00115 return matrix[i][j];
00116 }
00117
00118
00119 T& operator () (size_t row, size_t col);
00120 T operator () (size_t row, size_t col) const;
00121
00122
00123 Matrix<T> operator + ();
00124 Matrix<T> operator - ();
00125
00126 Matrix<T>& operator += (const T& c);
00127 Matrix<T>& operator -= (const T& c);
00128 Matrix<T>& operator *= (const T& c);
00129 Matrix<T>& operator /= (const T& c);
00130
00131 Matrix<T>& operator += (const Matrix<T>& m);
00132 Matrix<T>& operator -= (const Matrix<T>& m);
00133 Matrix<T>& operator *= (const Matrix<T>& m);
00134 Matrix<T>& operator /= (const Matrix<T>& m);
00135
00136 Matrix<T>& operator ^= (const int pow);
00137
00138 void SwapRows(size_t n1, size_t n2);
00139 void SwapCols(size_t m1, size_t m2);
00140
00141
00142 void Zeros ();
00143 void Ones ();
00144 void Identity ();
00145
00146
00147 Matrix<T> Inv ();
00148 T Det () const;
00149 T Norm () const;
00150 T TotalSum () const;
00151
00152
00153 bool IsSquare () const;
00154 bool IsSingular () const;
00155 bool IsDiagonal () const;
00156 bool IsScalar () const;
00157 bool IsIdentity () const;
00158 bool IsEmpty () const;
00159 bool IsNull () const;
00160 bool IsSymmetric () const;
00161 bool IsUpperTriangular () const;
00162 bool IsLowerTriangular () const;
00163
00164 void FreeMatrix()
00165 {
00166 if(matrix != NULL)
00167 {
00168
00169 for(size_t i=0; i<rows ;i++)
00170 {
00171 delete [] matrix[i];
00172 }
00173 delete [] matrix;
00174 matrix = NULL;
00175 }
00176 }
00177
00178 T** GetMatrix() {return matrix;}
00179 void transform(Matrix<T>& trans)
00180 {
00181 Matrix<T> temp(4,1);
00182 for(unsigned int i=0; i<rows; ++i)
00183 {
00184 temp[0][0] = matrix[i][0];
00185 temp[1][0] = matrix[i][1];
00186 temp[2][0] = matrix[i][2];
00187 temp[3][0] = 1;
00188 temp = trans * temp;
00189 matrix[i][0] = temp[0][0];
00190 matrix[i][1] = temp[1][0];
00191 matrix[i][2] = temp[2][0];
00192 }
00193 }
00194
00195 protected:
00196
00197 T** matrix;
00198 size_t rows;
00199 size_t cols;
00200
00201 void AllocateMatrix()
00202 {
00203 matrix = new T* [rows];
00204 for(size_t i=0;i<rows;i++)
00205 {
00206 matrix[i] = new T[cols];
00207 memset(matrix[i],0,sizeof(T)*cols/sizeof(char));
00208 }
00209 }
00210 };
00211
00212
00213 template <class T>
00214 class Row : public Matrix<T>
00215 {
00216 public:
00217
00218
00219 Row (size_t col = 0) : Matrix<T>(1,col)
00220 {
00221 }
00222
00223
00224 virtual ~Row()
00225 {
00226 }
00227
00228
00229 void SetRowSize(size_t col)
00230 {
00231 this->SetSize(1,col);
00232 }
00233
00234
00235 size_t GetSize() const
00236 {
00237 return this->cols;
00238 }
00239
00240 T& operator[](int idx)
00241 {
00242 return this->matrix[0][idx];
00243 }
00244
00245 T operator[](int idx) const
00246 {
00247 return this->matrix[0][idx];
00248 }
00249
00250 T& operator () (size_t idx)
00251 {
00252 return this->matrix[0][idx-1];
00253 }
00254
00255 T operator () (size_t idx) const
00256 {
00257 return this->matrix[0][idx-1];
00258 }
00259
00260
00261 Row<T>&
00262 operator = (const Row<T> &r)
00263 {
00264 SetRowSize(r.GetSize());
00265 for (size_t i=0;i<r.GetSize();i++)
00266 {
00267 this->matrix[0][i] = r[i];
00268 }
00269 return *this;
00270 }
00271 };
00272
00273
00274 template <class T>
00275 class Column : public Matrix<T>
00276 {
00277 public:
00278
00279
00280 Column (size_t n = 0) : Matrix<T>(n,1)
00281 {
00282 }
00283
00284
00285 virtual ~Column()
00286 {
00287 }
00288
00289
00290 void SetColumnSize(size_t n)
00291 {
00292 this->SetSize(n,1);
00293 }
00294
00295
00296 size_t GetSize() const
00297 {
00298 return this->rows;
00299 }
00300
00301 T& operator[](int idx)
00302 {
00303 return this->matrix[idx][0];
00304 }
00305
00306 T operator[](int idx) const
00307 {
00308 return this->matrix[idx][0];
00309 }
00310
00311 T& operator () (size_t idx)
00312 {
00313 return this->matrix[idx-1][0];
00314 }
00315
00316 T operator () (size_t idx) const
00317 {
00318 return this->matrix[idx-1][0];
00319 }
00320
00321
00322 Column<T>&
00323 operator = (const Column<T> &c)
00324 {
00325 SetColumnSize(c.GetSize());
00326 for (size_t i=0;i<c.GetSize();i++)
00327 {
00328 this->matrix[i][0] = c[i];
00329 }
00330 return *this;
00331 }
00332 };
00333
00334
00335 inline void MatrixError (const char* msg)
00336 {
00337 std::cout << msg << std::endl;
00338 exit(1);
00339 }
00340
00341 template <class T> inline void
00342 swap(T *a, T *b)
00343 {
00344 T temp = *b;
00345 *b = *a;
00346 *a = temp;
00347 }
00348
00349 template <class T> inline void
00350 Matrix<T>::SwapRows(size_t n1, size_t n2)
00351 {
00352 for (size_t i=0;i<cols;i++)
00353 {
00354 T temp = matrix[n2-1][i];
00355 matrix[n2-1][i] = matrix[n1-1][i];
00356 matrix[n1-1][i] = temp;
00357 }
00358 }
00359
00360 template <class T> inline void
00361 Matrix<T>::SwapCols(size_t m1, size_t m2)
00362 {
00363 for (size_t i=0;i<rows;i++)
00364 {
00365 T temp = matrix[i][m2-1];
00366 matrix[i][m2-1] = matrix[i][m1-1];
00367 matrix[i][m1-1] = temp;
00368 }
00369 }
00370
00371
00372 template <class T> inline T&
00373 Matrix<T>::operator () (size_t row, size_t col)
00374 {
00375 if (row == 0 || row >= rows + 1 || col == 0 || col >= cols + 1)
00376 {
00377 MatrixError( "Matrix<T>::operator (): Index out of range!");
00378 }
00379 return matrix[row-1][col-1];
00380 }
00381
00382
00383 template <class T> inline T
00384 Matrix<T>::operator () (size_t row, size_t col) const
00385 {
00386 if (row == 0 || row >= rows + 1 || col == 0 || col >= cols + 1)
00387 {
00388 MatrixError( "Matrix<T>::operator (): Index out of range!");
00389 }
00390 return matrix[row-1][col-1];
00391 }
00392
00393
00394 template <class T> inline Matrix<T>&
00395 Matrix<T>::operator += (const T& c)
00396 {
00397 for (size_t i=0; i < rows; i++)
00398 {
00399 for (size_t j=0; j < cols; j++)
00400 {
00401 matrix[i][j] += c;
00402 }
00403 }
00404 return *this;
00405 }
00406
00407
00408 template <class T> inline Matrix<T>&
00409 Matrix<T>::operator += (const Matrix<T>& m)
00410 {
00411 if (rows != m.GetRows() || cols != m.GetCols())
00412 {
00413 MatrixError( "Matrix<T>::operator += : Inconsistent Matrix sizes in addition!");
00414 }
00415
00416 for (size_t i=0; i < m.GetRows(); i++)
00417 {
00418 for (size_t j=0; j < m.GetCols(); j++)
00419 {
00420 matrix[i][j] += m[i][j];
00421 }
00422 }
00423 return *this;
00424 }
00425
00426
00427 template <class T> inline Matrix<T>&
00428 Matrix<T>::operator -= (const T& c)
00429 {
00430 for (size_t i=0; i < rows; i++)
00431 {
00432 for (size_t j=0; j < cols; j++)
00433 {
00434 matrix[i][j] -= c;
00435 }
00436 }
00437 return *this;
00438 }
00439
00440
00441 template <class T> inline Matrix<T>&
00442 Matrix<T>::operator -= (const Matrix<T>& m)
00443 {
00444 if (rows != m.GetRows() || cols != m.GetCols())
00445 {
00446 MatrixError( "Matrix<T>::operator -= : Inconsistent Matrix sizes in subtraction!");
00447 }
00448
00449 for (size_t i=0; i < m.GetRows(); i++)
00450 {
00451 for (size_t j=0; j < m.GetCols(); j++)
00452 {
00453 matrix[i][j] -= m[i][j];
00454 }
00455 }
00456 return *this;
00457 }
00458
00459
00460 template <class T> inline Matrix<T>&
00461 Matrix<T>::operator *= (const T& c)
00462 {
00463 for (size_t i=0; i < rows; i++)
00464 {
00465 for (size_t j=0; j < cols; j++)
00466 {
00467 matrix[i][j] *= c;
00468 }
00469 }
00470 return *this;
00471 }
00472
00473
00474 template <class T> inline Matrix<T>&
00475 Matrix<T>::operator *= (const Matrix<T>& m)
00476 {
00477 if (cols != m.GetRows())
00478 {
00479 MatrixError( "Matrix<T>::operator *= : Inconsistent Matrix sizes in multiplication!");
00480 }
00481
00482 Matrix<T> temp(rows,m.GetCols());
00483
00484 for (size_t i=0; i < rows; i++)
00485 {
00486 for (size_t j=0; j < m.GetCols(); j++)
00487 {
00488 temp[i][j] = T(0);
00489 for (size_t k=0; k < cols; k++)
00490 {
00491 temp[i][j] += matrix[i][k] * m[k][j];
00492 }
00493 }
00494 }
00495 *this = temp;
00496
00497 return *this;
00498 }
00499
00500
00501 template <class T> inline Matrix<T>&
00502 Matrix<T>::operator /= (const T& c)
00503 {
00504 for (size_t i=0; i < rows; i++)
00505 {
00506 for (size_t j=0; j < cols; j++)
00507 {
00508 matrix[i][j] /= c;
00509 }
00510 }
00511 return *this;
00512 }
00513
00514
00515 template <class T> inline Matrix<T>&
00516 Matrix<T>::operator ^= (const int pow)
00517 {
00518 Matrix<T> temp(*this);
00519 int i = 1;
00520 while (i < pow)
00521 {
00522 *this *= temp;
00523 i++;
00524 }
00525 return *this;
00526 }
00527
00528
00529 template <class T> inline Matrix<T>
00530 operator ~ (const Matrix<T>& m)
00531 {
00532 Matrix<T> temp(m.GetCols(),m.GetRows());
00533 for (size_t i=0; i < m.GetRows(); i++)
00534 {
00535 for (size_t j=0; j < m.GetCols(); j++)
00536 {
00537 temp[j][i] = m[i][j];
00538 }
00539 }
00540 return temp;
00541 }
00542
00543
00544 template <class T> inline Column<T>
00545 operator ~ (const Row<T>& r)
00546 {
00547 Column<T> temp(r.GetSize());
00548 for (size_t i=0; i < r.GetSize(); i++)
00549 {
00550 temp[i] = r[i];
00551 }
00552 return temp;
00553 }
00554
00555
00556 template <class T> inline Row<T>
00557 operator ~ (const Column<T>& c)
00558 {
00559 Row<T> temp(c.GetSize());
00560 for (size_t i=0; i < c.GetSize(); i++)
00561 {
00562 temp[i] = c[i];
00563 }
00564 return temp;
00565 }
00566
00567
00568 template <class T> inline Matrix<T>
00569 operator ! (const Matrix<T>& m)
00570 {
00571 Matrix<T> temp(m);
00572 temp.Inv();
00573 return temp;
00574 }
00575
00576
00577 template <class T> inline Matrix<T>
00578 Matrix<T>::operator + ()
00579 {
00580 return *this;
00581 }
00582
00583
00584 template <class T> inline Matrix<T>
00585 Matrix<T>::operator - ()
00586 {
00587 Matrix<T> temp(rows,cols);
00588 for (size_t i=0; i < rows; i++)
00589 {
00590 for (size_t j=0; j < cols; j++)
00591 {
00592 temp[i][j] = - matrix[i][j];
00593 }
00594 }
00595 return temp;
00596 }
00597
00598
00599 template <class T> inline Matrix<T>
00600 operator + (const Matrix<T>& m1, const Matrix<T>& m2)
00601 {
00602 if (m1.GetRows() != m2.GetRows() || m1.GetCols() != m2.GetCols())
00603 {
00604 MatrixError( "Matrix<T>::operator + : Inconsistent Matrix sizes in addition!");
00605 }
00606 Matrix<T> temp(m1);
00607 temp += m2;
00608 return temp;
00609 }
00610
00611
00612 template <class T> inline Matrix<T>
00613 operator + (const Matrix<T>& m, const T& c)
00614 {
00615 Matrix<T> temp(m);
00616 temp += c;
00617 return temp;
00618 }
00619
00620
00621 template <class T> inline Matrix<T>
00622 operator + (const T& c, const Matrix<T>& m)
00623 {
00624 return m + c;
00625 }
00626
00627
00628 template <class T> inline Matrix<T>
00629 operator - (const Matrix<T>& m1, const Matrix<T>& m2)
00630 {
00631 if (m1.GetRows() != m2.GetRows() || m1.GetCols() != m2.GetCols())
00632 {
00633 MatrixError( "Matrix<T>::operator - : Inconsistent Matrix sizes in subtraction!");
00634 }
00635
00636 Matrix<T> temp(m1);
00637 temp -= m2;
00638 return temp;
00639 }
00640
00641
00642 template <class T> inline Matrix<T>
00643 operator - (const Matrix<T>& m, const T& c)
00644 {
00645 Matrix<T> temp(m);
00646 temp -= c;
00647 return temp;
00648 }
00649
00650
00651 template <class T> inline Matrix<T>
00652 operator - (const T& c, const Matrix<T>& m)
00653 {
00654 Matrix<T> temp(m);
00655 for (size_t i=0; i < m.GetRows(); i++)
00656 {
00657 for (size_t j=0; j < m.GetCols(); j++)
00658 {
00659 temp[i][j] = c - m[i][j];
00660 }
00661 }
00662 return temp;
00663 }
00664
00665
00666 template <class T> inline Matrix<T>
00667 operator * (const Matrix<T>& m, const T& c)
00668 {
00669 Matrix<T> temp(m);
00670 temp *= c;
00671 return temp;
00672 }
00673
00674
00675 template <class T> inline Matrix<T>
00676 operator * (const T& c, const Matrix<T>& m)
00677 {
00678 return m * c;
00679 }
00680
00681
00682 template <class T> inline Matrix<T>
00683 operator * (const Matrix<T>& m1, const Matrix<T>& m2)
00684 {
00685 if (m1.GetCols() != m2.GetRows())
00686 {
00687 MatrixError( "operator * : Inconsistent Matrix sizes in multiplication!");
00688 }
00689 Matrix<T> temp(m1);
00690 temp *= m2;
00691 return temp;
00692 }
00693
00694
00695 template <class T> inline Matrix<T>
00696 operator / (const Matrix<T>& m, const T& c)
00697 {
00698 Matrix<T> temp(m);
00699 temp /= c;
00700 return temp;
00701 }
00702
00703
00704 template <class T> inline Matrix<T>
00705 operator / (const T& c, const Matrix<T>& m)
00706 {
00707 return (!m * c);
00708 }
00709
00710
00711 template <class T> inline Matrix<T>
00712 operator / (const Matrix<T>& m1, const Matrix<T>& m2)
00713 {
00714 return (m1 * !m2);
00715 }
00716
00717
00718 template <class T> inline Matrix<T>
00719 operator ^ (const Matrix<T>& m, const size_t& pow)
00720 {
00721 Matrix<T> temp(m);
00722 temp ^= pow;
00723 return temp;
00724 }
00725
00726
00727 template <class T> inline std::ostream&
00728 operator << (std::ostream& out, Matrix<T>& m)
00729 {
00730 for (size_t i=0; i < m.GetRows(); i++)
00731 {
00732 out << m[i][0];
00733 for (size_t j=1; j < m.GetCols(); j++)
00734 {
00735 out << "\t" << m[i][j] ;
00736 }
00737 out << std::endl;
00738 }
00739 return out;
00740 }
00741
00742
00743 template <class T> inline bool
00744 operator == (const Matrix<T>& m1, const Matrix<T>& m2)
00745 {
00746 if (m1.GetRows() != m2.GetRows() || m1.GetCols() != m2.GetCols())
00747 {
00748 return false;
00749 }
00750 for (size_t i=0; i < m1.GetRows(); i++)
00751 {
00752 for (size_t j=0; j < m1.GetCols(); j++)
00753 {
00754 if (m1[i][j] != m2[i][j])
00755 {
00756 return false;
00757 }
00758 }
00759 }
00760 return true;
00761 }
00762
00763
00764 template <class T> inline bool
00765 operator != (const Matrix<T>& m1, const Matrix<T>& m2)
00766 {
00767 return !(m1 == m2);
00768 }
00769
00770
00771 template <class T> inline T
00772 Matrix<T>::Norm () const
00773 {
00774 T retval = T(0);
00775 for (size_t i=0; i < rows; i++)
00776 {
00777 for (size_t j=0; j < cols; j++)
00778 {
00779 retval += matrix[i][j] * matrix[i][j];
00780 }
00781 }
00782 retval = sqrt(retval);
00783 return retval;
00784 }
00785
00786
00787 template <class T> inline bool
00788 Matrix<T>::IsSquare () const
00789 {
00790 return (rows == cols);
00791 }
00792
00793
00794 template <class T> inline bool
00795 Matrix<T>::IsSingular () const
00796 {
00797 if (rows != cols)
00798 {
00799 return false;
00800 }
00801 return (Det() == T(0));
00802 }
00803
00804
00805 template <class T> inline bool
00806 Matrix<T>::IsDiagonal () const
00807 {
00808 if (rows != cols)
00809 {
00810 return false;
00811 }
00812 for (size_t i=0; i < rows; i++)
00813 {
00814 for (size_t j=0; j < cols; j++)
00815 {
00816 if (i != j && matrix[i][j] != T(0))
00817 {
00818 return false;
00819 }
00820 }
00821 }
00822 return true;
00823 }
00824
00825
00826 template <class T> inline bool
00827 Matrix<T>::IsScalar () const
00828 {
00829 if (!IsDiagonal())
00830 {
00831 return false;
00832 }
00833 T val = matrix[0][0];
00834 for (size_t i=1; i < rows; i++)
00835 {
00836 if (matrix[i][i] != val)
00837 {
00838 return false;
00839 }
00840 }
00841 return true;
00842 }
00843
00844
00845 template <class T> inline bool
00846 Matrix<T>::IsIdentity () const
00847 {
00848 return (IsScalar() && matrix[0][0] == T(1));
00849 }
00850
00851
00852 template <class T> inline bool
00853 Matrix<T>::IsEmpty () const
00854 {
00855 return ( (rows == 0) || (cols == 0) );
00856 }
00857
00858
00859 template <class T> inline bool
00860 Matrix<T>::IsNull () const
00861 {
00862 for (size_t i=0; i < rows; i++)
00863 {
00864 for (size_t j=0; j < cols; j++)
00865 {
00866 if (matrix[i][j] != T(0))
00867 {
00868 return false;
00869 }
00870 }
00871 }
00872 return true;
00873 }
00874
00875
00876 template <class T> inline bool
00877 Matrix<T>::IsSymmetric () const
00878 {
00879 if (rows != cols)
00880 {
00881 return false;
00882 }
00883 for (size_t i=0; i < rows; i++)
00884 {
00885 for (size_t j=0; j < cols; j++)
00886 {
00887 if (abs(matrix[i][j] - matrix[j][i]) > 0.000001)
00888 {
00889 return false;
00890 }
00891 }
00892 }
00893 return true;
00894 }
00895
00896
00897 template <class T> inline bool
00898 Matrix<T>::IsUpperTriangular () const
00899 {
00900 if (rows != cols)
00901 {
00902 return false;
00903 }
00904 for (size_t i=1; i < rows; i++)
00905 {
00906 for (size_t j=0; j < i; j++)
00907 {
00908 if (matrix[i][j] != T(0))
00909 {
00910 return false;
00911 }
00912 }
00913 }
00914 return true;
00915 }
00916
00917
00918 template <class T> inline bool
00919 Matrix<T>::IsLowerTriangular () const
00920 {
00921 if (rows != cols)
00922 {
00923 return false;
00924 }
00925 for (size_t j=1; j < cols; j++)
00926 {
00927 for (size_t i=0; i < j; i++)
00928 {
00929 if (matrix[i][j] != T(0))
00930 {
00931 return false;
00932 }
00933 }
00934 }
00935 return true;
00936 }
00937
00938
00939 template <class T> inline Row<T>
00940 Sum (Matrix<T> mat)
00941 {
00942 size_t cols = mat.GetCols();
00943 size_t rows = mat.GetRows();
00944 Row<T> row(cols);
00945 for (size_t j=0; j < cols; j++)
00946 {
00947 T sum = T(0);
00948 for (size_t i=0; i < rows; i++)
00949 {
00950 sum += mat[i][j];
00951 }
00952 row[j] = sum;
00953 }
00954 return row;
00955 }
00956
00957
00958 template <class T> inline Row<T>
00959 Mean (Matrix<T> mat)
00960 {
00961 size_t cols = mat.GetCols();
00962 size_t rows = mat.GetRows();
00963 Row<T> row(cols);
00964 for (size_t j=0; j < cols; j++)
00965 {
00966 T sum = T(0);
00967 for (size_t i=0; i < rows; i++)
00968 {
00969 sum += mat[i][j];
00970 }
00971 row[j] = sum / rows;
00972 }
00973 return row;
00974 }
00975
00976
00977 template <class T> inline T
00978 Matrix<T>::TotalSum () const
00979 {
00980 T sum = T(0);
00981 for (size_t i=0; i < rows; i++)
00982 {
00983 for (size_t j=0; j < cols; j++)
00984 {
00985 sum += matrix[i][j];
00986 }
00987 }
00988 return sum;
00989 }
00990
00991
00992 template <class T> inline T
00993 TotalSum (Matrix<T> mat)
00994 {
00995 return mat.TotalSum();
00996 }
00997
00998
00999 template <class T> inline Row<T>
01000 CumSum (Matrix<T> mat)
01001 {
01002 size_t cols = mat->GetCols();
01003 size_t rows = mat->GetRows();
01004 Row<T> row(cols);
01005 T sum = T(0);
01006 for (size_t j=0; j < cols; j++)
01007 {
01008 for (size_t i=0; i < rows; i++)
01009 {
01010 sum += mat[i][j];
01011 }
01012 row[j] = sum;
01013 }
01014 return row;
01015 }
01016
01017
01018 template <class T> inline void
01019 Matrix<T>::Zeros ()
01020 {
01021 for (size_t i=0; i < rows; i++)
01022 {
01023 for (size_t j=0; j < cols; j++)
01024 {
01025 matrix[i][j] = T(0);
01026 }
01027 }
01028 }
01029
01030
01031 template <class T> inline void
01032 Matrix<T>::Ones ()
01033 {
01034 for (size_t i=0; i < rows; i++)
01035 {
01036 for (size_t j=0; j < cols; j++)
01037 {
01038 matrix[i][j] = T(1);
01039 }
01040 }
01041 }
01042
01043
01044 template <class T> inline void
01045 Matrix<T>::Identity ()
01046 {
01047 Zeros();
01048 for (size_t i=0; i < rows; i++)
01049 {
01050 matrix[i][i] = T(1);
01051 }
01052 return;
01053 }
01054
01055
01056 template <class T> inline Matrix<T>
01057 Matrix<T>::Inv ()
01058 {
01059 if (rows != cols)
01060 {
01061 MatrixError( "Matrix<T>::Inv(): Inversion of a non-square Matrix");
01062 }
01063
01064 size_t i,j,k;
01065 size_t n = rows;
01066
01067 Matrix<T> temp(*this);
01068
01069
01070 Matrix<T> inverse(n,n);
01071 inverse.Identity();
01072
01073 for (k=0;k<n;k++)
01074 {
01075
01076 size_t p = k;
01077 for (i = k+1; i < n; i++)
01078 {
01079 if (fabs(temp[i][k]) > fabs(temp[p][k]))
01080 {
01081 p = i;
01082 }
01083 }
01084
01085
01086 if (p != k)
01087 {
01088 temp.SwapRows(k+1,p+1);
01089 inverse.SwapRows(k+1,p+1);
01090 }
01091
01092 T x = temp[k][k];
01093 if (x == T(0))
01094 {
01095 MatrixError( "Matrix<T>::Inv(): Matrix is singular");
01096 }
01097 for (j=0;j<n;j++)
01098 {
01099 if (j>k)
01100 {
01101 temp[k][j] /= x;
01102 }
01103 inverse[k][j] /= x;
01104 }
01105 for (i=0;i<n;i++)
01106 {
01107
01108 if (i == k)
01109 {
01110 continue;
01111 }
01112 x = temp[i][k];
01113 for (j=0;j<n;j++)
01114 {
01115 if (j>k)
01116 {
01117 temp[i][j] -= temp[k][j] * x;
01118 }
01119 inverse[i][j] -= inverse[k][j] * x;
01120 }
01121 }
01122 }
01123
01124 *this = inverse;
01125
01126 return *this;
01127 }
01128
01129
01130 template <class T> inline T
01131 Matrix<T>::Det () const
01132 {
01133 if (rows != cols)
01134 {
01135 MatrixError( "Matrix<T>::Det(): Matrix must be square");
01136 }
01137
01138 size_t i,j,k;
01139 size_t n = rows;
01140 int sign = 1;
01141
01142 Matrix<T> temp(*this);
01143
01144 for (k=0;k<n;k++)
01145 {
01146
01147 size_t p = k;
01148 for (i = k+1; i < n; i++)
01149 {
01150 if (fabs(temp[i][k]) > fabs(temp[p][k]))
01151 {
01152 p = i;
01153 }
01154 }
01155
01156 if (p != k)
01157 {
01158 temp.SwapRows(k+1,p+1);
01159 sign = -sign;
01160 }
01161
01162 if (temp[k][k] == T(0))
01163 {
01164
01165 return T(0);
01166 }
01167 for (j=k+1;j<n;j++)
01168 {
01169 temp[k][j] /= temp[k][k];
01170 }
01171 for (i=0;i<n;i++)
01172 {
01173
01174 if (i == k)
01175 {
01176 continue;
01177 }
01178 for (j=k+1;j<n;j++)
01179 {
01180 temp[i][j] -= temp[k][j] * temp[i][k];
01181 }
01182 }
01183 }
01184
01185
01186 T det = T(1) * sign;
01187 for (i = 0; i < n; i++)
01188 {
01189 det *= temp[i][i];
01190 }
01191
01192 return det;
01193 }
01194
01195 template <class T> Matrix<T>
01196 Inverse (Matrix<T> m)
01197 {
01198 return m.Inv();
01199 }
01200
01201 template <class T> Matrix<T>
01202 Transpose (Matrix<T> m)
01203 {
01204 return ~m;
01205 }
01206
01207 template <class T> T
01208 Det (Matrix<T> m)
01209 {
01210 return m.Det();
01211 }
01212
01213
01214 #endif // MATRIX_H