Mat_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // 
00003 // Authors:
00004 // - Conrad Sanderson (conradsand at ieee dot org)
00005 // 
00006 // This file is part of the Armadillo C++ library.
00007 // It is provided without any warranty of fitness
00008 // for any purpose. You can redistribute this file
00009 // and/or modify it under the terms of the GNU
00010 // Lesser General Public License (LGPL) as published
00011 // by the Free Software Foundation, either version 3
00012 // of the License or (at your option) any later version.
00013 // (see http://www.opensource.org/licenses for more info)
00014 
00015 
00016 //! \addtogroup Mat
00017 //! @{
00018 
00019 
00020 template<typename eT>
00021 inline
00022 Mat<eT>::~Mat()
00023   {
00024   arma_extra_debug_sigprint_this(this);
00025   
00026   if(n_elem > sizeof(mem_local)/sizeof(eT) )
00027     {
00028     delete [] mem;
00029     }
00030     
00031   if(arma_config::debug == true)
00032     {
00033     access::rw(mem) = 0;  // try to expose buggy code that still accesses deleted objects
00034     }
00035   
00036   isnt_supported_elem_type<eT>::check();
00037   }
00038 
00039 
00040 
00041 template<typename eT>
00042 inline
00043 Mat<eT>::Mat()
00044   : n_rows(0)
00045   , n_cols(0)
00046   , n_elem(0)
00047   //, mem(0)
00048   , mem(mem)
00049   {
00050   arma_extra_debug_sigprint_this(this);
00051   }
00052 
00053 
00054 //! construct the matrix to have user specified dimensions
00055 template<typename eT>
00056 inline
00057 Mat<eT>::Mat(const u32 in_n_rows, const u32 in_n_cols)
00058   : n_rows(0)
00059   , n_cols(0)
00060   , n_elem(0)
00061   //, mem(0)
00062   , mem(mem)
00063   {
00064   arma_extra_debug_sigprint_this(this);
00065   
00066   init(in_n_rows, in_n_cols);
00067   }
00068 
00069 
00070 
00071 //! change the matrix to have user specified dimensions (data is not preserved)
00072 template<typename eT>
00073 inline
00074 void
00075 Mat<eT>::set_size(const u32 in_n_rows, const u32 in_n_cols)
00076   {
00077   arma_extra_debug_sigprint(arma_boost::format("in_n_rows = %d, in_n_cols = %d") % in_n_rows % in_n_cols);
00078   
00079   init(in_n_rows,in_n_cols);
00080   }
00081 
00082 
00083 //! internal matrix construction; if the requested size is small enough, memory from the stack is used. otherwise memory is allocated via 'new'
00084 template<typename eT>
00085 inline
00086 void
00087 Mat<eT>::init(const u32 in_n_rows, const u32 in_n_cols)
00088   {
00089   arma_extra_debug_sigprint( arma_boost::format("in_n_rows = %d, in_n_cols = %d") % in_n_rows % in_n_cols );
00090   
00091   
00092   const u32 new_n_elem = in_n_rows * in_n_cols;
00093 
00094   if(n_elem == new_n_elem)
00095     {
00096     access::rw(n_rows) = in_n_rows;
00097     access::rw(n_cols) = in_n_cols;
00098     }
00099   else
00100     {
00101     
00102     if(n_elem > sizeof(mem_local)/sizeof(eT) )
00103       {
00104       delete [] mem;
00105       }
00106     
00107     if(new_n_elem <= sizeof(mem_local)/sizeof(eT) )
00108       {
00109       access::rw(mem) = mem_local;
00110       }
00111     else
00112       {
00113       access::rw(mem) = new(std::nothrow) eT[new_n_elem];
00114       arma_check( (mem == 0), "Mat::init(): out of memory" );
00115       }
00116     
00117     access::rw(n_elem) = new_n_elem;
00118 
00119     if(new_n_elem == 0)
00120       {
00121       access::rw(n_rows) = 0;
00122       access::rw(n_cols) = 0;
00123       }
00124     else
00125       {
00126       access::rw(n_rows) = in_n_rows;
00127       access::rw(n_cols) = in_n_cols;
00128       }
00129     
00130     }
00131   
00132   }
00133 
00134 
00135 //! create the matrix from a textual description
00136 template<typename eT>
00137 inline
00138 Mat<eT>::Mat(const char* text)
00139   : n_rows(0)
00140   , n_cols(0)
00141   , n_elem(0)
00142   //, mem(0)
00143   , mem(mem)
00144   {
00145   arma_extra_debug_sigprint_this(this);
00146   
00147   init( std::string(text) );
00148   }
00149   
00150   
00151   
00152 //! create the matrix from a textual description
00153 template<typename eT>
00154 inline
00155 const Mat<eT>&
00156 Mat<eT>::operator=(const char* text)
00157   {
00158   arma_extra_debug_sigprint();
00159   
00160   init( std::string(text) );
00161   return *this;
00162   }
00163   
00164   
00165 
00166 //! create the matrix from a textual description
00167 template<typename eT>
00168 inline
00169 Mat<eT>::Mat(const std::string& text)
00170   : n_rows(0)
00171   , n_cols(0)
00172   , n_elem(0)
00173   //, mem(0)
00174   , mem(mem)
00175   {
00176   arma_extra_debug_sigprint_this(this);
00177   
00178   init(text);
00179   }
00180   
00181   
00182   
00183 //! create the matrix from a textual description
00184 template<typename eT>
00185 inline
00186 const Mat<eT>&
00187 Mat<eT>::operator=(const std::string& text)
00188   {
00189   arma_extra_debug_sigprint();
00190   
00191   init(text);
00192   return *this;
00193   }
00194 
00195 
00196 
00197 //! internal function to create the matrix from a textual description
00198 template<typename eT>
00199 inline 
00200 void
00201 Mat<eT>::init(const std::string& text)
00202   {
00203   arma_extra_debug_sigprint();
00204   
00205   //
00206   // work out the size
00207   
00208   u32 t_n_rows = 0;
00209   u32 t_n_cols = 0;
00210   
00211   bool t_n_cols_found = false;
00212   
00213   std::string token;
00214   
00215   std::string::size_type line_start = 0;
00216   std::string::size_type   line_end = 0;
00217   
00218   while( line_start < text.length() )
00219     {
00220     
00221     line_end = text.find(';', line_start);
00222     
00223     if(line_end == std::string::npos)
00224       line_end = text.length()-1;
00225     
00226     std::string::size_type line_len = line_end - line_start + 1;
00227     std::stringstream line_stream( text.substr(line_start,line_len) );
00228     
00229     
00230     u32 line_n_cols = 0;
00231     while(line_stream >> token)
00232       {
00233       ++line_n_cols;
00234       }
00235     
00236     
00237     if(line_n_cols > 0)
00238       {
00239       if(t_n_cols_found == false)
00240         {
00241         t_n_cols = line_n_cols;
00242         t_n_cols_found = true;
00243         }
00244       else
00245         arma_check( (line_n_cols != t_n_cols), "Mat::init(): inconsistent number of columns in given string");
00246       
00247       ++t_n_rows;
00248       }
00249     line_start = line_end+1;
00250     
00251     }
00252     
00253   Mat<eT> &x = *this;
00254   x.set_size(t_n_rows, t_n_cols);
00255   
00256   line_start = 0;
00257   line_end = 0;
00258   
00259   u32 row = 0;
00260   
00261   while( line_start < text.length() )
00262     {
00263     
00264     line_end = text.find(';', line_start);
00265     
00266     if(line_end == std::string::npos)
00267       line_end = text.length()-1;
00268     
00269     std::string::size_type line_len = line_end - line_start + 1;
00270     std::stringstream line_stream( text.substr(line_start,line_len) );
00271     
00272 //     u32 col = 0;
00273 //     while(line_stream >> token)
00274 //       {
00275 //       x.at(row,col) = strtod(token.c_str(), 0);
00276 //       ++col;
00277 //       }
00278     
00279     u32 col = 0;
00280     eT val;
00281     while(line_stream >> val)
00282       {
00283       x.at(row,col) = val;
00284       ++col;
00285       }
00286     
00287     ++row;
00288     line_start = line_end+1;
00289     }
00290   
00291   }
00292 
00293 
00294 
00295 //! Set the matrix to be equal to the specified scalar.
00296 //! NOTE: the size of the matrix will be 1x1
00297 template<typename eT>
00298 arma_inline
00299 const Mat<eT>&
00300 Mat<eT>::operator=(const eT val)
00301   {
00302   arma_extra_debug_sigprint();
00303   
00304   init(1,1);
00305   access::rw(mem[0]) = val;
00306   return *this;
00307   }
00308 
00309 
00310 
00311 //! In-place addition of a scalar to all elements of the matrix
00312 template<typename eT>
00313 arma_inline
00314 const Mat<eT>&
00315 Mat<eT>::operator+=(const eT val)
00316   {
00317   arma_extra_debug_sigprint();
00318   
00319   for(u32 i=0; i<n_elem; ++i)
00320     {
00321     access::rw(mem[i]) += val;
00322     }
00323   
00324   return *this;
00325   }
00326 
00327 
00328 
00329 //! In-place subtraction of a scalar from all elements of the matrix
00330 template<typename eT>
00331 arma_inline
00332 const Mat<eT>&
00333 Mat<eT>::operator-=(const eT val)
00334   {
00335   arma_extra_debug_sigprint();
00336   
00337   for(u32 i=0; i<n_elem; ++i)
00338     {
00339     access::rw(mem[i]) -= val;
00340     }
00341       
00342   return *this;
00343   }
00344 
00345 
00346 
00347 //! In-place multiplication of all elements of the matrix with a scalar
00348 template<typename eT>
00349 arma_inline
00350 const Mat<eT>&
00351 Mat<eT>::operator*=(const eT val)
00352   {
00353   arma_extra_debug_sigprint();
00354   
00355   for(u32 i=0; i<n_elem; ++i)
00356     {
00357     access::rw(mem[i]) *= val;
00358     }
00359   
00360   return *this;
00361   }
00362 
00363 
00364 
00365 //! In-place division of all elements of the matrix with a scalar
00366 template<typename eT>
00367 arma_inline
00368 const Mat<eT>&
00369 Mat<eT>::operator/=(const eT val)
00370   {
00371   arma_extra_debug_sigprint();
00372   
00373   for(u32 i=0; i<n_elem; ++i)
00374     {
00375     access::rw(mem[i]) /= val;
00376     }
00377   
00378   return *this;
00379   }
00380 
00381 
00382 
00383 //! construct a matrix from a given matrix
00384 template<typename eT>
00385 inline
00386 Mat<eT>::Mat(const Mat<eT> &in_mat)
00387   : n_rows(0)
00388   , n_cols(0)
00389   , n_elem(0)
00390   //, mem(0)
00391   , mem(mem)
00392   {
00393   arma_extra_debug_sigprint(arma_boost::format("this = %x   in_mat = %x") % this % &in_mat);
00394   
00395   init(in_mat);
00396   }
00397 
00398 
00399 
00400 //! construct a matrix from a given matrix
00401 template<typename eT>
00402 inline
00403 const Mat<eT>&
00404 Mat<eT>::operator=(const Mat<eT>& x)
00405   {
00406   arma_extra_debug_sigprint();
00407   
00408   init(x);
00409   return *this;
00410   }
00411 
00412 
00413 
00414 //! construct a matrix from a given matrix
00415 template<typename eT>
00416 inline
00417 void
00418 Mat<eT>::init(const Mat<eT> &x)
00419   {
00420   arma_extra_debug_sigprint();
00421   
00422   if(this != &x)
00423     {
00424     init(x.n_rows, x.n_cols);
00425     syslib::copy_elem( memptr(), x.mem, n_elem );
00426     }
00427   }
00428 
00429 
00430 
00431 //! construct a matrix from a given auxillary array of eTs
00432 template<typename eT>
00433 inline
00434 Mat<eT>::Mat(const eT* aux_mem, const u32 aux_n_rows, const u32 aux_n_cols)
00435   : n_rows(0)
00436   , n_cols(0)
00437   , n_elem(0)
00438   //, mem(0)
00439   , mem(mem)
00440   {
00441   arma_extra_debug_sigprint_this(this);
00442   
00443   init(aux_n_rows, aux_n_cols);
00444   syslib::copy_elem( memptr(), aux_mem, n_elem );
00445   }
00446 
00447 
00448 
00449 //! in-place matrix addition
00450 template<typename eT>
00451 inline
00452 const Mat<eT>&
00453 Mat<eT>::operator+=(const Mat<eT>& m)
00454   {
00455   arma_extra_debug_sigprint();
00456   
00457   glue_plus::apply_inplace(*this, m);
00458   return *this;
00459   }
00460 
00461 
00462 
00463 //! in-place matrix subtraction
00464 template<typename eT>
00465 inline
00466 const Mat<eT>&
00467 Mat<eT>::operator-=(const Mat<eT>& m)
00468   {
00469   arma_extra_debug_sigprint();
00470   
00471   glue_minus::apply_inplace(*this, m);
00472   return *this;
00473   }
00474 
00475 
00476 
00477 //! in-place matrix multiplication
00478 template<typename eT>
00479 inline
00480 const Mat<eT>&
00481 Mat<eT>::operator*=(const Mat<eT>& m)
00482   {
00483   arma_extra_debug_sigprint();
00484   
00485   glue_times::apply_inplace(*this, m);
00486   return *this;
00487   }
00488 
00489 
00490 
00491 //! in-place element-wise matrix multiplication
00492 template<typename eT>
00493 inline
00494 const Mat<eT>&
00495 Mat<eT>::operator%=(const Mat<eT>& m)
00496   {
00497   arma_extra_debug_sigprint();
00498   
00499   glue_schur::apply_inplace(*this, m);
00500   return *this;
00501   }
00502 
00503 
00504 
00505 //! in-place element-wise matrix division
00506 template<typename eT>
00507 inline
00508 const Mat<eT>&
00509 Mat<eT>::operator/=(const Mat<eT>& m)
00510   {
00511   arma_extra_debug_sigprint();
00512   
00513   glue_div::apply_inplace(*this, m);
00514   return *this;
00515   }
00516 
00517 
00518 
00519 //! for constructing a complex matrix out of two non-complex matrices
00520 template<typename eT>
00521 template<typename T1, typename T2>
00522 inline
00523 Mat<eT>::Mat
00524   (
00525   const Base<typename Mat<eT>::pod_type,T1>& A,
00526   const Base<typename Mat<eT>::pod_type,T2>& B
00527   )
00528   : n_rows(0)
00529   , n_cols(0)
00530   , n_elem(0)
00531   //, mem(0)
00532   , mem(mem)
00533   {
00534   arma_extra_debug_sigprint_this(this);
00535   
00536   arma_type_check< is_complex<eT>::value == false >::apply();   //!< compile-time abort if eT isn't std::complex
00537   
00538   typedef typename T1::elem_type T;
00539   arma_type_check< is_complex<T>::value == true >::apply();   //!< compile-time abort if T is std::complex
00540   
00541   isnt_same_type<std::complex<T>, eT>::check();   //!< compile-time abort if types are not compatible
00542   
00543   const unwrap<T1> tmp_A(A.get_ref());
00544   const unwrap<T2> tmp_B(B.get_ref());
00545   
00546   const Mat<T>& X = tmp_A.M;
00547   const Mat<T>& Y = tmp_B.M;
00548   
00549   arma_assert_same_size(X, Y, "Mat()");
00550   
00551   init(X.n_rows, Y.n_cols);
00552   
00553   const T* X_mem = X.mem;
00554   const T* Y_mem = Y.mem;
00555   
00556   for(u32 i=0; i<n_elem; ++i)
00557     {
00558     access::rw(mem[i]) = std::complex<T>(X_mem[i], Y_mem[i]);
00559     }
00560   }
00561 
00562 
00563 
00564 //! construct a matrix from subview (e.g. construct a matrix from a delayed submatrix operation)
00565 template<typename eT>
00566 inline
00567 Mat<eT>::Mat(const subview<eT>& X)
00568   : n_rows(0)
00569   , n_cols(0)
00570   , n_elem(0)
00571   //, mem(0)
00572   , mem(mem)
00573   {
00574   arma_extra_debug_sigprint_this(this);
00575   
00576   this->operator=(X);
00577   }
00578 
00579 
00580 
00581 //! construct a matrix from subview (e.g. construct a matrix from a delayed submatrix operation)
00582 template<typename eT>
00583 inline
00584 const Mat<eT>&
00585 Mat<eT>::operator=(const subview<eT>& X)
00586   {
00587   arma_extra_debug_sigprint();
00588   
00589   subview<eT>::extract(*this, X);
00590   return *this;
00591   }
00592 
00593 
00594 //! in-place matrix addition (using a submatrix on the right-hand-side)
00595 template<typename eT>
00596 inline
00597 const Mat<eT>&
00598 Mat<eT>::operator+=(const subview<eT>& X)
00599   {
00600   arma_extra_debug_sigprint();
00601   
00602   subview<eT>::plus_inplace(*this, X);
00603   return *this;
00604   }
00605 
00606 
00607 //! in-place matrix subtraction (using a submatrix on the right-hand-side)
00608 template<typename eT>
00609 inline
00610 const Mat<eT>&
00611 Mat<eT>::operator-=(const subview<eT>& X)
00612   {
00613   arma_extra_debug_sigprint();
00614   
00615   subview<eT>::minus_inplace(*this, X);
00616   return *this;
00617   }
00618 
00619 
00620 
00621 //! in-place matrix mutiplication (using a submatrix on the right-hand-side)
00622 template<typename eT>
00623 inline
00624 const Mat<eT>&
00625 Mat<eT>::operator*=(const subview<eT>& X)
00626   {
00627   arma_extra_debug_sigprint();
00628   
00629   subview<eT>::times_inplace(*this, X);
00630   return *this;
00631   }
00632 
00633 
00634 
00635 //! in-place element-wise matrix mutiplication (using a submatrix on the right-hand-side)
00636 template<typename eT>
00637 inline
00638 const Mat<eT>&
00639 Mat<eT>::operator%=(const subview<eT>& X)
00640   {
00641   arma_extra_debug_sigprint();
00642   
00643   subview<eT>::schur_inplace(*this, X);
00644   return *this;
00645   }
00646 
00647 
00648 
00649 //! in-place element-wise matrix division (using a submatrix on the right-hand-side)
00650 template<typename eT>
00651 inline
00652 const Mat<eT>&
00653 Mat<eT>::operator/=(const subview<eT>& X)
00654   {
00655   arma_extra_debug_sigprint();
00656   
00657   subview<eT>::div_inplace(*this, X);
00658   return *this;
00659   }
00660 
00661 
00662 
00663 //! construct a matrix from diagview (e.g. construct a matrix from a delayed diag operation)
00664 template<typename eT>
00665 inline
00666 Mat<eT>::Mat(const diagview<eT>& X)
00667   : n_rows(0)
00668   , n_cols(0)
00669   , n_elem(0)
00670   //, mem(0)
00671   , mem(mem)
00672   {
00673   arma_extra_debug_sigprint_this(this);
00674   
00675   this->operator=(X);
00676   }
00677 
00678 
00679 
00680 //! construct a matrix from diagview (e.g. construct a matrix from a delayed diag operation)
00681 template<typename eT>
00682 inline
00683 const Mat<eT>&
00684 Mat<eT>::operator=(const diagview<eT>& X)
00685   {
00686   arma_extra_debug_sigprint();
00687   
00688   diagview<eT>::extract(*this, X);
00689   return *this;
00690   }
00691 
00692 
00693 
00694 //! creation of subview (row vector)
00695 template<typename eT>
00696 arma_inline
00697 subview_row<eT>
00698 Mat<eT>::row(const u32 row_num)
00699   {
00700   arma_extra_debug_sigprint();
00701   
00702   arma_debug_check( row_num >= n_rows, "Mat::row(): row out of bounds" );
00703   
00704   return subview_row<eT>(*this, row_num);
00705   }
00706 
00707 
00708 
00709 //! creation of subview (row vector)
00710 template<typename eT>
00711 arma_inline
00712 const subview_row<eT>
00713 Mat<eT>::row(const u32 row_num) const
00714   {
00715   arma_extra_debug_sigprint();
00716   
00717   arma_debug_check( row_num >= n_rows, "Mat::row(): row out of bounds" );
00718   
00719   return subview_row<eT>(*this, row_num);
00720   }
00721 
00722 
00723 
00724 //! creation of subview (column vector)
00725 template<typename eT>
00726 arma_inline
00727 subview_col<eT>
00728 Mat<eT>::col(const u32 col_num)
00729   {
00730   arma_extra_debug_sigprint();
00731   
00732   arma_debug_check( col_num >= n_cols, "Mat::col(): out of bounds");
00733   
00734   return subview_col<eT>(*this, col_num);
00735   }
00736 
00737 
00738 
00739 //! creation of subview (column vector)
00740 template<typename eT>
00741 arma_inline
00742 const subview_col<eT>
00743 Mat<eT>::col(const u32 col_num) const
00744   {
00745   arma_extra_debug_sigprint();
00746   
00747   arma_debug_check( col_num >= n_cols, "Mat::col(): out of bounds");
00748   
00749   return subview_col<eT>(*this, col_num);
00750   }
00751 
00752 
00753 
00754 //! creation of subview (submatrix comprised of specified row vectors)
00755 template<typename eT>
00756 arma_inline
00757 subview<eT>
00758 Mat<eT>::rows(const u32 in_row1, const u32 in_row2)
00759   {
00760   arma_extra_debug_sigprint();
00761   
00762   arma_debug_check
00763     (
00764     (in_row1 > in_row2) || (in_row2 >= n_rows),
00765     "Mat::rows(): indices out of bounds or incorrectly used"
00766     );
00767   
00768   return subview<eT>(*this, in_row1, 0, in_row2, n_cols-1);
00769   }
00770 
00771 
00772 
00773 //! creation of subview (submatrix comprised of specified row vectors)
00774 template<typename eT>
00775 arma_inline
00776 const subview<eT>
00777 Mat<eT>::rows(const u32 in_row1, const u32 in_row2) const
00778   {
00779   arma_extra_debug_sigprint();
00780   
00781   arma_debug_check
00782     (
00783     (in_row1 > in_row2) || (in_row2 >= n_rows),
00784     "Mat::rows(): indices out of bounds or incorrectly used"
00785     );
00786   
00787   return subview<eT>(*this, in_row1, 0, in_row2, n_cols-1);
00788   }
00789 
00790 
00791 
00792 //! creation of subview (submatrix comprised of specified column vectors)
00793 template<typename eT>
00794 arma_inline
00795 subview<eT>
00796 Mat<eT>::cols(const u32 in_col1, const u32 in_col2)
00797   {
00798   arma_extra_debug_sigprint();
00799   
00800   arma_debug_check
00801     (
00802     (in_col1 > in_col2) || (in_col2 >= n_cols),
00803     "Mat::cols(): indices out of bounds or incorrectly used"
00804     );
00805   
00806   return subview<eT>(*this, 0, in_col1, n_rows-1, in_col2);
00807   }
00808 
00809 
00810 
00811 //! creation of subview (submatrix comprised of specified column vectors)
00812 template<typename eT>
00813 arma_inline
00814 const subview<eT>
00815 Mat<eT>::cols(const u32 in_col1, const u32 in_col2) const
00816   {
00817   arma_extra_debug_sigprint();
00818   
00819   arma_debug_check
00820     (
00821     (in_col1 > in_col2) || (in_col2 >= n_cols),
00822     "Mat::cols(): indices out of bounds or incorrectly used"
00823     );
00824   
00825   return subview<eT>(*this, 0, in_col1, n_rows-1, in_col2);
00826   }
00827 
00828 
00829 
00830 //! creation of subview (submatrix)
00831 template<typename eT>
00832 arma_inline
00833 subview<eT>
00834 Mat<eT>::submat(const u32 in_row1, const u32 in_col1, const u32 in_row2, const u32 in_col2)
00835   {
00836   arma_extra_debug_sigprint();
00837   
00838   arma_debug_check
00839     (
00840     (in_row1 > in_row2) || (in_col1 >  in_col2) || (in_row2 >= n_rows) || (in_col2 >= n_cols),
00841     "Mat::submat(): indices out of bounds or incorrectly used"
00842     );
00843   
00844   return subview<eT>(*this, in_row1, in_col1, in_row2, in_col2);
00845   }
00846 
00847 
00848 
00849 //! creation of subview (generic submatrix)
00850 template<typename eT>
00851 arma_inline
00852 const subview<eT>
00853 Mat<eT>::submat(const u32 in_row1, const u32 in_col1, const u32 in_row2, const u32 in_col2) const
00854   {
00855   arma_extra_debug_sigprint();
00856   
00857   arma_debug_check
00858     (
00859     (in_row1 > in_row2) || (in_col1 >  in_col2) || (in_row2 >= n_rows) || (in_col2 >= n_cols),
00860     "Mat::submat(): indices out of bounds or incorrectly used"
00861     );
00862     
00863   return subview<eT>(*this, in_row1, in_col1, in_row2, in_col2);
00864   }
00865 
00866 
00867 
00868 //! creation of diagview (diagonal)
00869 template<typename eT>
00870 arma_inline
00871 diagview<eT>
00872 Mat<eT>::diag(const s32 in_id)
00873   {
00874   arma_extra_debug_sigprint();
00875   
00876   const u32 row_offset = (in_id < 0) ? -in_id : 0;
00877   const u32 col_offset = (in_id > 0) ?  in_id : 0;
00878   
00879   arma_debug_check
00880     (
00881     (row_offset >= n_rows) || (col_offset >= n_cols),
00882     "Mat::diag(): out of bounds"
00883     );
00884   
00885   const u32 len = (std::min)(n_rows - row_offset, n_cols - col_offset);
00886   
00887   return diagview<eT>(*this, row_offset, col_offset, len);
00888   }
00889 
00890 
00891 
00892 //! creation of diagview (diagonal)
00893 template<typename eT>
00894 arma_inline
00895 const diagview<eT>
00896 Mat<eT>::diag(const s32 in_id) const
00897   {
00898   arma_extra_debug_sigprint();
00899   
00900   const u32 row_offset = (in_id < 0) ? -in_id : 0;
00901   const u32 col_offset = (in_id > 0) ?  in_id : 0;
00902   
00903   arma_debug_check
00904     (
00905     (row_offset >= n_rows) || (col_offset >= n_cols),
00906     "Mat::diag(): out of bounds"
00907     );
00908   
00909   
00910   const u32 len = (std::min)(n_rows - row_offset, n_cols - col_offset);
00911   
00912   return diagview<eT>(*this, row_offset, col_offset, len);
00913   }
00914 
00915 
00916 
00917 template<typename eT>
00918 inline
00919 void
00920 Mat<eT>::swap_rows(const u32 in_row1, const u32 in_row2)
00921   {
00922   arma_extra_debug_sigprint();
00923   
00924   arma_debug_check
00925     (
00926     (in_row1 >= n_rows) || (in_row2 >= n_rows),
00927     "Mat::swap_rows(): out of bounds"
00928     );
00929   
00930   for(u32 col=0; col<n_cols; ++col)
00931     {
00932     const u32 offset = col*n_rows;
00933     const u32 pos1   = in_row1 + offset;
00934     const u32 pos2   = in_row2 + offset;
00935     
00936     const eT tmp          = mem[pos1];
00937     access::rw(mem[pos1]) = mem[pos2];
00938     access::rw(mem[pos2]) = tmp;
00939     }
00940   
00941   }
00942 
00943 
00944 
00945 template<typename eT>
00946 inline
00947 void
00948 Mat<eT>::swap_cols(const u32 in_col1, const u32 in_col2)
00949   {
00950   arma_extra_debug_sigprint();
00951   
00952   arma_debug_check
00953     (
00954     (in_col1 >= n_cols) || (in_col2 >= n_cols),
00955     "Mat::swap_cols(): out of bounds"
00956     );
00957   
00958   eT* ptr1 = colptr(in_col1);
00959   eT* ptr2 = colptr(in_col2);
00960   
00961   for(u32 row=0; row<n_rows; ++row)
00962     {
00963     const eT tmp = ptr1[row];
00964     ptr1[row]    = ptr2[row];
00965     ptr2[row]    = tmp;
00966     }
00967   
00968   }
00969 
00970 
00971 
00972 //! create a matrix from Op, i.e. run the previously delayed unary operations
00973 template<typename eT>
00974 template<typename T1, typename op_type>
00975 inline
00976 Mat<eT>::Mat(const Op<T1, op_type>& X)
00977   : n_rows(0)
00978   , n_cols(0)
00979   , n_elem(0)
00980   //, mem(0)
00981   , mem(mem)
00982   {
00983   arma_extra_debug_sigprint_this(this);
00984 
00985   isnt_same_type<eT, typename T1::elem_type>::check();
00986   
00987   op_type::apply(*this, X);
00988   }
00989 
00990 
00991 
00992 //! create a matrix from Op, i.e. run the previously delayed unary operations
00993 template<typename eT>
00994 template<typename T1, typename op_type>
00995 inline
00996 const Mat<eT>&
00997 Mat<eT>::operator=(const Op<T1, op_type>& X)
00998   {
00999   arma_extra_debug_sigprint();
01000 
01001   isnt_same_type<eT, typename T1::elem_type>::check();
01002   
01003   op_type::apply(*this, X);
01004   
01005   return *this;
01006   }
01007 
01008 
01009 
01010 //! in-place matrix addition, with the right-hand-side operand having delayed operations
01011 template<typename eT>
01012 template<typename T1, typename op_type>
01013 inline
01014 const Mat<eT>&
01015 Mat<eT>::operator+=(const Op<T1, op_type>& X)
01016   {
01017   arma_extra_debug_sigprint();
01018   
01019   isnt_same_type<eT, typename T1::elem_type>::check();
01020   
01021   glue_plus::apply_inplace(*this, X);
01022   
01023   return *this;
01024   }
01025 
01026 
01027 
01028 //! in-place matrix subtraction, with the right-hand-side operand having delayed operations
01029 template<typename eT>
01030 template<typename T1, typename op_type>
01031 inline
01032 const Mat<eT>&
01033 Mat<eT>::operator-=(const Op<T1, op_type>& X)
01034   {
01035   arma_extra_debug_sigprint();
01036   
01037   isnt_same_type<eT, typename T1::elem_type>::check();
01038   
01039   glue_minus::apply_inplace(*this, X);
01040   
01041   return *this;
01042   }
01043 
01044 
01045 
01046 //! in-place matrix multiplication, with the right-hand-side operand having delayed operations
01047 template<typename eT>
01048 template<typename T1, typename op_type>
01049 inline
01050 const Mat<eT>&
01051 Mat<eT>::operator*=(const Op<T1, op_type>& X)
01052   {
01053   arma_extra_debug_sigprint();
01054   
01055   isnt_same_type<eT, typename T1::elem_type>::check();
01056   
01057   glue_times::apply_inplace(*this, X);
01058   
01059   return *this;
01060   }
01061 
01062 
01063 
01064 //! in-place matrix element-wise multiplication, with the right-hand-side operand having delayed operations
01065 template<typename eT>
01066 template<typename T1, typename op_type>
01067 inline
01068 const Mat<eT>&
01069 Mat<eT>::operator%=(const Op<T1, op_type>& X)
01070   {
01071   arma_extra_debug_sigprint();
01072   
01073   isnt_same_type<eT, typename T1::elem_type>::check();
01074   glue_schur::apply_inplace(*this, X);
01075   
01076   return *this;
01077   }
01078 
01079 
01080 
01081 //! in-place matrix element-wise division, with the right-hand-side operand having delayed operations
01082 template<typename eT>
01083 template<typename T1, typename op_type>
01084 inline
01085 const Mat<eT>&
01086 Mat<eT>::operator/=(const Op<T1, op_type>& X)
01087   {
01088   arma_extra_debug_sigprint();
01089   
01090   isnt_same_type<eT, typename T1::elem_type>::check();
01091   glue_div::apply_inplace(*this, X);
01092   
01093   return *this;
01094   }
01095 
01096 
01097 
01098 //! create a matrix from Glue, i.e. run the previously delayed binary operations
01099 template<typename eT>
01100 template<typename T1, typename T2, typename glue_type>
01101 inline
01102 Mat<eT>::Mat(const Glue<T1, T2, glue_type>& X)
01103   : n_rows(0)
01104   , n_cols(0)
01105   , n_elem(0)
01106   //, mem(0)
01107   , mem(mem)
01108   {
01109   arma_extra_debug_sigprint_this(this);
01110   this->operator=(X);
01111   }
01112 
01113 
01114 
01115 //! create a matrix from Glue, i.e. run the previously delayed binary operations
01116 template<typename eT>
01117 template<typename T1, typename T2, typename glue_type>
01118 inline
01119 const Mat<eT>&
01120 Mat<eT>::operator=(const Glue<T1, T2, glue_type>& X)
01121   {
01122   arma_extra_debug_sigprint();
01123   
01124   // TODO:
01125   // it may be simpler to pass the two objects (currently wrapped in X)
01126   // directly to the apply function.
01127   // (many adjustments throughout the source code will be required)
01128   
01129   isnt_same_type<eT, typename T1::elem_type>::check();
01130   isnt_same_type<eT, typename T2::elem_type>::check();
01131   
01132   glue_type::apply(*this, X);
01133   
01134   return *this;
01135   }
01136 
01137 
01138 //! in-place matrix addition, with the right-hand-side operands having delayed operations
01139 template<typename eT>
01140 template<typename T1, typename T2, typename glue_type>
01141 inline
01142 const Mat<eT>&
01143 Mat<eT>::operator+=(const Glue<T1, T2, glue_type>& X)
01144   {
01145   arma_extra_debug_sigprint();
01146   
01147   isnt_same_type<eT, typename T1::elem_type>::check();
01148   isnt_same_type<eT, typename T2::elem_type>::check();
01149   
01150   glue_plus::apply_inplace(*this, X);
01151   
01152   return *this;
01153   }
01154 
01155 
01156 
01157 //! in-place matrix subtraction, with the right-hand-side operands having delayed operations
01158 template<typename eT>
01159 template<typename T1, typename T2, typename glue_type>
01160 inline
01161 const Mat<eT>&
01162 Mat<eT>::operator-=(const Glue<T1, T2, glue_type>& X)
01163   {
01164   arma_extra_debug_sigprint();
01165   
01166   isnt_same_type<eT, typename T1::elem_type>::check();
01167   isnt_same_type<eT, typename T2::elem_type>::check();
01168   
01169   glue_minus::apply_inplace(*this, X);
01170   
01171   return *this;
01172   }
01173 
01174 
01175 
01176 //! in-place matrix multiplications, with the right-hand-side operands having delayed operations
01177 template<typename eT>
01178 template<typename T1, typename T2, typename glue_type>
01179 inline
01180 const Mat<eT>&
01181 Mat<eT>::operator*=(const Glue<T1, T2, glue_type>& X)
01182   {
01183   arma_extra_debug_sigprint();
01184   
01185   isnt_same_type<eT, typename T1::elem_type>::check();
01186   isnt_same_type<eT, typename T2::elem_type>::check();
01187   
01188   glue_times::apply_inplace(*this, X);
01189   return *this;
01190   }
01191 
01192 
01193 
01194 //! in-place matrix element-wise multiplication, with the right-hand-side operands having delayed operations
01195 template<typename eT>
01196 template<typename T1, typename T2, typename glue_type>
01197 inline
01198 const Mat<eT>&
01199 Mat<eT>::operator%=(const Glue<T1, T2, glue_type>& X)
01200   {
01201   arma_extra_debug_sigprint();
01202   
01203   isnt_same_type<eT, typename T1::elem_type>::check();
01204   isnt_same_type<eT, typename T2::elem_type>::check();
01205   
01206   glue_schur::apply_inplace(*this, X);
01207   return *this;
01208   }
01209 
01210 
01211 
01212 //! in-place matrix element-wise division, with the right-hand-side operands having delayed operations
01213 template<typename eT>
01214 template<typename T1, typename T2, typename glue_type>
01215 inline
01216 const Mat<eT>&
01217 Mat<eT>::operator/=(const Glue<T1, T2, glue_type>& X)
01218   {
01219   arma_extra_debug_sigprint();
01220   
01221   isnt_same_type<eT, typename T1::elem_type>::check();
01222   isnt_same_type<eT, typename T2::elem_type>::check();
01223   
01224   glue_div::apply_inplace(*this, X);
01225   return *this;
01226   }
01227 
01228 
01229 
01230 //! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
01231 template<typename eT>
01232 arma_inline
01233 eT&
01234 Mat<eT>::operator() (const u32 i)
01235   {
01236   arma_debug_check( (i >= n_elem), "Mat::operator(): index out of bounds");
01237   return access::rw(mem[i]);
01238   }
01239 
01240 
01241 
01242 //! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
01243 template<typename eT>
01244 arma_inline
01245 eT
01246 Mat<eT>::operator() (const u32 i) const
01247   {
01248   arma_debug_check( (i >= n_elem), "Mat::operator(): index out of bounds");
01249   return mem[i];
01250   }
01251 
01252 
01253 //! linear element accessor (treats the matrix as a vector); no bounds check.  
01254 template<typename eT>
01255 arma_inline
01256 eT&
01257 Mat<eT>::operator[] (const u32 i)
01258   {
01259   return access::rw(mem[i]);
01260   }
01261 
01262 
01263 
01264 //! linear element accessor (treats the matrix as a vector); no bounds check
01265 template<typename eT>
01266 arma_inline
01267 eT
01268 Mat<eT>::operator[] (const u32 i) const
01269   {
01270   return mem[i];
01271   }
01272 
01273 
01274 
01275 //! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
01276 template<typename eT>
01277 arma_inline
01278 eT&
01279 Mat<eT>::operator() (const u32 in_row, const u32 in_col)
01280   {
01281   arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "Mat::operator(): index out of bounds");
01282   return access::rw(mem[in_row + in_col*n_rows]);
01283   }
01284 
01285 
01286 
01287 //! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
01288 template<typename eT>
01289 arma_inline
01290 eT
01291 Mat<eT>::operator() (const u32 in_row, const u32 in_col) const
01292   {
01293   arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "Mat::operator(): index out of bounds");
01294   return mem[in_row + in_col*n_rows];
01295   }
01296 
01297 
01298 
01299 //! element accessor; no bounds check
01300 template<typename eT>
01301 arma_inline
01302 eT&
01303 Mat<eT>::at(const u32 in_row, const u32 in_col)
01304   {
01305   return access::rw( mem[in_row + in_col*n_rows] );
01306   }
01307 
01308 
01309 
01310 //! element accessor; no bounds check
01311 template<typename eT>
01312 arma_inline
01313 eT
01314 Mat<eT>::at(const u32 in_row, const u32 in_col) const
01315   {
01316   return mem[in_row + in_col*n_rows];
01317   }
01318 
01319 
01320 
01321 //! prefix ++
01322 template<typename eT>
01323 arma_inline
01324 const Mat<eT>&
01325 Mat<eT>::operator++()
01326   {
01327   Mat_aux::prefix_pp(*this);
01328   return *this;
01329   }
01330 
01331 
01332 
01333 //! postfix ++  (must not return the object by reference)
01334 template<typename eT>
01335 arma_inline
01336 void
01337 Mat<eT>::operator++(int)
01338   {
01339   Mat_aux::postfix_pp(*this);
01340   }
01341 
01342 
01343 
01344 //! prefix --
01345 template<typename eT>
01346 arma_inline
01347 const Mat<eT>&
01348 Mat<eT>::operator--()
01349   {
01350   Mat_aux::prefix_mm(*this);
01351   return *this;
01352   }
01353 
01354 
01355 
01356 //! postfix --  (must not return the object by reference)
01357 template<typename eT>
01358 arma_inline
01359 void
01360 Mat<eT>::operator--(int)
01361   {
01362   Mat_aux::postfix_mm(*this);
01363   }
01364 
01365 
01366 
01367 //! returns true if the object can be interpreted as a column or row vector
01368 template<typename eT>
01369 arma_inline
01370 bool
01371 Mat<eT>::is_vec() const
01372   {
01373   return ( (n_rows == 1) || (n_cols == 1) );
01374   }
01375 
01376 
01377 
01378 //! returns true if the object has the same number of non-zero rows and columnns
01379 template<typename eT>
01380 arma_inline
01381 bool
01382 Mat<eT>::is_square() const
01383   {
01384   return ( (n_rows == n_cols) && (n_elem > 0) );
01385   }
01386 
01387 
01388 
01389 //! returns true if all of the elements are finite
01390 template<typename eT>
01391 arma_inline
01392 bool
01393 Mat<eT>::is_finite() const
01394   {
01395   for(u32 i=0; i<n_elem; ++i)
01396     {
01397     if(arma_isfinite(mem[i]) == false)
01398       {
01399       return false;
01400       }
01401     }
01402 
01403   return true;
01404   }
01405 
01406 
01407 
01408 //! returns a pointer to array of eTs for a specified column; no bounds check
01409 template<typename eT>
01410 arma_inline
01411 eT*
01412 Mat<eT>::colptr(const u32 in_col)
01413   {
01414   return & access::rw(mem[in_col*n_rows]);
01415   }
01416 
01417 
01418 
01419 //! returns a pointer to array of eTs for a specified column; no bounds check
01420 template<typename eT>
01421 arma_inline
01422 const eT*
01423 Mat<eT>::colptr(const u32 in_col) const
01424   {
01425   return & mem[in_col*n_rows];
01426   }
01427 
01428 
01429 
01430 //! returns a pointer to array of eTs used by the matrix
01431 template<typename eT>
01432 arma_inline
01433 eT*
01434 Mat<eT>::memptr()
01435   {
01436   return const_cast<eT*>(mem);
01437   }
01438 
01439 
01440 
01441 //! returns a pointer to array of eTs used by the matrix
01442 template<typename eT>
01443 arma_inline
01444 const eT*
01445 Mat<eT>::memptr() const
01446   {
01447   return mem;
01448   }
01449 
01450 
01451 
01452 //! print contents of the matrix, optionally preceding with a user specified line of text
01453 template<typename eT>
01454 inline
01455 void
01456 Mat<eT>::print(const std::string extra_text) const
01457   {
01458   arma_extra_debug_sigprint();
01459   
01460   if(extra_text.length() != 0)
01461     {
01462     cout << extra_text << '\n';
01463     }
01464   
01465   cout << *this << '\n';
01466   }
01467 
01468 
01469 
01470 //! fill the matrix with the specified value
01471 template<typename eT>
01472 inline
01473 void
01474 Mat<eT>::fill(const eT val)
01475   {
01476   arma_extra_debug_sigprint();
01477   
01478   for(u32 i=0; i<n_elem; ++i)
01479     {
01480     access::rw(mem[i]) = val;
01481     }
01482   }
01483 
01484 
01485 
01486 template<typename eT>
01487 inline
01488 void
01489 Mat<eT>::zeros()
01490   {
01491   arma_extra_debug_sigprint();
01492   
01493   fill(eT(0));
01494   }
01495 
01496 
01497 
01498 template<typename eT>
01499 inline
01500 void
01501 Mat<eT>::zeros(const u32 in_rows, const u32 in_cols)
01502   {
01503   arma_extra_debug_sigprint( arma_boost::format("in_rows = %d, in_cols = %d") % in_rows % in_cols );
01504 
01505   set_size(in_rows, in_cols);
01506   fill(eT(0));
01507   }
01508 
01509 
01510 
01511 template<typename eT>
01512 inline
01513 void
01514 Mat<eT>::reset()
01515   {
01516   arma_extra_debug_sigprint();
01517   
01518   init(0,0);
01519   }
01520 
01521 
01522 
01523 //! save the matrix to a file
01524 template<typename eT>
01525 inline
01526 void
01527 Mat<eT>::save(const std::string name, const file_type type) const
01528   {
01529   arma_extra_debug_sigprint();
01530   
01531   switch(type)
01532     {
01533     case raw_ascii:
01534       diskio::save_raw_ascii(*this, name);
01535       break;
01536     
01537     case arma_ascii:
01538       diskio::save_arma_ascii(*this, name);
01539       break;
01540     
01541     case arma_binary:
01542       diskio::save_arma_binary(*this, name);
01543       break;
01544       
01545     case pgm_binary:
01546       diskio::save_pgm_binary(*this, name);
01547       break;
01548     
01549     default:
01550       arma_stop("Mat::save(): unsupported type");
01551     }
01552   
01553   }
01554 
01555 
01556 
01557 //! load a matrix from a file
01558 template<typename eT>
01559 inline
01560 void
01561 Mat<eT>::load(const std::string name, const file_type type)
01562   {
01563   arma_extra_debug_sigprint();
01564   
01565   switch(type)
01566     {
01567     case auto_detect:
01568       diskio::load_auto_detect(*this, name);
01569       break;
01570     
01571     case raw_ascii:
01572       diskio::load_raw_ascii(*this, name);
01573       break;
01574     
01575     case arma_ascii:
01576       diskio::load_arma_ascii(*this, name);
01577       break;
01578     
01579     case arma_binary:
01580       diskio::load_arma_binary(*this, name);
01581       break;
01582       
01583     case pgm_binary:
01584       diskio::load_pgm_binary(*this, name);
01585       break;
01586     
01587     default:
01588       arma_stop("Mat::load(): unsupported type");
01589     }
01590   
01591   }
01592 
01593 
01594 
01595 //! prefix ++
01596 template<typename eT>
01597 arma_inline
01598 void
01599 Mat_aux::prefix_pp(Mat<eT>& x)
01600   {
01601         eT* memptr = x.memptr();
01602   const u32 n_elem = x.n_elem;
01603 
01604   for(u32 i=0; i<n_elem; ++i)
01605     {
01606     ++(memptr[i]);
01607     }
01608   }
01609 
01610 
01611 
01612 //! prefix ++ for complex numbers (work around for limitations of the std::complex class)
01613 template<typename T>
01614 arma_inline
01615 void
01616 Mat_aux::prefix_pp(Mat< std::complex<T> >& x)
01617   {
01618   x += T(1);
01619   }
01620 
01621 
01622 
01623 //! postfix ++
01624 template<typename eT>
01625 arma_inline
01626 void
01627 Mat_aux::postfix_pp(Mat<eT>& x)
01628   {
01629         eT* memptr = x.memptr();
01630   const u32 n_elem = x.n_elem;
01631 
01632   for(u32 i=0; i<n_elem; ++i)
01633     {
01634     (memptr[i])++;
01635     }
01636   }
01637 
01638 
01639 
01640 //! postfix ++ for complex numbers (work around for limitations of the std::complex class)
01641 template<typename T>
01642 arma_inline
01643 void
01644 Mat_aux::postfix_pp(Mat< std::complex<T> >& x)
01645   {
01646   x += T(1);
01647   }
01648 
01649 
01650 
01651 //! prefix --
01652 template<typename eT>
01653 arma_inline
01654 void
01655 Mat_aux::prefix_mm(Mat<eT>& x)
01656   {
01657         eT* memptr = x.memptr();
01658   const u32 n_elem = x.n_elem;
01659 
01660   for(u32 i=0; i<n_elem; ++i)
01661     {
01662     --(memptr[i]);
01663     }
01664   }
01665 
01666 
01667 
01668 //! prefix -- for complex numbers (work around for limitations of the std::complex class)
01669 template<typename T>
01670 arma_inline
01671 void
01672 Mat_aux::prefix_mm(Mat< std::complex<T> >& x)
01673   {
01674   x -= T(1);
01675   }
01676 
01677 
01678 
01679 //! postfix --
01680 template<typename eT>
01681 arma_inline
01682 void
01683 Mat_aux::postfix_mm(Mat<eT>& x)
01684   {
01685         eT* memptr = x.memptr();
01686   const u32 n_elem = x.n_elem;
01687 
01688   for(u32 i=0; i<n_elem; ++i)
01689     {
01690     (memptr[i])--;
01691     }
01692   }
01693 
01694 
01695 
01696 //! postfix ++ for complex numbers (work around for limitations of the std::complex class)
01697 template<typename T>
01698 arma_inline
01699 void
01700 Mat_aux::postfix_mm(Mat< std::complex<T> >& x)
01701   {
01702   x -= T(1);
01703   }
01704 
01705 
01706 
01707 //! @}