home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / matrix.hpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  12.3 KB  |  392 lines
  1. //===================================================================
  2. // matrix.hpp
  3. //
  4. // Version 1.0
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1999 WordenWare
  12. //
  13. // Created:  April 10, 1999
  14. // Revised:  
  15. //===================================================================
  16.  
  17. #ifndef _MATRIX_HPP_
  18. #define _MATRIX_HPP_
  19.  
  20. #include <algorithm>
  21. #include <cstdlib>
  22. #include <functional>
  23. #include <numeric>
  24. #include <cmath>
  25.  
  26. #include "funcobj.hpp"
  27. #include "vector.hpp"
  28.  
  29. NUM_BEGIN
  30.  
  31. template<class T>
  32. class Matrix
  33. //-------------------------------------------------------------------
  34. // mathematical matrix object
  35. //-------------------------------------------------------------------
  36. {
  37. public:
  38.     typedef Vector<T>* iterator;
  39.     typedef const Vector<T>* const_iterator;
  40.     typedef size_t size_type;
  41.  
  42.     Matrix()
  43.     //---------------------------------------------------------------
  44.     // Create an empty matrix.
  45.     //---------------------------------------------------------------
  46.     : _row_capacity(0), _col_capacity(0), _rows(0), _cols(0), _data(NULL) {}
  47.  
  48.     Matrix(size_type m, size_type n)
  49.     //---------------------------------------------------------------
  50.     // Create a m x n matrix
  51.     //---------------------------------------------------------------
  52.     : _row_capacity(0), _col_capacity(0), _rows(m), _cols(n)
  53.     {
  54.         _data = _allocate(m, n, T(0));
  55.     }
  56.  
  57.     Matrix(const Matrix<T>& other)
  58.     //---------------------------------------------------------------
  59.     // Create a matrix by copying other.
  60.     //---------------------------------------------------------------
  61.     {
  62.         size_type r = other.rows();
  63.         size_type c = other.columns();
  64.         _row_capacity = r;
  65.         _col_capacity = c,
  66.         _rows = r;
  67.         _cols = c,
  68.         _data = _allocate(rowCapacity(), columnCapacity(), T(0));
  69.         std::copy(other.row_begin(), other.row_end(), _data);
  70.     }
  71.  
  72.     Matrix(const_iterator first, const_iterator last)
  73.     //---------------------------------------------------------------
  74.     // Create a matrix with initial rows [first, last).
  75.     //---------------------------------------------------------------
  76.     : _rows(_length(first, last)), _cols(first->size())
  77.     {
  78.         size_type r = _length(first, last);
  79.         size_type c = first->size();
  80.         _row_capacity = r;
  81.         _col_capacity = c,
  82.         _rows = r;
  83.         _cols = c,
  84.         _data = _allocate(rowCapacity(), columnCapacity(), T(0));
  85.         std::copy(first, last, _data);
  86.     }
  87.  
  88.     ~Matrix()
  89.     //---------------------------------------------------------------
  90.     // Destroy this matrix.
  91.     //---------------------------------------------------------------
  92.     { _destroy(); }
  93.  
  94.     const_iterator row_begin() const
  95.     //---------------------------------------------------------------
  96.     // An iterator positioned at the first row of this matrix.
  97.     //---------------------------------------------------------------
  98.     { return (const_iterator)_data; }
  99.  
  100.     const_iterator row_end() const
  101.     //---------------------------------------------------------------
  102.     // An interator positioned one position past the last row of this
  103.     // matrix.
  104.     //---------------------------------------------------------------
  105.     { return (const_iterator)(_data + _rows); }
  106.  
  107.     iterator row_begin()
  108.     //---------------------------------------------------------------
  109.     // An iterator positioned at the first row of this matrix.
  110.     //---------------------------------------------------------------
  111.     { return _data; }
  112.  
  113.     iterator row_end()
  114.     //---------------------------------------------------------------
  115.     // An interator positioned one position past the last row of this
  116.     // matrix.
  117.     //---------------------------------------------------------------
  118.     { return (_data + _rows); }
  119.  
  120.     size_type rows() const
  121.     //---------------------------------------------------------------
  122.     // Access the number of rows of this matrix.
  123.     //---------------------------------------------------------------
  124.     { return _rows; }
  125.  
  126.     size_type columns() const
  127.     //---------------------------------------------------------------
  128.     // Access the number of columns of this matrix.
  129.     //---------------------------------------------------------------
  130.     { return _cols; }
  131.  
  132.     size_type rowCapacity() const
  133.     //---------------------------------------------------------------
  134.     // Access the capcity of rows of this matrix.
  135.     //---------------------------------------------------------------
  136.     { return _row_capacity; }
  137.  
  138.     size_type columnCapacity() const
  139.     //---------------------------------------------------------------
  140.     // Access the number of columns of this matrix.
  141.     //---------------------------------------------------------------
  142.     { return _col_capacity; }
  143.  
  144.     const Matrix<T>& operator=(const Matrix<T>& rhs)
  145.     //---------------------------------------------------------------
  146.     // Assign rhs to this matrix.
  147.     //---------------------------------------------------------------
  148.     { 
  149.         if(this != &rhs){
  150.             size_type r = rhs.rows();
  151.             size_type c = rhs.columns();
  152.             if(r <= rowCapacity()){
  153.                 std::copy(rhs.row_begin(), rhs.row_end(), _data);
  154.             } else {
  155.                 _destroy();
  156.                 _data = _allocate(r, c, T(0));
  157.                 std::copy(rhs.row_begin(), rhs.row_end(), _data);
  158.                 _row_capacity = r;
  159.                 _col_capacity = c;
  160.             }
  161.             _rows = r;
  162.             _cols = c;
  163.         }
  164.         return *this;
  165.     }
  166.  
  167.     const Matrix<T>& operator+=(const Matrix<T>& rhs)
  168.     //---------------------------------------------------------------
  169.     // Assign the matrix sum of this matrix and rhs to this matrix.
  170.     //---------------------------------------------------------------
  171.     {
  172.         if(rows() != rhs.rows() || columns() != rhs.columns()){
  173.             throw Exception("Matrix<T>::operator+=", "Incompatible matrix sizes");
  174.         }
  175.         std::transform(row_begin(), row_end(), rhs.row_begin(), row_begin(), std::plus<Vector<T> >());
  176.         return *this;
  177.     }
  178.  
  179.     const Matrix<T>& operator-=(const Matrix<T>& rhs)
  180.     //---------------------------------------------------------------
  181.     // Assign the matrix difference of this matrix and rhs to this
  182.     // matrix.
  183.     //---------------------------------------------------------------
  184.     {
  185.         if(rows() != rhs.rows() || columns() != rhs.columns()){
  186.             throw Exception("Matrix<T>::operator-=", "Incompatible matrix sizes");
  187.         }
  188.         std::transform(row_begin(), row_end(), rhs.row_begin(), row_begin(), std::minus<Vector<T> >());
  189.         return *this;
  190.     }
  191.  
  192.     const Matrix<T>& operator*=(T scale)
  193.     //---------------------------------------------------------------
  194.     // Assign the scalar procuct of this matrix and scale to this
  195.     // matrix.
  196.     //---------------------------------------------------------------
  197.     {
  198.         std::transform(row_begin(), row_end(), row_begin(), ScaleValue<Vector<T> >(scale));
  199.         return *this;
  200.     }
  201.  
  202.     Vector<T>& operator[](size_type index)
  203.     //---------------------------------------------------------------
  204.     // Access the index-th row of this matrix.
  205.     //---------------------------------------------------------------
  206.     {
  207.         if(index < 0 || index >= rows()){
  208.             throw Exception("Matrix<T>::operator[]", "invalid index parameter");
  209.         }
  210.         return *(_data + index);
  211.     }
  212.  
  213.     void resize(size_type r, size_type c)
  214.     //---------------------------------------------------------------
  215.     // Resize this matrix to r x c.
  216.     //---------------------------------------------------------------
  217.     {
  218.         if(r > rowCapacity() || c > columnCapacity()){
  219.             iterator tmp = _allocate(r, c, T(0));
  220.             std::copy(_data, _data + rows(), tmp);
  221.             _destroy();
  222.             _data = tmp;
  223.             _row_capacity = r;
  224.             _col_capacity = c;
  225.         }
  226.         _rows = r;
  227.         _cols = c;
  228.     }
  229.  
  230. private:
  231.     iterator _allocate(size_type r, size_type c, const T& init = T(0))
  232.     //---------------------------------------------------------------
  233.     // Create r x c elements all initialized to init.
  234.     //---------------------------------------------------------------
  235.     {
  236.         iterator ret = new Vector<T>[r];
  237.         for(iterator p = ret; p != ret + r; ++p){
  238.             p->resize(c);
  239.             std::fill(p->begin(), p->end(), init);
  240.         }
  241.         return ret;
  242.     }
  243.  
  244.     void _destroy()
  245.     //---------------------------------------------------------------
  246.     // Destroy all the elements in this Vector.
  247.     //---------------------------------------------------------------
  248.     {
  249.         delete [] _data;
  250.         _row_capacity = 0;
  251.         _rows = 0;
  252.         _col_capacity = 0;
  253.         _cols = 0;
  254.         _data = NULL;
  255.     }
  256.  
  257.     size_type _length(const_iterator first, const_iterator last)
  258.     //---------------------------------------------------------------
  259.     // Return the number of elements in [first, last).
  260.     //---------------------------------------------------------------
  261.     {
  262.         size_type ret;
  263.         for(ret = 0; first != last; ++ret, ++first);
  264.         return ret;
  265.     }
  266.  
  267.     size_type _row_capacity;
  268.     size_type _col_capacity;
  269.     size_type _rows;         // number of rows
  270.     size_type _cols;         // number of columns
  271.     Vector<T>* _data;  // matrix rows
  272. };
  273.  
  274. template<class T>
  275. Matrix<T> operator+(const Matrix<T>& lhs, const Matrix<T>& rhs)
  276. //-------------------------------------------------------------------
  277. // Matrix sum of lhs and rhs.
  278. //-------------------------------------------------------------------
  279. {
  280.     if((lhs.columns() != rhs.columns()) ||
  281.        (lhs.rows() != rhs.rows())){
  282.         throw exception("operator+(Matrix<T>, Vector<T>)", "Incompatible sizes");
  283.     }
  284.     Matrix<T> ret(lhs);
  285.     return ret += rhs;
  286. }
  287.  
  288. template<class T>
  289. Matrix<T> operator-(const Matrix<T>& lhs, const Matrix<T>& rhs)
  290. //-------------------------------------------------------------------
  291. // Matrix sum of lhs and rhs.
  292. //-------------------------------------------------------------------
  293. {
  294.     if((lhs.columns() != rhs.columns()) ||
  295.        (lhs.rows() != rhs.rows())){
  296.         throw exception("operator-(Matrix<T>, Vector<T>)", "Incompatible sizes");
  297.     }
  298.     Matrix<T> ret(lhs);
  299.     return ret -= rhs;
  300. }
  301.  
  302. template<class T>
  303. Matrix<T> operator*(const Matrix<T>& lhs, const T& rhs)
  304. //-------------------------------------------------------------------
  305. // Scalar product of lhs and rhs.
  306. //-------------------------------------------------------------------
  307. {
  308.     Matrix<T> ret(lhs);
  309.     return ret *= rhs;
  310. }
  311.  
  312. template<class T>
  313. Matrix<T> operator*(const T& lhs, const Matrix<T>& rhs)
  314. //-------------------------------------------------------------------
  315. // Scalar product of lhs and rhs.
  316. //-------------------------------------------------------------------
  317. {
  318.     return rhs * lhs;
  319. }
  320.  
  321. template<class T>
  322. Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs)
  323. //-------------------------------------------------------------------
  324. // Dot product of lhs and rhs.
  325. //-------------------------------------------------------------------
  326. {
  327.     if(lhs.columns() != rhs.rows()){
  328.         throw exception("operator*(Matrix<T>, Matrix<T>)", "Incompatible matrix sizes");
  329.     }
  330.     Matrix<T>::size_type r = lhs.rows();
  331.     Matrix<T>::size_type c = rhs.columns();
  332.     Matrix<T> ret(r, c);
  333.     Matrix<T> t;
  334.     rhs.transpose(t);
  335.     Matrix<T>::iterator liter = lhs.row_begin();
  336.     Matrix<T>::iterator riter = ret.row_begin();
  337.     while(liter != lhs.row_end()){
  338.         Matrix<T>::iterator titer = t.row_begin();
  339.         Vector<T>::iterator iter = riter.begin();
  340.         while(titer != t.row_end()){
  341.             *iter = (*liter) * (*titer);
  342.             ++iter;
  343.             ++titer;
  344.         }
  345.         ++riter;
  346.         ++liter;
  347.     }
  348.     return ret;
  349. }
  350.  
  351. template<class T>
  352. Vector<T> operator*(const Vector<T>& lhs, const Matrix<T>& rhs)
  353. //-------------------------------------------------------------------
  354. // Dot product of lhs and rhs.  lhs is treated as a row vector.
  355. //-------------------------------------------------------------------
  356. {
  357.     if(lhs.size() != rhs.rows()){
  358.         throw exception("operator*(Vector<T>, Matrix<T>)", "Incompatible sizes");
  359.     }
  360.     return rhs.transpose() * lhs;
  361. }
  362.  
  363. template<class T>
  364. Vector<T> operator*(const Matrix<T>& lhs, const Vector<T>& rhs)
  365. //-------------------------------------------------------------------
  366. // Dot product of lhs and rhs.  rhs is treated as a column vector.
  367. //-------------------------------------------------------------------
  368. {
  369.     if(lhs.columns() != rhs.size()){
  370.         throw exception("operator*(Matrix<T>, Vector<T>)", "Incompatible sizes");
  371.     }
  372.     Vector<T> ret(lhs.rows());
  373.     Matrix<T>::iterator liter = lhs.row_begin();
  374.     Vector<T>::iterator riter = ret.begin();
  375.     while(riter != ret.end()){
  376.         *riter = rhs * (*liter);
  377.         ++riter;
  378.         ++liter;
  379.     }
  380.     return ret;
  381. }
  382.  
  383. NUM_END
  384.  
  385. #endif
  386.  
  387. //===================================================================
  388. // Revision History
  389. //
  390. // Version 1.0 - 04/10/1999 - New.
  391. //===================================================================
  392.