开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
matrix.hpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：24k
源码类别：

生物技术

开发平台：

C/C++

matrix.hpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: matrix.hpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/04/16 15:28:41  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [CATCHUP_003] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. #ifndef UTIL_MATH___MATRIX__HPP
  11. #define UTIL_MATH___MATRIX__HPP
  13. /*  $Id: matrix.hpp,v 1000.1 2004/04/16 15:28:41 gouriano Exp $
  14.  * ===========================================================================
  15.  *
  16.  *                            PUBLIC DOMAIN NOTICE
  17.  *               National Center for Biotechnology Information
  18.  *
  19.  *  This software/database is a "United States Government Work" under the
  20.  *  terms of the United States Copyright Act.  It was written as part of
  21.  *  the author's official duties as a United States Government employee and
  22.  *  thus cannot be copyrighted.  This software/database is freely available
  23.  *  to the public for use. The National Library of Medicine and the U.S.
  24.  *  Government have not placed any restriction on its use or reproduction.
  25.  *
  26.  *  Although all reasonable efforts have been taken to ensure the accuracy
  27.  *  and reliability of the software and data, the NLM and the U.S.
  28.  *  Government do not and cannot warrant the performance or results that
  29.  *  may be obtained by using this software or data. The NLM and the U.S.
  30.  *  Government disclaim all warranties, express or implied, including
  31.  *  warranties of performance, merchantability or fitness for any particular
  32.  *  purpose.
  33.  *
  34.  *  Please cite the author in any work or product based on this material.
  35.  *
  36.  * ===========================================================================
  37.  *
  38.  * Authors:  Mike DiCuccio
  39.  *
  40.  * File Description:
  41.  *
  42.  */
  44. #include <corelib/ncbistr.hpp>
  45. #include <util/math/promote.hpp>
  46. #include <vector>
  49. BEGIN_NCBI_SCOPE
  51. template <class T> class CNcbiMatrix;
  54. template <class T>
  55. class CNcbiMatrix
  56. {
  57. public:
  58.     typedef vector<T>                      TData;
  59.     typedef typename TData::iterator       iterator;
  60.     typedef typename TData::const_iterator const_iterator;
  62.     CNcbiMatrix();
  63.     CNcbiMatrix(size_t r, size_t c, T val = T(0));
  65.     // make this matrix an identity matrix of a given size
  66.     void Identity(size_t size);
  68.     // make this matrix an identity matrix
  69.     void Identity();
  71.     // make this matrix a diagonal matrix of a given size, with a given value
  72.     // on the diagonal
  73.     void Diagonal(size_t size, T val);
  75.     // make this matrix a diagonal matrix, with a given value on the diagonal
  76.     void Diagonal(T val);
  78.     // transpose this matrix
  79.     void Transpose();
  81.     // resize this matrix, filling the empty cells with a known value
  82.     void Resize(size_t i, size_t j, T val = T(0));
  84.     // swap two rows in the matrix
  85.     void SwapRows(size_t i, size_t j);
  87.     // get the number of rows in this matrix
  88.     size_t GetRows() const;
  90.     // get the number of columns in this matrix
  91.     size_t GetCols() const;
  93.     // retrieve the data associated with this matrix
  94.     TData&         GetData();
  95.     const TData&   GetData() const;
  97.     // iterators
  98.     iterator begin();
  99.     iterator end();
  100.     const_iterator begin() const;
  101.     const_iterator end() const;
  103.     // set all values in the matrix to a known value
  104.     void Set(T val);
  106.     // swap two matrices efficiently
  107.     void Swap(CNcbiMatrix<T>& M);
  109.     // operator[] for raw data indexing
  110.     const T&    operator[] (size_t i) const;
  111.     T&          operator[] (size_t i);
  113.     // operator() for row/column indexing
  114.     const T&    operator() (size_t i, size_t j) const;
  115.     T&          operator() (size_t i, size_t j);
  117. protected:
  119.     // the data strip we use
  120.     TData  m_Data;
  122.     // size of this matrix
  123.     size_t m_Rows;
  124.     size_t m_Cols;
  125. };
  129. ///
  130. /// global operators
  131. ///
  133. ///
  134. /// stream input.
  135. ///
  136. /// One line per row, with entries readable by successive calls
  137. /// to operator>>.  For doubles, this just means they're separated
  138. /// by whitespace.
  139. template <class T>
  140. inline CNcbiIstream&
  141. operator>> (CNcbiIstream& is, CNcbiMatrix<T>& M);
  143. ///
  144. /// stream output.
  145. ///
  146. /// One line per row, with entries separated by a single space
  147. template <class T>
  148. inline CNcbiOstream&
  149. operator<< (CNcbiOstream& os, const CNcbiMatrix<T>& M);
  151. ///
  152. /// global addition: matrix + matrix
  153. ///
  154. template <class T, class U>
  155. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  156. operator+ (const CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  158. ///
  159. /// global subtraction: matrix - matrix
  160. ///
  161. template <class T, class U>
  162. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  163. operator- (const CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  166. ///
  167. /// global multiplication: matrix * scalar
  168. /// this is a hack, as MSVC doesn't support partial template specialization
  169. /// we provide implementations for a number of popular types
  170. ///
  171. template <class T>
  172. inline CNcbiMatrix< NCBI_PROMOTE(T, size_t) >
  173. operator* (const CNcbiMatrix<T>&, size_t);
  175. template <class U>
  176. inline CNcbiMatrix< NCBI_PROMOTE(size_t, U) >
  177. operator* (size_t, const CNcbiMatrix<U>&);
  180. template <class T>
  181. inline CNcbiMatrix< NCBI_PROMOTE(T, float) >
  182. operator* (const CNcbiMatrix<T>&, float);
  184. template <class U>
  185. inline CNcbiMatrix< NCBI_PROMOTE(float, U) >
  186. operator* (float, const CNcbiMatrix<U>&);
  189. template <class T>
  190. inline CNcbiMatrix< NCBI_PROMOTE(T, double) >
  191. operator* (const CNcbiMatrix<T>&, double);
  193. template <class U>
  194. inline CNcbiMatrix< NCBI_PROMOTE(double, U) >
  195. operator* (double, const CNcbiMatrix<U>&);
  198. ///
  199. /// global multiplication: matrix * matrix
  200. ///
  201. template <class T, class U>
  202. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  203. operator* (const CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  205. ///
  206. /// global multiplication: matrix * vector
  207. ///
  208. template <class T, class U>
  209. inline vector< NCBI_PROMOTE(T,U) >
  210. operator* (const CNcbiMatrix<T>&, const vector<U>&);
  212. ///
  213. /// global multiplication: vector * matrix
  214. ///
  215. template <class T, class U>
  216. inline vector< NCBI_PROMOTE(T,U) >
  217. operator* (const vector<T>&, const CNcbiMatrix<U>&);
  219. ///
  220. /// global addition: matrix += matrix
  221. ///
  222. template <class T, class U>
  223. inline CNcbiMatrix<T>&
  224. operator+= (CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  226. ///
  227. /// global subtraction: matrix -= matrix
  228. ///
  229. template <class T, class U>
  230. inline CNcbiMatrix<T>&
  231. operator-= (CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  233. ///
  234. /// global multiplication: matrix *= matrix
  235. ///
  236. template <class T, class U>
  237. inline CNcbiMatrix<T>&
  238. operator*= (CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  240. ///
  241. /// global multiplication: matrix *= vector
  242. ///
  243. template <class T, class U>
  244. inline vector<T>&
  245. operator*= (CNcbiMatrix<T>&, const vector<U>&);
  247. ///
  248. /// global multiplication: matrix *= scalar
  249. ///
  250. template <class T>
  251. inline CNcbiMatrix<T>&
  252. operator*= (CNcbiMatrix<T>&, T);
  254. ///
  255. /// global division: matrix /= matrix
  256. ///
  257. template <class T>
  258. inline CNcbiMatrix<T>&
  259. operator/= (CNcbiMatrix<T>&, T);
  261. ///
  262. /// global comparison: matrix == matrix
  263. ///
  264. template <class T, class U>
  265. inline bool
  266. operator== (const CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  268. ///
  269. /// global comparison: matrix != matrix
  270. ///
  271. template <class T, class U>
  272. inline bool
  273. operator!= (const CNcbiMatrix<T>&, const CNcbiMatrix<U>&);
  276. /////////////////////////////////////////////////////////////////////////////
  277. //
  278. // Inline Methods
  279. //
  281. ///
  282. /// default ctor
  283. template <class T>
  284. inline CNcbiMatrix<T>::CNcbiMatrix()
  285.     : m_Rows(0),
  286.       m_Cols(0)
  287. {
  288. }
  291. template <class T>
  292. inline CNcbiMatrix<T>::CNcbiMatrix(size_t r, size_t c, T val)
  293.     : m_Rows(r),
  294.       m_Cols(c)
  295. {
  296.     m_Data.resize(r*c, val);
  297. }
  300. template <class T>
  301. inline size_t CNcbiMatrix<T>::GetRows() const
  302. {
  303.     return m_Rows;
  304. }
  307. template <class T>
  308. inline size_t CNcbiMatrix<T>::GetCols() const
  309. {
  310.     return m_Cols;
  311. }
  314. template <class T>
  315. inline typename CNcbiMatrix<T>::TData& CNcbiMatrix<T>::GetData()
  316. {
  317.     return m_Data;
  318. }
  321. template <class T>
  322. inline const typename CNcbiMatrix<T>::TData& CNcbiMatrix<T>::GetData() const
  323. {
  324.     return m_Data;
  325. }
  327. template <class T>
  328. inline typename CNcbiMatrix<T>::iterator CNcbiMatrix<T>::begin()
  329. {
  330.     return m_Data.begin();
  331. }
  334. template <class T>
  335. inline typename CNcbiMatrix<T>::iterator CNcbiMatrix<T>::end()
  336. {
  337.     return m_Data.end();
  338. }
  341. template <class T>
  342. inline typename CNcbiMatrix<T>::const_iterator CNcbiMatrix<T>::begin() const
  343. {
  344.     return m_Data.begin();
  345. }
  348. template <class T>
  349. inline typename CNcbiMatrix<T>::const_iterator CNcbiMatrix<T>::end() const
  350. {
  351.     return m_Data.end();
  352. }
  355. template <class T>
  356. inline const T&
  357. CNcbiMatrix<T>::operator[] (size_t i) const
  358. {
  359.     _ASSERT(i >= 0  &&  i < m_Data.size());
  360.     return m_Data[i];
  361. }
  364. template <class T>
  365. inline T& CNcbiMatrix<T>::operator[] (size_t i)
  366. {
  367.     _ASSERT(i >= 0  &&  i < m_Data.size());
  368.     return m_Data[i];
  369. }
  372. template <class T>
  373. inline const T& CNcbiMatrix<T>::operator() (size_t i, size_t j) const
  374. {
  375.     _ASSERT(i >= 0  &&  i < m_Rows);
  376.     _ASSERT(j >= 0  &&  j < m_Cols);
  378.     return m_Data[i * m_Cols + j];
  379. } 
  382. template <class T>
  383. inline T& CNcbiMatrix<T>::operator() (size_t i, size_t j)
  384. {
  385.     _ASSERT(i >= 0  &&  i < m_Rows);
  386.     _ASSERT(j >= 0  &&  j < m_Cols);
  388.     return m_Data[i * m_Cols + j];
  389. }
  392. template <class T>
  393. inline void CNcbiMatrix<T>::Resize(size_t new_rows, size_t new_cols, T val)
  394. {
  396.     if (new_cols == m_Cols && new_rows >= m_Rows) {
  397.         // common special case that can easily be handled efficiently
  398.         m_Data.resize(new_rows * new_cols, val);
  399.     } else {
  400.         /// hack: we just make a new strip and copy things correctly
  401.         /// there is a faster way to do this
  402.         TData new_data(new_rows * new_cols, val);
  403.         size_t i = min(new_rows, m_Rows);
  404.         size_t j = min(new_cols, m_Cols);
  405.         
  406.         for (size_t r = 0;  r < i;  ++r) {
  407.             for (size_t c = 0;  c < j;  ++c) {
  408.                 new_data[r * new_cols + c] = m_Data[r * m_Cols + c];
  409.             }
  410.         }
  412.         new_data.swap(m_Data);
  413.     }
  414.     m_Rows = new_rows;
  415.     m_Cols = new_cols;
  416. }
  419. template <class T>
  420. inline void CNcbiMatrix<T>::Set(T val)
  421. {
  422.     m_Data.clear();
  423.     m_Data.resize(m_Rows * m_Cols, val);
  424. }
  427. template <class T>
  428. inline void CNcbiMatrix<T>::Identity()
  429. {
  430.     _ASSERT(m_Rows == m_Cols);
  431.     Set(T(0));
  433.     for (size_t i = 0;  i < m_Rows;  ++i) {
  434.         m_Data[i * m_Cols + i] = T(1);
  435.     }
  436. }
  439. template <class T>
  440. inline void CNcbiMatrix<T>::Identity(size_t size)
  441. {
  442.     Resize(size, size);
  443.     Set(T(0));
  445.     for (size_t i = 0;  i < m_Rows;  ++i) {
  446.         m_Data[i * m_Cols + i] = T(1);
  447.     }
  448. }
  451. template <class T>
  452. inline void CNcbiMatrix<T>::Diagonal(T val)
  453. {
  454.     _ASSERT(m_Rows == GetCols());
  455.     Set(T(0));
  457.     for (size_t i = 0;  i < m_Rows;  ++i) {
  458.         m_Data[i * m_Cols + i] = val;
  459.     }
  460. }
  463. template <class T>
  464. inline void CNcbiMatrix<T>::Diagonal(size_t size, T val)
  465. {
  466.     Resize(size, size);
  467.     Set(T(0));
  469.     for (size_t i = 0;  i < m_Rows;  ++i) {
  470.         m_Data[i * m_Cols + i] = val;
  471.     }
  472. }
  475. template <class T>
  476. inline void CNcbiMatrix<T>::Transpose()
  477. {
  478.     TData new_data(m_Cols * m_Rows);
  480.     for (size_t i = 0;  i < m_Rows;  ++i) {
  481.         for (size_t j = 0;  j < m_Cols;  ++j) {
  482.             new_data[j * m_Cols + i] = m_Data[i * m_Cols + j];
  483.         }
  484.     }
  486.     m_Data.swap(new_data);
  487.     swap(m_Rows, m_Cols);
  488. }
  491. template <class T>
  492. inline void CNcbiMatrix<T>::SwapRows(size_t i, size_t j)
  493. {
  494.     size_t i_offs = i * m_Cols;
  495.     size_t j_offs = j * m_Cols;
  497.     for (size_t c = 0;  c < m_Cols;  ++c) {
  498.         swap(m_Data[i_offs + c], m_Data[j_offs + c] );
  499.     }
  500. }
  503. template <class T>
  504. inline void CNcbiMatrix<T>::Swap(CNcbiMatrix<T>& M)
  505. {
  506.     m_Data.swap(M.m_Data);
  507.     swap(m_Cols, M.m_Cols);
  508.     swap(m_Rows, M.m_Rows);
  509. }
  512. ///
  513. /// global operators
  514. ///
  517. ///
  518. /// addition
  519. ///
  521. template <class T, class U>
  522. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  523. operator+ (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  524. {
  525.     _ASSERT(m0.GetRows() == m1.GetRows());
  526.     _ASSERT(m0.GetCols() == m1.GetCols());
  528.     CNcbiMatrix<NCBI_PROMOTE(T,U)> res(m0.GetRows(), m0.GetCols());
  530.     typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_iter = res.begin();
  531.     typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_end  = res.end();
  532.     typename CNcbiMatrix<T>::const_iterator           m0_iter  = m0.begin();
  533.     typename CNcbiMatrix<U>::const_iterator           m1_iter  = m1.begin();
  535.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter, ++m1_iter) {
  536.         *res_iter = *m0_iter + *m1_iter;
  537.     }
  539.     return res;
  540. }
  543. template <class T, class U>
  544. inline CNcbiMatrix<T>&
  545. operator+= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  546. {
  547.     _ASSERT(m0.GetRows() == m1.GetRows());
  548.     _ASSERT(m0.GetCols() == m1.GetCols());
  550.     typename CNcbiMatrix<T>::iterator           m0_iter  = m0.begin();
  551.     typename CNcbiMatrix<T>::iterator           m0_end   = m0.end();
  552.     typename CNcbiMatrix<U>::const_iterator     m1_iter  = m1.begin();
  554.     size_t mod = (m0.GetRows() * m0.GetCols()) % 4;
  555.     for ( ;  mod;  --mod, ++m0_iter, ++m1_iter) {
  556.         *m0_iter += *m1_iter;
  557.     }
  559.     for ( ;  m0_iter != m0_end;  ) {
  560.         *m0_iter++ += *m1_iter++;
  561.         *m0_iter++ += *m1_iter++;
  562.         *m0_iter++ += *m1_iter++;
  563.         *m0_iter++ += *m1_iter++;
  564.     }
  566.     return m0;
  567. }
  570. ///
  571. /// subtraction
  572. ///
  574. template <class T, class U>
  575. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  576. operator- (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  577. {
  578.     _ASSERT(m0.GetRows() == m1.GetRows());
  579.     _ASSERT(m0.GetCols() == m1.GetCols());
  581.     CNcbiMatrix< NCBI_PROMOTE(T,U) > res(m0.GetRows(), m0.GetCols());
  582.     typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_iter = res.begin();
  583.     typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_end  = res.end();
  584.     typename CNcbiMatrix<T>::const_iterator           m0_iter  = m0.begin();
  585.     typename CNcbiMatrix<U>::const_iterator           m1_iter  = m1.begin();
  587.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter, ++m1_iter) {
  588.         *res_iter = *m0_iter - *m1_iter;
  589.     }
  591.     return res;
  592. }
  595. template <class T, class U>
  596. inline CNcbiMatrix<T>&
  597. operator-= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  598. {
  599.     _ASSERT(m0.GetRows() == m1.GetRows());
  600.     _ASSERT(m0.GetCols() == m1.GetCols());
  602.     typename CNcbiMatrix<T>::iterator           m0_iter  = m0.begin();
  603.     typename CNcbiMatrix<T>::iterator           m0_end   = m0.end();
  604.     typename CNcbiMatrix<U>::const_iterator     m1_iter  = m1.begin();
  606.     for ( ;  m0_iter != m0_end;  ++m0_iter, ++m1_iter) {
  607.         *m0_iter -= *m1_iter;
  608.     }
  610.     return m0;
  611. }
  614. ///
  615. /// multiplication
  616. ///
  617. template <class T>
  618. inline CNcbiMatrix< NCBI_PROMOTE(T, size_t) >
  619. operator* (const CNcbiMatrix<T>& m0, size_t val)
  620. {
  621.     CNcbiMatrix< NCBI_PROMOTE(T, size_t) > res(m0.GetRows(), m0.GetCols());
  623.     typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
  624.     typename CNcbiMatrix<T>::iterator       res_end  = res.end();
  625.     typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();
  627.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  628.         *res_iter = *m0_iter * val;
  629.     }
  631.     return res;
  632. }
  635. template <class U>
  636. inline CNcbiMatrix< NCBI_PROMOTE(size_t, U) >
  637. operator* (size_t val, const CNcbiMatrix<U>& m0)
  638. {
  639.     CNcbiMatrix< NCBI_PROMOTE(U, size_t) > res(m0.GetRows(), m0.GetCols());
  641.     typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
  642.     typename CNcbiMatrix<U>::iterator       res_end  = res.end();
  643.     typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();
  645.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  646.         *res_iter = *m0_iter * val;
  647.     }
  649.     return res;
  650. }
  653. template <class T>
  654. inline CNcbiMatrix< NCBI_PROMOTE(T, float) >
  655. operator* (const CNcbiMatrix<T>& m0, float val)
  656. {
  657.     CNcbiMatrix< NCBI_PROMOTE(T, float) > res(m0.GetRows(), m0.GetCols());
  659.     typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
  660.     typename CNcbiMatrix<T>::iterator       res_end  = res.end();
  661.     typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();
  663.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  664.         *res_iter = *m0_iter * val;
  665.     }
  667.     return res;
  668. }
  671. template <class U>
  672. inline CNcbiMatrix< NCBI_PROMOTE(float, U) >
  673. operator* (float val, const CNcbiMatrix<U>& m0)
  674. {
  675.     CNcbiMatrix< NCBI_PROMOTE(U, float) > res(m0.GetRows(), m0.GetCols());
  677.     typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
  678.     typename CNcbiMatrix<U>::iterator       res_end  = res.end();
  679.     typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();
  681.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  682.         *res_iter = *m0_iter * val;
  683.     }
  685.     return res;
  686. }
  689. template <class T>
  690. inline CNcbiMatrix< NCBI_PROMOTE(T, double) >
  691. operator* (const CNcbiMatrix<T>& m0, double val)
  692. {
  693.     CNcbiMatrix< NCBI_PROMOTE(T, double) > res(m0.GetRows(), m0.GetCols());
  695.     typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
  696.     typename CNcbiMatrix<T>::iterator       res_end  = res.end();
  697.     typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();
  699.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  700.         *res_iter = *m0_iter * val;
  701.     }
  703.     return res;
  704. }
  707. template <class U>
  708. inline CNcbiMatrix< NCBI_PROMOTE(double, U) >
  709. operator* (double val, const CNcbiMatrix<U>& m0)
  710. {
  711.     CNcbiMatrix< NCBI_PROMOTE(U, double) > res(m0.GetRows(), m0.GetCols());
  713.     typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
  714.     typename CNcbiMatrix<U>::iterator       res_end  = res.end();
  715.     typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();
  717.     for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
  718.         *res_iter = *m0_iter * val;
  719.     }
  721.     return res;
  722. }
  725. template <class T>
  726. inline CNcbiMatrix<T>&
  727. operator*= (CNcbiMatrix<T>& m0, T val)
  728. {
  729.     typename CNcbiMatrix<T>::iterator m0_iter  = m0.begin();
  730.     typename CNcbiMatrix<T>::iterator m0_end  = m0.end();
  732.     for ( ;  m0_iter != m0_end;  ++m0_iter) {
  733.         *m0_iter *= val;
  734.     }
  736.     return m0;
  737. }
  740. template <class T>
  741. inline CNcbiMatrix<T>&
  742. operator/= (CNcbiMatrix<T>& m0, T val)
  743. {
  744.     val = T(1/val);
  746.     typename CNcbiMatrix<T>::iterator m0_iter  = m0.begin();
  747.     typename CNcbiMatrix<T>::iterator m0_end  = m0.end();
  749.     for ( ;  m0_iter != m0_end;  ++m0_iter) {
  750.         *m0_iter *= val;
  751.     }
  753.     return m0;
  754. }
  757. template <class T, class U>
  758. inline vector< NCBI_PROMOTE(T,U) >
  759. operator* (const CNcbiMatrix<T>& m, const vector<U>& v)
  760. {
  761.     _ASSERT(m.GetCols() == v.size());
  763.     typedef NCBI_PROMOTE(T,U) TPromote;
  764.     vector< TPromote > res(m.GetRows(), TPromote(0));
  766.     for (size_t r = 0;  r < m.GetRows();  ++r) {
  767.         for (size_t i = 0;  i < m.GetCols();  ++i) {
  768.             res[r] += m(r,i) * v[i];
  769.         }
  770.     }
  772.     return res;
  773. }
  776. template <class T, class U>
  777. inline vector< NCBI_PROMOTE(T,U) >
  778. operator* (const vector<T>& v, const CNcbiMatrix<U>& m)
  779. {
  780.     _ASSERT(m.GetRows() == v.size());
  782.     typedef NCBI_PROMOTE(T,U) TPromote;
  783.     vector<TPromote> res(m.GetCols(), TPromote(0));
  785.     for (size_t c = 0;  c < m.GetCols();  ++c) {
  786.         for (size_t i = 0;  i < m.GetRows();  ++i) {
  787.             res[c] += m(i,c) * v[i];
  788.         }
  789.     }
  791.     return res;
  792. }
  795. template <class T, class U>
  796. inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
  797. operator* (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  798. {
  799.     _ASSERT(m0.GetCols() == m1.GetRows());
  801.     typedef NCBI_PROMOTE(T,U) TPromote;
  802.     CNcbiMatrix< TPromote > res(m0.GetRows(), m1.GetCols(), TPromote(0));
  804. #if 1
  805.     for (size_t r = 0; r < m0.GetRows();  ++r) {
  806.         for (size_t c = 0;  c < m1.GetCols();  ++c) {
  807.             for (size_t i = 0;  i < m0.GetCols();  ++i) {
  808.                 res(r,c) += m0(r,i) * m1(i,c);
  809.             }
  810.         }
  811.     }
  812. #endif
  815. #if 0
  816.     const vector<T>& lhs = m0.GetData();
  817.     const vector<U>& rhs = m1.GetData();
  818.     size_t col_mod = m0.GetCols() % 4;
  820.     for (size_t r = 0;  r < m0.GetRows();  ++r) {
  821.         size_t r_offs = r * m0.GetCols();
  823.         for (size_t c = 0;  c < m1.GetCols();  ++c) {
  824.             size_t i=0;
  825.             T t0 = 0;
  826.             T t1 = 0;
  827.             T t2 = 0;
  828.             T t3 = 0;
  830.             switch(col_mod) {
  831.             default:
  832.             case 0:
  833.                 break;
  835.             case 3:
  836.                 t0 += m0.GetData()[r_offs + 2] * m1(2,c);
  837.                 ++i;
  838.             case 2:
  839.                 t0 += m0.GetData()[r_offs + 1] * m1(1,c);
  840.                 ++i;
  841.             case 1:
  842.                 t0 += m0.GetData()[r_offs + 0] * m1(0,c);
  843.                 ++i;
  844.                 break;
  845.             }
  847.             for (;  i < m0.GetCols();  i += 2) {
  848.                 t0 += lhs[r_offs + i + 0] * m1(i + 0, c);
  849.                 t1 += lhs[r_offs + i + 1] * m1(i + 1, c);
  850.                 //t2 += lhs[r_offs + i + 2] * m1(i + 2, c);
  851.                 //t3 += lhs[r_offs + i + 3] * m1(i + 3, c);
  852.             }
  853.             res(r,c) = t0 + t1 + t2 + t3;
  854.         }
  855.     }
  856. #endif
  858.     return res;
  859. }
  862. template <class T, class U>
  863. inline CNcbiMatrix<T>&
  864. operator*= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  865. {
  866.     _ASSERT(m0.GetCols() == m1.GetRows());
  868.     m0 = m0 * m1;
  869.     return m0;
  870. }
  873. ///
  874. /// comparators
  875. ///
  877. template <class T, class U>
  878. inline bool
  879. operator== (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  880. {
  881.     if (m0.GetRows() != m1.GetRows()) {
  882.         return false;
  883.     }
  884.     if (m0.GetCols() != m1.GetCols()) {
  885.         return false;
  886.     }
  888.     for (size_t i = 0;  i < m0.GetData().size();  ++i) {
  889.         if (m0[i] != m1[i]) {
  890.             return false;
  891.         }
  892.     }
  894.     return true;
  895. }
  898. template <class T, class U>
  899. inline bool
  900. operator!= (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
  901. {
  902.     return !(m0 == m1);
  903. }
  906. ///
  907. /// stream output
  908. template <class T>
  909. inline CNcbiOstream& operator<<(CNcbiOstream& os, const CNcbiMatrix<T>& M)
  910. {
  911.     for (size_t i = 0;  i < M.GetRows();  ++i) {
  912.         for (size_t j = 0;  j < M.GetCols();  ++j) {
  913.             if (j > 0) {
  914.                 os << " ";
  915.             }
  916.             os << M(i, j);
  917.         }
  918.         os << NcbiEndl;
  919.     }
  920.     return os;
  921. }
  924. ///
  925. /// stream input
  926. template <class T>
  927. inline CNcbiIstream& operator>>(CNcbiIstream& is, CNcbiMatrix<T>& M)
  928. {
  929.     CNcbiMatrix<T> A;
  930.     string line;
  931.     vector<T> row;
  932.     T entry;
  933.     int linenum = 0;
  934.     while(getline(is, line)) {
  935.         linenum++;
  936.         if (line.empty() || line[0] == '#') {
  937.             continue;
  938.         }
  939.         CNcbiIstrstream iss(line.c_str(), line.size());
  940.         row.clear();
  941.         while(1) {
  942.             iss >> entry;
  943.             if (!iss) {
  944.                 break;
  945.             }
  946.             row.push_back(entry);
  947.         }
  948.         if (A.GetCols() == 0) {
  949.             A.Resize(A.GetCols(), row.size());
  950.         }
  951.         if (row.size() != A.GetCols()) {
  952.             NCBI_THROW(CException, eUnknown,
  953.                        "error at line " +
  954.                        NStr::IntToString(linenum) + ": expected " +
  955.                        NStr::IntToString(A.GetCols()) + " columns; got" +
  956.                        NStr::IntToString(row.size()));
  957.         }
  958.         A.Resize(A.GetRows() + 1, A.GetCols());
  959.         for (int i = 0;  i < A.GetCols();  ++i) {
  960.             A(A.GetRows() - 1, i) = row[i];
  961.         }
  962.     }
  963.     M.Swap(A);
  964.     return is;
  965. }
  969. END_NCBI_SCOPE
  971. /*
  972.  * ===========================================================================
  973.  * $Log: matrix.hpp,v $
  974.  * Revision 1000.1  2004/04/16 15:28:41  gouriano
  975.  * PRODUCTION: UPGRADED [CATCHUP_003] Dev-tree R1.6
  976.  *
  977.  * Revision 1.6  2004/03/11 17:31:42  jcherry
  978.  * Sped up adding of rows (greatly speeds reading from streams)
  979.  *
  980.  * Revision 1.5  2004/03/10 14:13:29  dicuccio
  981.  * Corrected include - use ncbistr instead of ncbistre
  982.  *
  983.  * Revision 1.4  2004/02/10 21:17:54  jcherry
  984.  * Made operator>> work for a matrix of any contained type that defines
  985.  * operator>> (not just doubles)
  986.  *
  987.  * Revision 1.3  2004/02/10 14:29:51  jcherry
  988.  * Removed wayward use of protected data in operator>>
  989.  *
  990.  * Revision 1.2  2004/02/10 14:00:40  dicuccio
  991.  * Added comments.  Fixed spacing issues.  Changed extraction operator to be
  992.  * non-friend function
  993.  *
  994.  * Revision 1.1  2004/02/09 17:34:11  dicuccio
  995.  * Moved from gui/math
  996.  *
  997.  * Revision 1.4  2004/02/02 15:55:49  jcherry
  998.  * Added Swap() method and stream reader
  999.  *
  1000.  * Revision 1.3  2004/01/29 22:38:03  jcherry
  1001.  * Fixed typo
  1002.  *
  1003.  * Revision 1.2  2004/01/29 21:30:17  jcherry
  1004.  * Added stream output operator
  1005.  *
  1006.  * Revision 1.1  2004/01/23 13:54:15  dicuccio
  1007.  * Initial revision
  1008.  *
  1009.  * ===========================================================================
  1010.  */
  1012. #endif  // UTIL_MATH___MATRIX__HPP