Chris@49: // Copyright (C) 2011-2013 Ryan Curtin
Chris@49: // Copyright (C) 2012-2013 Conrad Sanderson
Chris@49: // Copyright (C) 2011 Matthew Amidon
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: //! \addtogroup SpMat
Chris@49: //! @{
Chris@49: 
Chris@49: /**
Chris@49:  * Initialize a sparse matrix with size 0x0 (empty).
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat()
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(memory::acquire_chunked<eT>(1))
Chris@49:   , row_indices(memory::acquire_chunked<uword>(1))
Chris@49:   , col_ptrs(memory::acquire<uword>(2))
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49:   
Chris@49:   access::rw(values[0]) = 0;
Chris@49:   access::rw(row_indices[0]) = 0;
Chris@49:   
Chris@49:   access::rw(col_ptrs[0]) = 0; // No elements.
Chris@49:   access::rw(col_ptrs[1]) = std::numeric_limits<uword>::max();
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Clean up the memory of a sparse matrix and destruct it.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::~SpMat()
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   // If necessary, release the memory.
Chris@49:   if (values)
Chris@49:     {
Chris@49:     // values being non-NULL implies row_indices is non-NULL.
Chris@49:     memory::release(access::rw(values));
Chris@49:     memory::release(access::rw(row_indices));
Chris@49:     }
Chris@49: 
Chris@49:   // Column pointers always must be deleted.
Chris@49:   memory::release(access::rw(col_ptrs));
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Constructor with size given.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const uword in_rows, const uword in_cols)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   init(in_rows, in_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Assemble from text.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const char* text)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   init(std::string(text));
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const char* text)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(std::string(text));
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const std::string& text)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(text);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const std::string& text)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(text);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const SpMat<eT>& x)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   init(x);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Insert a large number of values at once.
Chris@49: //! locations.row[0] should be row indices, locations.row[1] should be column indices,
Chris@49: //! and values should be the corresponding values.
Chris@49: //! If sort_locations is false, then it is assumed that the locations and values
Chris@49: //! are already sorted in column-major ordering.
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const Base<uword,T1>& locations_expr, const Base<eT,T2>& vals_expr, const bool sort_locations)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49:   
Chris@49:   const unwrap<T1>         locs_tmp( locations_expr.get_ref() );
Chris@49:   const Mat<uword>& locs = locs_tmp.M;
Chris@49:   
Chris@49:   const unwrap<T2> vals_tmp( vals_expr.get_ref() );
Chris@49:   const Mat<eT>& vals = vals_tmp.M;
Chris@49:   
Chris@49:   arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'values' object is not a vector" );
Chris@49:   
Chris@49:   arma_debug_check((locs.n_cols != vals.n_elem), "SpMat::SpMat(): number of locations is different than number of values");
Chris@49:   
Chris@49:   // If there are no elements in the list, max() will fail.
Chris@49:   if (locs.n_cols == 0)
Chris@49:     {
Chris@49:     init(0, 0);
Chris@49:     return;
Chris@49:     }
Chris@49: 
Chris@49:   arma_debug_check((locs.n_rows != 2), "SpMat::SpMat(): locations matrix must have two rows");
Chris@49: 
Chris@49:   // Automatically determine size (and check if it's sorted).
Chris@49:   uvec bounds = arma::max(locs, 1);
Chris@49:   init(bounds[0] + 1, bounds[1] + 1);
Chris@49: 
Chris@49:   // Resize to correct number of elements.
Chris@49:   mem_resize(vals.n_elem);
Chris@49: 
Chris@49:   // Reset column pointers to zero.
Chris@49:   arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1);
Chris@49: 
Chris@49:   bool actually_sorted = true;
Chris@49:   if(sort_locations == true)
Chris@49:     {
Chris@49:     // sort_index() uses std::sort() which may use quicksort... so we better
Chris@49:     // make sure it's not already sorted before taking an O(N^2) sort penalty.
Chris@49:     for (uword i = 1; i < locs.n_cols; ++i)
Chris@49:       {
Chris@49:       if ((locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) <= locs.at(0, i - 1)))
Chris@49:         {
Chris@49:         actually_sorted = false;
Chris@49:         break;
Chris@49:         }
Chris@49:       }
Chris@49:     
Chris@49:     if(actually_sorted == false)
Chris@49:       {
Chris@49:       // This may not be the fastest possible implementation but it maximizes code reuse.
Chris@49:       Col<uword> abslocs(locs.n_cols);
Chris@49:       
Chris@49:       for (uword i = 0; i < locs.n_cols; ++i)
Chris@49:         {
Chris@49:         abslocs[i] = locs.at(1, i) * n_rows + locs.at(0, i);
Chris@49:         }
Chris@49:       
Chris@49:       // Now we will sort with sort_index().
Chris@49:       uvec sorted_indices = sort_index(abslocs); // Ascending sort.
Chris@49:       
Chris@49:       // Now we add the elements in this sorted order.
Chris@49:       for (uword i = 0; i < sorted_indices.n_elem; ++i)
Chris@49:         {
Chris@49:         arma_debug_check((locs.at(0, sorted_indices[i]) >= n_rows), "SpMat::SpMat(): invalid row index");
Chris@49:         arma_debug_check((locs.at(1, sorted_indices[i]) >= n_cols), "SpMat::SpMat(): invalid column index");
Chris@49:         
Chris@49:         access::rw(values[i])      = vals[sorted_indices[i]];
Chris@49:         access::rw(row_indices[i]) = locs.at(0, sorted_indices[i]);
Chris@49:         
Chris@49:         access::rw(col_ptrs[locs.at(1, sorted_indices[i]) + 1])++;
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   if( (sort_locations == false) || (actually_sorted == true) )
Chris@49:     {
Chris@49:     // Now set the values and row indices correctly.
Chris@49:     // Increment the column pointers in each column (so they are column "counts").
Chris@49:     for (uword i = 0; i < vals.n_elem; ++i)
Chris@49:       {
Chris@49:       arma_debug_check((locs.at(0, i) >= n_rows), "SpMat::SpMat(): invalid row index");
Chris@49:       arma_debug_check((locs.at(1, i) >= n_cols), "SpMat::SpMat(): invalid column index");
Chris@49:       
Chris@49:       // Check ordering in debug mode.
Chris@49:       if(i > 0)
Chris@49:         {
Chris@49:         arma_debug_check
Chris@49:           (
Chris@49:           ( (locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) < locs.at(0, i - 1)) ),
Chris@49:           "SpMat::SpMat(): out of order points; either pass sort_locations = true, or sort points in column-major ordering"
Chris@49:           );
Chris@49:         arma_debug_check((locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) == locs.at(0, i - 1)), "SpMat::SpMat(): two identical point locations in list");
Chris@49:         }
Chris@49:       
Chris@49:       access::rw(values[i])      = vals[i];
Chris@49:       access::rw(row_indices[i]) = locs.at(0, i);
Chris@49:       
Chris@49:       access::rw(col_ptrs[locs.at(1, i) + 1])++;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   // Now fix the column pointers.
Chris@49:   for (uword i = 0; i <= n_cols; ++i)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[i + 1]) += col_ptrs[i];
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Insert a large number of values at once.
Chris@49: //! locations.row[0] should be row indices, locations.row[1] should be column indices,
Chris@49: //! and values should be the corresponding values.
Chris@49: //! If sort_locations is false, then it is assumed that the locations and values
Chris@49: //! are already sorted in column-major ordering.
Chris@49: //! In this constructor the size is explicitly given.
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const Base<uword,T1>& locations_expr, const Base<eT,T2>& vals_expr, const uword in_n_rows, const uword in_n_cols, const bool sort_locations)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL)
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49:   
Chris@49:   init(in_n_rows, in_n_cols);
Chris@49:   
Chris@49:   const unwrap<T1>         locs_tmp( locations_expr.get_ref() );
Chris@49:   const Mat<uword>& locs = locs_tmp.M;
Chris@49:   
Chris@49:   const unwrap<T2> vals_tmp( vals_expr.get_ref() );
Chris@49:   const Mat<eT>& vals = vals_tmp.M;
Chris@49:   
Chris@49:   arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'values' object is not a vector" );
Chris@49:   
Chris@49:   arma_debug_check((locs.n_rows != 2), "SpMat::SpMat(): locations matrix must have two rows");
Chris@49:   
Chris@49:   arma_debug_check((locs.n_cols != vals.n_elem), "SpMat::SpMat(): number of locations is different than number of values");
Chris@49:   
Chris@49:   // Resize to correct number of elements.
Chris@49:   mem_resize(vals.n_elem);
Chris@49:   
Chris@49:   // Reset column pointers to zero.
Chris@49:   arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1);
Chris@49:   
Chris@49:   bool actually_sorted = true;
Chris@49:   if(sort_locations == true)
Chris@49:     {
Chris@49:     // sort_index() uses std::sort() which may use quicksort... so we better
Chris@49:     // make sure it's not already sorted before taking an O(N^2) sort penalty.
Chris@49:     for (uword i = 1; i < locs.n_cols; ++i)
Chris@49:       {
Chris@49:       if ((locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) <= locs.at(0, i - 1)))
Chris@49:         {
Chris@49:         actually_sorted = false;
Chris@49:         break;
Chris@49:         }
Chris@49:       }
Chris@49:     
Chris@49:     if(actually_sorted == false)
Chris@49:       {
Chris@49:       // This may not be the fastest possible implementation but it maximizes code reuse.
Chris@49:       Col<uword> abslocs(locs.n_cols);
Chris@49:       
Chris@49:       for (uword i = 0; i < locs.n_cols; ++i)
Chris@49:         {
Chris@49:         abslocs[i] = locs.at(1, i) * n_rows + locs.at(0, i);
Chris@49:         }
Chris@49: 
Chris@49:       // Now we will sort with sort_index().
Chris@49:       uvec sorted_indices = sort_index(abslocs); // Ascending sort.
Chris@49: 
Chris@49:       // Now we add the elements in this sorted order.
Chris@49:       for (uword i = 0; i < sorted_indices.n_elem; ++i)
Chris@49:         {
Chris@49:         arma_debug_check((locs.at(0, sorted_indices[i]) >= n_rows), "SpMat::SpMat(): invalid row index");
Chris@49:         arma_debug_check((locs.at(1, sorted_indices[i]) >= n_cols), "SpMat::SpMat(): invalid column index");
Chris@49: 
Chris@49:         access::rw(values[i])      = vals[sorted_indices[i]];
Chris@49:         access::rw(row_indices[i]) = locs.at(0, sorted_indices[i]);
Chris@49: 
Chris@49:         access::rw(col_ptrs[locs.at(1, sorted_indices[i]) + 1])++;
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   if( (sort_locations == false) || (actually_sorted == true) )
Chris@49:     {
Chris@49:     // Now set the values and row indices correctly.
Chris@49:     // Increment the column pointers in each column (so they are column "counts").
Chris@49:     for (uword i = 0; i < vals.n_elem; ++i)
Chris@49:       {
Chris@49:       arma_debug_check((locs.at(0, i) >= n_rows), "SpMat::SpMat(): invalid row index");
Chris@49:       arma_debug_check((locs.at(1, i) >= n_cols), "SpMat::SpMat(): invalid column index");
Chris@49: 
Chris@49:       // Check ordering in debug mode.
Chris@49:       if(i > 0)
Chris@49:         {
Chris@49:         arma_debug_check
Chris@49:           (
Chris@49:           ( (locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) < locs.at(0, i - 1)) ),
Chris@49:           "SpMat::SpMat(): out of order points; either pass sort_locations = true or sort points in column-major ordering"
Chris@49:           );
Chris@49:         arma_debug_check((locs.at(1, i) == locs.at(1, i - 1) && locs.at(0, i) == locs.at(0, i - 1)), "SpMat::SpMat(): two identical point locations in list");
Chris@49:         }
Chris@49: 
Chris@49:       access::rw(values[i])      = vals[i];
Chris@49:       access::rw(row_indices[i]) = locs.at(0, i);
Chris@49: 
Chris@49:       access::rw(col_ptrs[locs.at(1, i) + 1])++;
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Now fix the column pointers.
Chris@49:   for (uword i = 0; i <= n_cols; ++i)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[i + 1]) += col_ptrs[i];
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Simple operators with plain values.  These operate on every value in the
Chris@49:  * matrix, so a sparse matrix += 1 will turn all those zeroes into ones.  Be
Chris@49:  * careful and make sure that's what you really want!
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const eT val)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   // Resize to 1x1 then set that to the right value.
Chris@49:   init(1, 1); // Sets col_ptrs to 0.
Chris@49:   mem_resize(1); // One element.
Chris@49: 
Chris@49:   // Manually set element.
Chris@49:   access::rw(values[0]) = val;
Chris@49:   access::rw(row_indices[0]) = 0;
Chris@49:   access::rw(col_ptrs[1]) = 1;
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const eT val)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   if(val == eT(0))
Chris@49:     {
Chris@49:     // Everything will be zero.
Chris@49:     init(n_rows, n_cols);
Chris@49:     return *this;
Chris@49:     }
Chris@49:   
Chris@49:   arrayops::inplace_mul( access::rwp(values), val, n_nonzero );
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const eT val)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_check( (val == eT(0)), "element-wise division: division by zero" );
Chris@49:   
Chris@49:   arrayops::inplace_div( access::rwp(values), val, n_nonzero );
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   init(x);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "addition");
Chris@49:   
Chris@49:   // Iterate over nonzero values of other matrix.
Chris@49:   for (const_iterator it = x.begin(); it != x.end(); it++)
Chris@49:     {
Chris@49:     get_value(it.row(), it.col()) += *it;
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "subtraction");
Chris@49:   
Chris@49:   // Iterate over nonzero values of other matrix.
Chris@49:   for (const_iterator it = x.begin(); it != x.end(); it++)
Chris@49:     {
Chris@49:     get_value(it.row(), it.col()) -= *it;
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const SpMat<eT>& y)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_mul_size(n_rows, n_cols, y.n_rows, y.n_cols, "matrix multiplication");
Chris@49:   
Chris@49:   SpMat<eT> z;
Chris@49:   z = (*this) * y;
Chris@49:   steal_mem(z);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: // This is in-place element-wise matrix multiplication.
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise multiplication");
Chris@49:   
Chris@49:   // We can do this with two iterators rather simply.
Chris@49:         iterator it    = begin();
Chris@49:   const_iterator x_it  = x.begin();
Chris@49:   
Chris@49:   while (it != end() && x_it != x.end())
Chris@49:     {
Chris@49:     // One of these will be further advanced than the other (or they will be at the same place).
Chris@49:     if ((it.row() == x_it.row()) && (it.col() == x_it.col()))
Chris@49:       {
Chris@49:       // There is an element at this place in both matrices.  Multiply.
Chris@49:       (*it) *= (*x_it);
Chris@49: 
Chris@49:       // Now move on to the next position.
Chris@49:       it++;
Chris@49:       x_it++;
Chris@49:       }
Chris@49: 
Chris@49:     else if ((it.col() < x_it.col()) || ((it.col() == x_it.col()) && (it.row() < x_it.row())))
Chris@49:       {
Chris@49:       // This case is when our matrix has an element which the other matrix does not.
Chris@49:       // So we must delete this element.
Chris@49:       (*it) = 0;
Chris@49: 
Chris@49:       // Because we have deleted the element, we now have to manually set the position...
Chris@49:       it.internal_pos--;
Chris@49: 
Chris@49:       // Now we can increment our iterator.
Chris@49:       it++;
Chris@49:       }
Chris@49: 
Chris@49:     else /* if our iterator is ahead of the other matrix */
Chris@49:       {
Chris@49:       // In this case we don't need to set anything to 0; our element is already 0.
Chris@49:       // We can just increment the iterator of the other matrix.
Chris@49:       x_it++;
Chris@49:       }
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   // If we are not at the end of our matrix, then we must terminate the remaining elements.
Chris@49:   while (it != end())
Chris@49:     {
Chris@49:     (*it) = 0;
Chris@49: 
Chris@49:     // Hack to manually set the position right...
Chris@49:     it.internal_pos--;
Chris@49:     it++; // ...and then an increment.
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: // Construct a complex matrix out of two non-complex matrices
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat
Chris@49:   (
Chris@49:   const SpBase<typename SpMat<eT>::pod_type, T1>& A,
Chris@49:   const SpBase<typename SpMat<eT>::pod_type, T2>& B
Chris@49:   )
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra element is set when mem_resize is called
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type T;
Chris@49:   
Chris@49:   // Make sure eT is complex and T is not (compile-time check).
Chris@49:   arma_type_check(( is_complex<eT>::value == false ));
Chris@49:   arma_type_check(( is_complex< T>::value == true  ));
Chris@49:   
Chris@49:   // Compile-time abort if types are not compatible.
Chris@49:   arma_type_check(( is_same_type< std::complex<T>, eT >::value == false ));
Chris@49:   
Chris@49:   const unwrap_spmat<T1> tmp1(A.get_ref());
Chris@49:   const unwrap_spmat<T2> tmp2(B.get_ref());
Chris@49:   
Chris@49:   const SpMat<T>& X = tmp1.M;
Chris@49:   const SpMat<T>& Y = tmp2.M;
Chris@49:   
Chris@49:   arma_debug_assert_same_size(X.n_rows, X.n_cols, Y.n_rows, Y.n_cols, "SpMat()");
Chris@49:   
Chris@49:   const uword l_n_rows = X.n_rows;
Chris@49:   const uword l_n_cols = X.n_cols;
Chris@49:   
Chris@49:   // Set size of matrix correctly.
Chris@49:   init(l_n_rows, l_n_cols);
Chris@49:   mem_resize(n_unique(X, Y, op_n_unique_count()));
Chris@49:   
Chris@49:   // Now on a second iteration, fill it.
Chris@49:   typename SpMat<T>::const_iterator x_it  = X.begin();
Chris@49:   typename SpMat<T>::const_iterator x_end = X.end();
Chris@49:   
Chris@49:   typename SpMat<T>::const_iterator y_it  = Y.begin();
Chris@49:   typename SpMat<T>::const_iterator y_end = Y.end();
Chris@49:   
Chris@49:   uword cur_pos = 0;
Chris@49:   
Chris@49:   while ((x_it != x_end) || (y_it != y_end))
Chris@49:     {
Chris@49:     if(x_it == y_it) // if we are at the same place
Chris@49:       {
Chris@49:       access::rw(values[cur_pos]) = std::complex<T>((T) *x_it, (T) *y_it);
Chris@49:       access::rw(row_indices[cur_pos]) = x_it.row();
Chris@49:       ++access::rw(col_ptrs[x_it.col() + 1]);
Chris@49:       
Chris@49:       ++x_it;
Chris@49:       ++y_it;
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       if((x_it.col() < y_it.col()) || ((x_it.col() == y_it.col()) && (x_it.row() < y_it.row()))) // if y is closer to the end
Chris@49:         {
Chris@49:         access::rw(values[cur_pos]) = std::complex<T>((T) *x_it, T(0));
Chris@49:         access::rw(row_indices[cur_pos]) = x_it.row();
Chris@49:         ++access::rw(col_ptrs[x_it.col() + 1]);
Chris@49:         
Chris@49:         ++x_it;
Chris@49:         }
Chris@49:       else // x is closer to the end
Chris@49:         {
Chris@49:         access::rw(values[cur_pos]) = std::complex<T>(T(0), (T) *y_it);
Chris@49:         access::rw(row_indices[cur_pos]) = y_it.row();
Chris@49:         ++access::rw(col_ptrs[y_it.col() + 1]);
Chris@49:         
Chris@49:         ++y_it;
Chris@49:         }
Chris@49:       }
Chris@49:     
Chris@49:     ++cur_pos;
Chris@49:     }
Chris@49:   
Chris@49:   // Now fix the column pointers; they are supposed to be a sum.
Chris@49:   for (uword c = 1; c <= n_cols; ++c)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[c]) += col_ptrs[c - 1];
Chris@49:     }
Chris@49:   
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise division");
Chris@49:   
Chris@49:   // If you use this method, you are probably stupid or misguided, but for compatibility with Mat, we have implemented it anyway.
Chris@49:   // We have to loop over every element, which is not good.  In fact, it makes me physically sad to write this.
Chris@49:   for(uword c = 0; c < n_cols; ++c)
Chris@49:     {
Chris@49:     for(uword r = 0; r < n_rows; ++r)
Chris@49:       {
Chris@49:       at(r, c) /= x.at(r, c);
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const Base<eT, T1>& x)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra element is set when mem_resize is called in operator=()
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   (*this).operator=(x);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const Base<eT, T1>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const Proxy<T1> p(x.get_ref());
Chris@49:   
Chris@49:   const uword x_n_rows = p.get_n_rows();
Chris@49:   const uword x_n_cols = p.get_n_cols();
Chris@49:   const uword x_n_elem = p.get_n_elem();
Chris@49: 
Chris@49:   init(x_n_rows, x_n_cols);
Chris@49: 
Chris@49:   // Count number of nonzero elements in base object.
Chris@49:   uword n = 0;
Chris@49:   if(Proxy<T1>::prefer_at_accessor == true)
Chris@49:     {
Chris@49:     for(uword j = 0; j < x_n_cols; ++j)
Chris@49:     for(uword i = 0; i < x_n_rows; ++i)
Chris@49:       {
Chris@49:       if(p.at(i, j) != eT(0)) { ++n; }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     for(uword i = 0; i < x_n_elem; ++i)
Chris@49:       {
Chris@49:       if(p[i] != eT(0)) { ++n; }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   mem_resize(n);
Chris@49: 
Chris@49:   // Now the memory is resized correctly; add nonzero elements.
Chris@49:   n = 0;
Chris@49:   for(uword j = 0; j < x_n_cols; ++j)
Chris@49:   for(uword i = 0; i < x_n_rows; ++i)
Chris@49:     {
Chris@49:     const eT val = p.at(i, j);
Chris@49:     
Chris@49:     if(val != eT(0))
Chris@49:       {
Chris@49:       access::rw(values[n])      = val;
Chris@49:       access::rw(row_indices[n]) = i;
Chris@49:       access::rw(col_ptrs[j + 1])++;
Chris@49:       ++n;
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Sum column counts to be column pointers.
Chris@49:   for(uword c = 1; c <= n_cols; ++c)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[c]) += col_ptrs[c - 1];
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const Base<eT, T1>& y)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const Proxy<T1> p(y.get_ref());
Chris@49: 
Chris@49:   arma_debug_assert_mul_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "matrix multiplication");
Chris@49: 
Chris@49:   // We assume the matrix structure is such that we will end up with a sparse
Chris@49:   // matrix.  Assuming that every entry in the dense matrix is nonzero (which is
Chris@49:   // a fairly valid assumption), each row with any nonzero elements in it (in this
Chris@49:   // matrix) implies an entire nonzero column.  Therefore, we iterate over all
Chris@49:   // the row_indices and count the number of rows with any elements in them
Chris@49:   // (using the quasi-linked-list idea from SYMBMM -- see operator_times.hpp).
Chris@49:   podarray<uword> index(n_rows);
Chris@49:   index.fill(n_rows); // Fill with invalid links.
Chris@49: 
Chris@49:   uword last_index = n_rows + 1;
Chris@49:   for(uword i = 0; i < n_nonzero; ++i)
Chris@49:     {
Chris@49:     if(index[row_indices[i]] == n_rows)
Chris@49:       {
Chris@49:       index[row_indices[i]] = last_index;
Chris@49:       last_index = row_indices[i];
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Now count the number of rows which have nonzero elements.
Chris@49:   uword nonzero_rows = 0;
Chris@49:   while(last_index != n_rows + 1)
Chris@49:     {
Chris@49:     ++nonzero_rows;
Chris@49:     last_index = index[last_index];
Chris@49:     }
Chris@49: 
Chris@49:   SpMat<eT> z(n_rows, p.get_n_cols());
Chris@49: 
Chris@49:   z.mem_resize(nonzero_rows * p.get_n_cols()); // upper bound on size
Chris@49: 
Chris@49:   // Now we have to fill all the elements using a modification of the NUMBMM algorithm.
Chris@49:   uword cur_pos = 0;
Chris@49: 
Chris@49:   podarray<eT> partial_sums(n_rows);
Chris@49:   partial_sums.zeros();
Chris@49: 
Chris@49:   for(uword lcol = 0; lcol < n_cols; ++lcol)
Chris@49:     {
Chris@49:     const_iterator it = begin();
Chris@49: 
Chris@49:     while(it != end())
Chris@49:       {
Chris@49:       const eT value = (*it);
Chris@49: 
Chris@49:       partial_sums[it.row()] += (value * p.at(it.col(), lcol));
Chris@49: 
Chris@49:       ++it;
Chris@49:       }
Chris@49: 
Chris@49:     // Now add all partial sums to the matrix.
Chris@49:     for(uword i = 0; i < n_rows; ++i)
Chris@49:       {
Chris@49:       if(partial_sums[i] != eT(0))
Chris@49:         {
Chris@49:         access::rw(z.values[cur_pos]) = partial_sums[i];
Chris@49:         access::rw(z.row_indices[cur_pos]) = i;
Chris@49:         ++access::rw(z.col_ptrs[lcol + 1]);
Chris@49:         //printf("colptr %d now %d\n", lcol + 1, z.col_ptrs[lcol + 1]);
Chris@49:         ++cur_pos;
Chris@49:         partial_sums[i] = 0; // Would it be faster to do this in batch later?
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Now fix the column pointers.
Chris@49:   for(uword c = 1; c <= z.n_cols; ++c)
Chris@49:     {
Chris@49:     access::rw(z.col_ptrs[c]) += z.col_ptrs[c - 1];
Chris@49:     }
Chris@49: 
Chris@49:   // Resize to final correct size.
Chris@49:   z.mem_resize(z.col_ptrs[z.n_cols]);
Chris@49:   
Chris@49:   // Now take the memory of the temporary matrix.
Chris@49:   steal_mem(z);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Don't use this function.  It's not mathematically well-defined and wastes
Chris@49:  * cycles to trash all your data.  This is dumb.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const Base<eT, T1>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   SpMat<eT> tmp = (*this) / x.get_ref();
Chris@49:   
Chris@49:   steal_mem(tmp);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const Base<eT, T1>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const Proxy<T1> p(x.get_ref());
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "element-wise multiplication");
Chris@49:   
Chris@49:   // Count the number of elements we will need.
Chris@49:   SpMat<eT> tmp(n_rows, n_cols);
Chris@49:   const_iterator it = begin();
Chris@49:   uword new_n_nonzero = 0;
Chris@49: 
Chris@49:   while(it != end())
Chris@49:     {
Chris@49:     // prefer_at_accessor == false can't save us any work here
Chris@49:     if(((*it) * p.at(it.row(), it.col())) != eT(0))
Chris@49:       {
Chris@49:       ++new_n_nonzero;
Chris@49:       }
Chris@49:     ++it;
Chris@49:     }
Chris@49: 
Chris@49:   // Resize.
Chris@49:   tmp.mem_resize(new_n_nonzero);
Chris@49: 
Chris@49:   const_iterator c_it = begin();
Chris@49:   uword cur_pos = 0;
Chris@49:   while(c_it != end())
Chris@49:     {
Chris@49:     // prefer_at_accessor == false can't save us any work here
Chris@49:     const eT val = (*c_it) * p.at(c_it.row(), c_it.col());
Chris@49:     if(val != eT(0))
Chris@49:       {
Chris@49:       access::rw(tmp.values[cur_pos]) = val;
Chris@49:       access::rw(tmp.row_indices[cur_pos]) = c_it.row();
Chris@49:       ++access::rw(tmp.col_ptrs[c_it.col() + 1]);
Chris@49:       ++cur_pos;
Chris@49:       }
Chris@49: 
Chris@49:     ++c_it;
Chris@49:     }
Chris@49: 
Chris@49:   // Fix column pointers.
Chris@49:   for(uword c = 1; c <= n_cols; ++c)
Chris@49:     {
Chris@49:     access::rw(tmp.col_ptrs[c]) += tmp.col_ptrs[c - 1];
Chris@49:     }
Chris@49: 
Chris@49:   steal_mem(tmp);
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Functions on subviews.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const SpSubview<eT>& X)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra element added when mem_resize is called
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   (*this).operator=(X);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const SpSubview<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword in_n_cols = X.n_cols;
Chris@49:   const uword in_n_rows = X.n_rows;
Chris@49:   
Chris@49:   const bool alias = (this == &(X.m));
Chris@49: 
Chris@49:   if(alias == false)
Chris@49:     {
Chris@49:     init(in_n_rows, in_n_cols);
Chris@49: 
Chris@49:     const uword x_n_nonzero = X.n_nonzero;
Chris@49: 
Chris@49:     mem_resize(x_n_nonzero);
Chris@49: 
Chris@49:     typename SpSubview<eT>::const_iterator it = X.begin();
Chris@49: 
Chris@49:     while(it != X.end())
Chris@49:       {
Chris@49:       access::rw(row_indices[it.pos()]) = it.row();
Chris@49:       access::rw(values[it.pos()]) = (*it);
Chris@49:       ++access::rw(col_ptrs[it.col() + 1]);
Chris@49:       ++it;
Chris@49:       }
Chris@49: 
Chris@49:     // Now sum column pointers.
Chris@49:     for(uword c = 1; c <= n_cols; ++c)
Chris@49:       {
Chris@49:       access::rw(col_ptrs[c]) += col_ptrs[c - 1];
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // Create it in a temporary.
Chris@49:     SpMat<eT> tmp(X);
Chris@49: 
Chris@49:     steal_mem(tmp);
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const SpSubview<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition");
Chris@49:   
Chris@49:   typename SpSubview<eT>::const_iterator it = X.begin();
Chris@49: 
Chris@49:   while(it != X.end())
Chris@49:     {
Chris@49:     at(it.row(), it.col()) += (*it);
Chris@49:     ++it;
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const SpSubview<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtraction");
Chris@49:   
Chris@49:   typename SpSubview<eT>::const_iterator it = X.begin();
Chris@49:   
Chris@49:   while(it != X.end())
Chris@49:     {
Chris@49:     at(it.row(), it.col()) -= (*it);
Chris@49:     ++it;
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const SpSubview<eT>& y)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_mul_size(n_rows, n_cols, y.n_rows, y.n_cols, "matrix multiplication");
Chris@49:   
Chris@49:   // Cannot be done in-place (easily).
Chris@49:   SpMat<eT> z = (*this) * y;
Chris@49:   steal_mem(z);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const SpSubview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise multiplication");
Chris@49: 
Chris@49:   iterator it = begin();
Chris@49:   typename SpSubview<eT>::const_iterator xit = x.begin();
Chris@49: 
Chris@49:   while((it != end()) || (xit != x.end()))
Chris@49:     {
Chris@49:     if((xit.row() == it.row()) && (xit.col() == it.col()))
Chris@49:       {
Chris@49:       (*it) *= (*xit);
Chris@49:       ++it;
Chris@49:       ++xit;
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       if((xit.col() > it.col()) || ((xit.col() == it.col()) && (xit.row() > it.row())))
Chris@49:         {
Chris@49:         // xit is "ahead"
Chris@49:         (*it) = eT(0); // erase element; x has a zero here
Chris@49:         it.internal_pos--; // update iterator so it still works
Chris@49:         ++it;
Chris@49:         }
Chris@49:       else
Chris@49:         {
Chris@49:         // it is "ahead"
Chris@49:         ++xit;
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const SpSubview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise division");
Chris@49:   
Chris@49:   // There is no pretty way to do this.
Chris@49:   for(uword elem = 0; elem < n_elem; elem++)
Chris@49:     {
Chris@49:     at(elem) /= x(elem);
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Operators on regular subviews.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const subview<eT>& x)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra value set in operator=()
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   (*this).operator=(x);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const subview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword x_n_rows = x.n_rows;
Chris@49:   const uword x_n_cols = x.n_cols;
Chris@49:   
Chris@49:   // Set the size correctly.
Chris@49:   init(x_n_rows, x_n_cols);
Chris@49: 
Chris@49:   // Count number of nonzero elements.
Chris@49:   uword n = 0;
Chris@49:   for(uword c = 0; c < x_n_cols; ++c)
Chris@49:     {
Chris@49:     for(uword r = 0; r < x_n_rows; ++r)
Chris@49:       {
Chris@49:       if(x.at(r, c) != eT(0))
Chris@49:         {
Chris@49:         ++n;
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Resize memory appropriately.
Chris@49:   mem_resize(n);
Chris@49: 
Chris@49:   n = 0;
Chris@49:   for(uword c = 0; c < x_n_cols; ++c)
Chris@49:     {
Chris@49:     for(uword r = 0; r < x_n_rows; ++r)
Chris@49:       {
Chris@49:       const eT val = x.at(r, c);
Chris@49:       
Chris@49:       if(val != eT(0))
Chris@49:         {
Chris@49:         access::rw(values[n]) = val;
Chris@49:         access::rw(row_indices[n]) = r;
Chris@49:         ++access::rw(col_ptrs[c + 1]);
Chris@49:         ++n;
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Fix column counts into column pointers.
Chris@49:   for(uword c = 1; c <= n_cols; ++c)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[c]) += col_ptrs[c - 1];
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const subview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "addition");
Chris@49:   
Chris@49:   // Loop over every element.  This could probably be written in a more
Chris@49:   // efficient way, by calculating the number of nonzero elements the output
Chris@49:   // matrix will have, allocating the memory correctly, and then filling the
Chris@49:   // matrix correctly.  However... for now, this works okay.
Chris@49:   for(uword lcol = 0; lcol < n_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < n_rows; ++lrow)
Chris@49:     {
Chris@49:     at(lrow, lcol) += x.at(lrow, lcol);
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const subview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "subtraction");
Chris@49:   
Chris@49:   // Loop over every element.
Chris@49:   for(uword lcol = 0; lcol < n_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < n_rows; ++lrow)
Chris@49:     {
Chris@49:     at(lrow, lcol) -= x.at(lrow, lcol);
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const subview<eT>& y)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_assert_mul_size(n_rows, n_cols, y.n_rows, y.n_cols, "matrix multiplication");
Chris@49: 
Chris@49:   SpMat<eT> z(n_rows, y.n_cols);
Chris@49: 
Chris@49:   // Performed in the same fashion as operator*=(SpMat).
Chris@49:   for (const_row_iterator x_row_it = begin_row(); x_row_it.pos() < n_nonzero; ++x_row_it)
Chris@49:     {
Chris@49:     for (uword lcol = 0; lcol < y.n_cols; ++lcol)
Chris@49:       {
Chris@49:       // At this moment in the loop, we are calculating anything that is contributed to by *x_row_it and *y_col_it.
Chris@49:       // Given that our position is x_ab and y_bc, there will only be a contribution if x.col == y.row, and that
Chris@49:       // contribution will be in location z_ac.
Chris@49:       z.at(x_row_it.row, lcol) += (*x_row_it) * y.at(x_row_it.col, lcol);
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   steal_mem(z);
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const subview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise multiplication");
Chris@49:   
Chris@49:   // Loop over every element.
Chris@49:   for(uword lcol = 0; lcol < n_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < n_rows; ++lrow)
Chris@49:     {
Chris@49:     at(lrow, lcol) *= x.at(lrow, lcol);
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const subview<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise division");
Chris@49:   
Chris@49:   // Loop over every element.
Chris@49:   for(uword lcol = 0; lcol < n_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < n_rows; ++lrow)
Chris@49:     {
Chris@49:     at(lrow, lcol) /= x.at(lrow, lcol);
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const SpOp<T1, spop_type>& X)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // set in application of sparse operation
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   spop_type::apply(*this, X);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   spop_type::apply(*this, X);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator+=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator-=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator*=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator%=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const SpOp<T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator/=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra element set in application of sparse glue
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   spglue_type::apply(*this, X);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: SpMat<eT>::SpMat(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   : n_rows(0)
Chris@49:   , n_cols(0)
Chris@49:   , n_elem(0)
Chris@49:   , n_nonzero(0)
Chris@49:   , vec_state(0)
Chris@49:   , values(NULL) // extra element set in application of sparse glue
Chris@49:   , row_indices(NULL)
Chris@49:   , col_ptrs(NULL)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint_this(this);
Chris@49: 
Chris@49:   spop_type::apply(*this, X);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   spop_type::apply(*this, X);
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const SpMat<eT> m(X);
Chris@49: 
Chris@49:   return (*this).operator+=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const SpMat<eT> m(X);
Chris@49: 
Chris@49:   return (*this).operator-=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const SpMat<eT> m(X);
Chris@49: 
Chris@49:   return (*this).operator*=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const SpMat<eT> m(X);
Chris@49: 
Chris@49:   return (*this).operator%=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename spop_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const mtSpOp<eT, T1, spop_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const SpMat<eT> m(X);
Chris@49: 
Chris@49:   return (*this).operator/=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   spglue_type::apply(*this, X);
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator+=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator+=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator-=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator-=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator*=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator*=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator%=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator%=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename T2, typename spglue_type>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::operator/=(const SpGlue<T1, T2, spglue_type>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_type_check(( is_same_type< eT, typename T1::elem_type >::value == false ));
Chris@49:   
Chris@49:   const SpMat<eT> m(X);
Chris@49:   
Chris@49:   return (*this).operator/=(m);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::row(const uword row_num)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check(row_num >= n_rows, "SpMat::row(): out of bounds");
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, row_num, 0, 1, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::row(const uword row_num) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check(row_num >= n_rows, "SpMat::row(): out of bounds");
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, row_num, 0, 1, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const uword row_num, const span& col_span)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const bool col_all = col_span.whole;
Chris@49:   
Chris@49:   const uword local_n_cols = n_cols;
Chris@49:   
Chris@49:   const uword in_col1       = col_all ? 0            : col_span.a;
Chris@49:   const uword in_col2       =                          col_span.b;
Chris@49:   const uword submat_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1;
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (row_num >= n_rows)
Chris@49:     ||
Chris@49:     ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) )
Chris@49:     ,
Chris@49:     "SpMat::operator(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, row_num, in_col1, 1, submat_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const uword row_num, const span& col_span) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const bool col_all = col_span.whole;
Chris@49:   
Chris@49:   const uword local_n_cols = n_cols;
Chris@49:   
Chris@49:   const uword in_col1       = col_all ? 0            : col_span.a;
Chris@49:   const uword in_col2       =                          col_span.b;
Chris@49:   const uword submat_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1;
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (row_num >= n_rows)
Chris@49:     ||
Chris@49:     ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) )
Chris@49:     ,
Chris@49:     "SpMat::operator(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, row_num, in_col1, 1, submat_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::col(const uword col_num)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check(col_num >= n_cols, "SpMat::col(): out of bounds");
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, 0, col_num, n_rows, 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::col(const uword col_num) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check(col_num >= n_cols, "SpMat::col(): out of bounds");
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, 0, col_num, n_rows, 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const span& row_span, const uword col_num)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const bool row_all = row_span.whole;
Chris@49:   
Chris@49:   const uword local_n_rows = n_rows;
Chris@49:   
Chris@49:   const uword in_row1       = row_all ? 0            : row_span.a;
Chris@49:   const uword in_row2       =                          row_span.b;
Chris@49:   const uword submat_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1;
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (col_num >= n_cols)
Chris@49:     ||
Chris@49:     ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) )
Chris@49:     ,
Chris@49:     "SpMat::operator(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, in_row1, col_num, submat_n_rows, 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const span& row_span, const uword col_num) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const bool row_all = row_span.whole;
Chris@49:   
Chris@49:   const uword local_n_rows = n_rows;
Chris@49:   
Chris@49:   const uword in_row1       = row_all ? 0            : row_span.a;
Chris@49:   const uword in_row2       =                          row_span.b;
Chris@49:   const uword submat_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1;
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (col_num >= n_cols)
Chris@49:     ||
Chris@49:     ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) )
Chris@49:     ,
Chris@49:     "SpMat::operator(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, in_row1, col_num, submat_n_rows, 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Swap in_row1 with in_row2.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::swap_rows(const uword in_row1, const uword in_row2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 >= n_rows) || (in_row2 >= n_rows),
Chris@49:     "SpMat::swap_rows(): out of bounds"
Chris@49:     );
Chris@49: 
Chris@49:   // Sanity check.
Chris@49:   if (in_row1 == in_row2)
Chris@49:     {
Chris@49:     return;
Chris@49:     }
Chris@49: 
Chris@49:   // The easier way to do this, instead of collecting all the elements in one row and then swapping with the other, will be
Chris@49:   // to iterate over each column of the matrix (since we store in column-major format) and then swap the two elements in the two rows at that time.
Chris@49:   // We will try to avoid using the at() call since it is expensive, instead preferring to use an iterator to track our position.
Chris@49:   uword col1 = (in_row1 < in_row2) ? in_row1 : in_row2;
Chris@49:   uword col2 = (in_row1 < in_row2) ? in_row2 : in_row1;
Chris@49: 
Chris@49:   for (uword lcol = 0; lcol < n_cols; lcol++)
Chris@49:     {
Chris@49:     // If there is nothing in this column we can ignore it.
Chris@49:     if (col_ptrs[lcol] == col_ptrs[lcol + 1])
Chris@49:       {
Chris@49:       continue;
Chris@49:       }
Chris@49: 
Chris@49:     // These will represent the positions of the items themselves.
Chris@49:     uword loc1 = n_nonzero + 1;
Chris@49:     uword loc2 = n_nonzero + 1;
Chris@49: 
Chris@49:     for (uword search_pos = col_ptrs[lcol]; search_pos < col_ptrs[lcol + 1]; search_pos++)
Chris@49:       {
Chris@49:       if (row_indices[search_pos] == col1)
Chris@49:         {
Chris@49:         loc1 = search_pos;
Chris@49:         }
Chris@49: 
Chris@49:       if (row_indices[search_pos] == col2)
Chris@49:         {
Chris@49:         loc2 = search_pos;
Chris@49:         break; // No need to look any further.
Chris@49:         }
Chris@49:       }
Chris@49: 
Chris@49:     // There are four cases: we found both elements; we found one element (loc1); we found one element (loc2); we found zero elements.
Chris@49:     // If we found zero elements no work needs to be done and we can continue to the next column.
Chris@49:     if ((loc1 != (n_nonzero + 1)) && (loc2 != (n_nonzero + 1)))
Chris@49:       {
Chris@49:       // This is an easy case: just swap the values.  No index modifying necessary.
Chris@49:       eT tmp = values[loc1];
Chris@49:       access::rw(values[loc1]) = values[loc2];
Chris@49:       access::rw(values[loc2]) = tmp;
Chris@49:       }
Chris@49:     else if (loc1 != (n_nonzero + 1)) // We only found loc1 and not loc2.
Chris@49:       {
Chris@49:       // We need to find the correct place to move our value to.  It will be forward (not backwards) because in_row2 > in_row1.
Chris@49:       // Each iteration of the loop swaps the current value (loc1) with (loc1 + 1); in this manner we move our value down to where it should be.
Chris@49:       while (((loc1 + 1) < col_ptrs[lcol + 1]) && (row_indices[loc1 + 1] < in_row2))
Chris@49:         {
Chris@49:         // Swap both the values and the indices.  The column should not change.
Chris@49:         eT tmp = values[loc1];
Chris@49:         access::rw(values[loc1]) = values[loc1 + 1];
Chris@49:         access::rw(values[loc1 + 1]) = tmp;
Chris@49: 
Chris@49:         uword tmp_index = row_indices[loc1];
Chris@49:         access::rw(row_indices[loc1]) = row_indices[loc1 + 1];
Chris@49:         access::rw(row_indices[loc1 + 1]) = tmp_index;
Chris@49: 
Chris@49:         loc1++; // And increment the counter.
Chris@49:         }
Chris@49: 
Chris@49:       // Now set the row index correctly.
Chris@49:       access::rw(row_indices[loc1]) = in_row2;
Chris@49: 
Chris@49:       }
Chris@49:     else if (loc2 != (n_nonzero + 1))
Chris@49:       {
Chris@49:       // We need to find the correct place to move our value to.  It will be backwards (not forwards) because in_row1 < in_row2.
Chris@49:       // Each iteration of the loop swaps the current value (loc2) with (loc2 - 1); in this manner we move our value up to where it should be.
Chris@49:       while (((loc2 - 1) >= col_ptrs[lcol]) && (row_indices[loc2 - 1] > in_row1))
Chris@49:         {
Chris@49:         // Swap both the values and the indices.  The column should not change.
Chris@49:         eT tmp = values[loc2];
Chris@49:         access::rw(values[loc2]) = values[loc2 - 1];
Chris@49:         access::rw(values[loc2 - 1]) = tmp;
Chris@49: 
Chris@49:         uword tmp_index = row_indices[loc2];
Chris@49:         access::rw(row_indices[loc2]) = row_indices[loc2 - 1];
Chris@49:         access::rw(row_indices[loc2 - 1]) = tmp_index;
Chris@49: 
Chris@49:         loc2--; // And decrement the counter.
Chris@49:         }
Chris@49: 
Chris@49:       // Now set the row index correctly.
Chris@49:       access::rw(row_indices[loc2]) = in_row1;
Chris@49: 
Chris@49:       }
Chris@49:     /* else: no need to swap anything; both values are zero */
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Swap in_col1 with in_col2.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::swap_cols(const uword in_col1, const uword in_col2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   // slow but works
Chris@49:   for(uword lrow = 0; lrow < n_rows; ++lrow)
Chris@49:     {
Chris@49:     eT tmp = at(lrow, in_col1);
Chris@49:     at(lrow, in_col1) = at(lrow, in_col2);
Chris@49:     at(lrow, in_col2) = tmp;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Remove the row row_num.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::shed_row(const uword row_num)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_debug_check (row_num >= n_rows, "SpMat::shed_row(): out of bounds");
Chris@49: 
Chris@49:   shed_rows (row_num, row_num);
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Remove the column col_num.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::shed_col(const uword col_num)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_debug_check (col_num >= n_cols, "SpMat::shed_col(): out of bounds");
Chris@49: 
Chris@49:   shed_cols(col_num, col_num);
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Remove all rows between (and including) in_row1 and in_row2.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::shed_rows(const uword in_row1, const uword in_row2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 > in_row2) || (in_row2 >= n_rows),
Chris@49:     "SpMat::shed_rows(): indices out of bounds or incorectly used"
Chris@49:     );
Chris@49: 
Chris@49:   uword i, j;
Chris@49:   // Store the length of values
Chris@49:   uword vlength = n_nonzero;
Chris@49:   // Store the length of col_ptrs
Chris@49:   uword clength = n_cols + 1;
Chris@49: 
Chris@49:   // This is O(n * n_cols) and inplace, there may be a faster way, though.
Chris@49:   for (i = 0, j = 0; i < vlength; ++i)
Chris@49:     {
Chris@49:     // Store the row of the ith element.
Chris@49:     const uword lrow = row_indices[i];
Chris@49:     // Is the ith element in the range of rows we want to remove?
Chris@49:     if (lrow >= in_row1 && lrow <= in_row2)
Chris@49:       {
Chris@49:       // Increment our "removed elements" counter.
Chris@49:       ++j;
Chris@49: 
Chris@49:       // Adjust the values of col_ptrs each time we remove an element.
Chris@49:       // Basically, the length of one column reduces by one, and everything to
Chris@49:       // its right gets reduced by one to represent all the elements being
Chris@49:       // shifted to the left by one.
Chris@49:       for(uword k = 0; k < clength; ++k)
Chris@49:         {
Chris@49:         if (col_ptrs[k] > (i - j + 1))
Chris@49:           {
Chris@49:           --access::rw(col_ptrs[k]);
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       // We shift the element we checked to the left by how many elements
Chris@49:       // we have removed.
Chris@49:       // j = 0 until we remove the first element.
Chris@49:       if (j != 0)
Chris@49:         {
Chris@49:         access::rw(row_indices[i - j]) = (lrow > in_row2) ? (lrow - (in_row2 - in_row1 + 1)) : lrow;
Chris@49:         access::rw(values[i - j]) = values[i];
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // j is the number of elements removed.
Chris@49: 
Chris@49:   // Shrink the vectors.  This will copy the memory.
Chris@49:   mem_resize(n_nonzero - j);
Chris@49: 
Chris@49:   // Adjust row and element counts.
Chris@49:   access::rw(n_rows)    = n_rows - (in_row2 - in_row1) - 1;
Chris@49:   access::rw(n_elem)    = n_rows * n_cols;
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Remove all columns between (and including) in_col1 and in_col2.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::shed_cols(const uword in_col1, const uword in_col2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_col1 > in_col2) || (in_col2 >= n_cols),
Chris@49:     "SpMat::shed_cols(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   // First we find the locations in values and row_indices for the column entries.
Chris@49:   uword col_beg = col_ptrs[in_col1];
Chris@49:   uword col_end = col_ptrs[in_col2 + 1];
Chris@49: 
Chris@49:   // Then we find the number of entries in the column.
Chris@49:   uword diff = col_end - col_beg;
Chris@49: 
Chris@49:   if (diff > 0)
Chris@49:     {
Chris@49:     eT*    new_values      = memory::acquire_chunked<eT>   (n_nonzero - diff);
Chris@49:     uword* new_row_indices = memory::acquire_chunked<uword>(n_nonzero - diff);
Chris@49: 
Chris@49:     // Copy first part.
Chris@49:     if (col_beg != 0)
Chris@49:       {
Chris@49:       arrayops::copy(new_values, values, col_beg);
Chris@49:       arrayops::copy(new_row_indices, row_indices, col_beg);
Chris@49:       }
Chris@49: 
Chris@49:     // Copy second part.
Chris@49:     if (col_end != n_nonzero)
Chris@49:       {
Chris@49:       arrayops::copy(new_values + col_beg, values + col_end, n_nonzero - col_end);
Chris@49:       arrayops::copy(new_row_indices + col_beg, row_indices + col_end, n_nonzero - col_end);
Chris@49:       }
Chris@49: 
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49: 
Chris@49:     access::rw(values)      = new_values;
Chris@49:     access::rw(row_indices) = new_row_indices;
Chris@49: 
Chris@49:     // Update counts and such.
Chris@49:     access::rw(n_nonzero) -= diff;
Chris@49:     }
Chris@49:   
Chris@49:   // Update column pointers.
Chris@49:   const uword new_n_cols = n_cols - ((in_col2 - in_col1) + 1);
Chris@49:   
Chris@49:   uword* new_col_ptrs = memory::acquire<uword>(new_n_cols + 2);
Chris@49:   new_col_ptrs[new_n_cols + 1] = std::numeric_limits<uword>::max();
Chris@49:   
Chris@49:   // Copy first set of columns (no manipulation required).
Chris@49:   if (in_col1 != 0)
Chris@49:     {
Chris@49:     arrayops::copy(new_col_ptrs, col_ptrs, in_col1);
Chris@49:     }
Chris@49:   
Chris@49:   // Copy second set of columns (manipulation required).
Chris@49:   uword cur_col = in_col1;
Chris@49:   for (uword i = in_col2 + 1; i <= n_cols; ++i, ++cur_col)
Chris@49:     {
Chris@49:     new_col_ptrs[cur_col] = col_ptrs[i] - diff;
Chris@49:     }
Chris@49:   
Chris@49:   memory::release(col_ptrs);
Chris@49:   access::rw(col_ptrs) = new_col_ptrs;
Chris@49:   
Chris@49:   // We update the element and column counts, and we're done.
Chris@49:   access::rw(n_cols) = new_n_cols;
Chris@49:   access::rw(n_elem) = n_cols * n_rows;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::rows(const uword in_row1, const uword in_row2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 > in_row2) || (in_row2 >= n_rows),
Chris@49:     "SpMat::rows(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_rows = in_row2 - in_row1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, in_row1, 0, subview_n_rows, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::rows(const uword in_row1, const uword in_row2) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 > in_row2) || (in_row2 >= n_rows),
Chris@49:     "SpMat::rows(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_rows = in_row2 - in_row1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, in_row1, 0, subview_n_rows, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::cols(const uword in_col1, const uword in_col2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_col1 > in_col2) || (in_col2 >= n_cols),
Chris@49:     "SpMat::cols(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_cols = in_col2 - in_col1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, 0, in_col1, n_rows, subview_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::cols(const uword in_col1, const uword in_col2) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_col1 > in_col2) || (in_col2 >= n_cols),
Chris@49:     "SpMat::cols(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_cols = in_col2 - in_col1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, 0, in_col1, n_rows, subview_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::submat(const uword in_row1, const uword in_col1, const uword in_row2, const uword in_col2)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 > in_row2) || (in_col1 >  in_col2) || (in_row2 >= n_rows) || (in_col2 >= n_cols),
Chris@49:     "SpMat::submat(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_rows = in_row2 - in_row1 + 1;
Chris@49:   const uword subview_n_cols = in_col2 - in_col1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, in_row1, in_col1, subview_n_rows, subview_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::submat(const uword in_row1, const uword in_col1, const uword in_row2, const uword in_col2) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (in_row1 > in_row2) || (in_col1 >  in_col2) || (in_row2 >= n_rows) || (in_col2 >= n_cols),
Chris@49:     "SpMat::submat(): indices out of bounds or incorrectly used"
Chris@49:     );
Chris@49: 
Chris@49:   const uword subview_n_rows = in_row2 - in_row1 + 1;
Chris@49:   const uword subview_n_cols = in_col2 - in_col1 + 1;
Chris@49: 
Chris@49:   return SpSubview<eT>(*this, in_row1, in_col1, subview_n_rows, subview_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::submat    (const span& row_span, const span& col_span)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const bool row_all = row_span.whole;
Chris@49:   const bool col_all = col_span.whole;
Chris@49:   
Chris@49:   const uword local_n_rows = n_rows;
Chris@49:   const uword local_n_cols = n_cols;
Chris@49:   
Chris@49:   const uword in_row1       = row_all ? 0            : row_span.a;
Chris@49:   const uword in_row2       =                          row_span.b;
Chris@49:   const uword submat_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; 
Chris@49:   
Chris@49:   const uword in_col1       = col_all ? 0            : col_span.a;
Chris@49:   const uword in_col2       =                          col_span.b;
Chris@49:   const uword submat_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1; 
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (    
Chris@49:     ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) )
Chris@49:     ||   
Chris@49:     ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) )
Chris@49:     ,    
Chris@49:     "SpMat::submat(): indices out of bounds or incorrectly used"
Chris@49:     );   
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, in_row1, in_col1, submat_n_rows, submat_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::submat    (const span& row_span, const span& col_span) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const bool row_all = row_span.whole;
Chris@49:   const bool col_all = col_span.whole;
Chris@49:   
Chris@49:   const uword local_n_rows = n_rows;
Chris@49:   const uword local_n_cols = n_cols;
Chris@49:   
Chris@49:   const uword in_row1       = row_all ? 0            : row_span.a;
Chris@49:   const uword in_row2       =                          row_span.b;
Chris@49:   const uword submat_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; 
Chris@49:   
Chris@49:   const uword in_col1       = col_all ? 0            : col_span.a;
Chris@49:   const uword in_col2       =                          col_span.b;
Chris@49:   const uword submat_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1; 
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (    
Chris@49:     ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) )
Chris@49:     ||   
Chris@49:     ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) )
Chris@49:     ,    
Chris@49:     "SpMat::submat(): indices out of bounds or incorrectly used"
Chris@49:     );   
Chris@49:   
Chris@49:   return SpSubview<eT>(*this, in_row1, in_col1, submat_n_rows, submat_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const span& row_span, const span& col_span)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   return submat(row_span, col_span);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpSubview<eT>
Chris@49: SpMat<eT>::operator()(const span& row_span, const span& col_span) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   return submat(row_span, col_span);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Element access; acces the i'th element (works identically to the Mat accessors).
Chris@49:  * If there is nothing at element i, 0 is returned.
Chris@49:  *
Chris@49:  * @param i Element to access.
Chris@49:  */
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::operator[](const uword i)
Chris@49:   {
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::operator[](const uword i) const
Chris@49:   {
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::at(const uword i)
Chris@49:   {
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::at(const uword i) const
Chris@49:   {
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::operator()(const uword i)
Chris@49:   {
Chris@49:   arma_debug_check( (i >= n_elem), "SpMat::operator(): out of bounds");
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::operator()(const uword i) const
Chris@49:   {
Chris@49:   arma_debug_check( (i >= n_elem), "SpMat::operator(): out of bounds");
Chris@49:   return get_value(i);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Element access; access the element at row in_rows and column in_col.
Chris@49:  * If there is nothing at that position, 0 is returned.
Chris@49:  */
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::at(const uword in_row, const uword in_col)
Chris@49:   {
Chris@49:   return get_value(in_row, in_col);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::at(const uword in_row, const uword in_col) const
Chris@49:   {
Chris@49:   return get_value(in_row, in_col);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::operator()(const uword in_row, const uword in_col)
Chris@49:   {
Chris@49:   arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "SpMat::operator(): out of bounds");
Chris@49:   return get_value(in_row, in_col);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::operator()(const uword in_row, const uword in_col) const
Chris@49:   {
Chris@49:   arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "SpMat::operator(): out of bounds");
Chris@49:   return get_value(in_row, in_col);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Check if matrix is empty (no size, no values).
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_empty() const
Chris@49:   {
Chris@49:   return(n_elem == 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the object can be interpreted as a column or row vector
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_vec() const
Chris@49:   {
Chris@49:   return ( (n_rows == 1) || (n_cols == 1) );
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the object can be interpreted as a row vector
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_rowvec() const
Chris@49:   {
Chris@49:   return (n_rows == 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the object can be interpreted as a column vector
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_colvec() const
Chris@49:   {
Chris@49:   return (n_cols == 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the object has the same number of non-zero rows and columnns
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_square() const
Chris@49:   {
Chris@49:   return (n_rows == n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if all of the elements are finite
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::is_finite() const
Chris@49:   {
Chris@49:   for(uword i = 0; i < n_nonzero; i++)
Chris@49:     {
Chris@49:     if(arma_isfinite(values[i]) == false)
Chris@49:       {
Chris@49:       return false;
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   return true; // No infinite values.
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the given index is currently in range
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const uword i) const
Chris@49:   {
Chris@49:   return (i < n_elem);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the given start and end indices are currently in range
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const span& x) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(x.whole == true)
Chris@49:     {
Chris@49:     return true;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword a = x.a;
Chris@49:     const uword b = x.b;
Chris@49: 
Chris@49:     return ( (a <= b) && (b < n_elem) );
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! returns true if the given location is currently in range
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const uword in_row, const uword in_col) const
Chris@49:   {
Chris@49:   return ( (in_row < n_rows) && (in_col < n_cols) );
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const span& row_span, const uword in_col) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(row_span.whole == true)
Chris@49:     {
Chris@49:     return (in_col < n_cols);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword in_row1 = row_span.a;
Chris@49:     const uword in_row2 = row_span.b;
Chris@49: 
Chris@49:     return ( (in_row1 <= in_row2) && (in_row2 < n_rows) && (in_col < n_cols) );
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const uword in_row, const span& col_span) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(col_span.whole == true)
Chris@49:     {
Chris@49:     return (in_row < n_rows);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword in_col1 = col_span.a;
Chris@49:     const uword in_col2 = col_span.b;
Chris@49: 
Chris@49:     return ( (in_row < n_rows) && (in_col1 <= in_col2) && (in_col2 < n_cols) );
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: bool
Chris@49: SpMat<eT>::in_range(const span& row_span, const span& col_span) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   const uword in_row1 = row_span.a;
Chris@49:   const uword in_row2 = row_span.b;
Chris@49: 
Chris@49:   const uword in_col1 = col_span.a;
Chris@49:   const uword in_col2 = col_span.b;
Chris@49: 
Chris@49:   const bool rows_ok = row_span.whole ? true : ( (in_row1 <= in_row2) && (in_row2 < n_rows) );
Chris@49:   const bool cols_ok = col_span.whole ? true : ( (in_col1 <= in_col2) && (in_col2 < n_cols) );
Chris@49: 
Chris@49:   return ( (rows_ok == true) && (cols_ok == true) );
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_print(const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = ARMA_DEFAULT_OSTREAM.width();
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM << extra_text << '\n';
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print(ARMA_DEFAULT_OSTREAM, *this, true);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_print(std::ostream& user_stream, const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = user_stream.width();
Chris@49: 
Chris@49:     user_stream << extra_text << '\n';
Chris@49: 
Chris@49:     user_stream.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print(user_stream, *this, true);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_raw_print(const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = ARMA_DEFAULT_OSTREAM.width();
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM << extra_text << '\n';
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print(ARMA_DEFAULT_OSTREAM, *this, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_raw_print(std::ostream& user_stream, const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = user_stream.width();
Chris@49: 
Chris@49:     user_stream << extra_text << '\n';
Chris@49: 
Chris@49:     user_stream.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print(user_stream, *this, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Matrix printing, prepends supplied text.
Chris@49:  * Prints 0 wherever no element exists.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_print_dense(const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = ARMA_DEFAULT_OSTREAM.width();
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM << extra_text << '\n';
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print_dense(ARMA_DEFAULT_OSTREAM, *this, true);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_print_dense(std::ostream& user_stream, const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = user_stream.width();
Chris@49: 
Chris@49:     user_stream << extra_text << '\n';
Chris@49: 
Chris@49:     user_stream.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print_dense(user_stream, *this, true);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_raw_print_dense(const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = ARMA_DEFAULT_OSTREAM.width();
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM << extra_text << '\n';
Chris@49: 
Chris@49:     ARMA_DEFAULT_OSTREAM.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print_dense(ARMA_DEFAULT_OSTREAM, *this, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::impl_raw_print_dense(std::ostream& user_stream, const std::string& extra_text) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(extra_text.length() != 0)
Chris@49:     {
Chris@49:     const std::streamsize orig_width = user_stream.width();
Chris@49: 
Chris@49:     user_stream << extra_text << '\n';
Chris@49: 
Chris@49:     user_stream.width(orig_width);
Chris@49:     }
Chris@49: 
Chris@49:   arma_ostream::print_dense(user_stream, *this, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Set the size to the size of another matrix.
Chris@49: template<typename eT>
Chris@49: template<typename eT2>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::copy_size(const SpMat<eT2>& m)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(m.n_rows, m.n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename eT2>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::copy_size(const Mat<eT2>& m)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(m.n_rows, m.n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Resize the matrix to a given size.  The matrix will be resized to be a column vector (i.e. in_elem columns, 1 row).
Chris@49:  *
Chris@49:  * @param in_elem Number of elements to allow.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::set_size(const uword in_elem)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   // If this is a row vector, we resize to a row vector.
Chris@49:   if(vec_state == 2)
Chris@49:     {
Chris@49:     init(1, in_elem);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     init(in_elem, 1);
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Resize the matrix to a given size.
Chris@49:  *
Chris@49:  * @param in_rows Number of rows to allow.
Chris@49:  * @param in_cols Number of columns to allow.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::set_size(const uword in_rows, const uword in_cols)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(in_rows, in_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::reshape(const uword in_rows, const uword in_cols, const uword dim)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if (dim == 0)
Chris@49:     {
Chris@49:     // We have to modify all of the relevant row indices and the relevant column pointers.
Chris@49:     // Iterate over all the points to do this.  We won't be deleting any points, but we will be modifying
Chris@49:     // columns and rows. We'll have to store a new set of column vectors.
Chris@49:     uword* new_col_ptrs    = memory::acquire<uword>(in_cols + 2);
Chris@49:     new_col_ptrs[in_cols + 1] = std::numeric_limits<uword>::max();
Chris@49:     
Chris@49:     uword* new_row_indices = memory::acquire_chunked<uword>(n_nonzero + 1);
Chris@49:     access::rw(new_row_indices[n_nonzero]) = 0;
Chris@49: 
Chris@49:     arrayops::inplace_set(new_col_ptrs, uword(0), in_cols + 1);
Chris@49: 
Chris@49:     for(const_iterator it = begin(); it != end(); it++)
Chris@49:       {
Chris@49:       uword vector_position = (it.col() * n_rows) + it.row();
Chris@49:       new_row_indices[it.pos()] = vector_position % in_rows;
Chris@49:       ++new_col_ptrs[vector_position / in_rows + 1];
Chris@49:       }
Chris@49: 
Chris@49:     // Now sum the column counts to get the new column pointers.
Chris@49:     for(uword i = 1; i <= in_cols; i++)
Chris@49:       {
Chris@49:       access::rw(new_col_ptrs[i]) += new_col_ptrs[i - 1];
Chris@49:       }
Chris@49: 
Chris@49:     // Copy the new row indices.
Chris@49:     memory::release(row_indices);
Chris@49:     access::rw(row_indices) = new_row_indices;
Chris@49: 
Chris@49:     memory::release(col_ptrs);
Chris@49:     access::rw(col_ptrs) = new_col_ptrs;
Chris@49: 
Chris@49:     // Now set the size.
Chris@49:     access::rw(n_rows) = in_rows;
Chris@49:     access::rw(n_cols) = in_cols;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // Row-wise reshaping.  This is more tedious and we will use a separate sparse matrix to do it.
Chris@49:     SpMat<eT> tmp(in_rows, in_cols);
Chris@49: 
Chris@49:     for(const_row_iterator it = begin_row(); it.pos() < n_nonzero; it++)
Chris@49:       {
Chris@49:       uword vector_position = (it.row() * n_cols) + it.col();
Chris@49: 
Chris@49:       tmp((vector_position / in_cols), (vector_position % in_cols)) = (*it);
Chris@49:       }
Chris@49: 
Chris@49:     (*this).operator=(tmp);
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::zeros()
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if (n_nonzero > 0)
Chris@49:     {
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49: 
Chris@49:     access::rw(values)      = memory::acquire_chunked<eT>(1);
Chris@49:     access::rw(row_indices) = memory::acquire_chunked<uword>(1);
Chris@49: 
Chris@49:     access::rw(values[0]) = 0;
Chris@49:     access::rw(row_indices[0]) = 0;
Chris@49:     }
Chris@49: 
Chris@49:   access::rw(n_nonzero) = 0;
Chris@49:   arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1);
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::zeros(const uword in_elem)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(vec_state == 2)
Chris@49:     {
Chris@49:     init(1, in_elem); // Row vector
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     init(in_elem, 1);
Chris@49:     }
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::zeros(const uword in_rows, const uword in_cols)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   init(in_rows, in_cols);
Chris@49: 
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::eye()
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   return (*this).eye(n_rows, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::eye(const uword in_rows, const uword in_cols)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword N = (std::min)(in_rows, in_cols);
Chris@49:   
Chris@49:   init(in_rows, in_cols);
Chris@49:   
Chris@49:   mem_resize(N);
Chris@49:   
Chris@49:   arrayops::inplace_set(access::rwp(values), eT(1), N);
Chris@49:   
Chris@49:   for(uword i = 0; i <  N; ++i) { access::rw(row_indices[i]) = i; }
Chris@49:   
Chris@49:   for(uword i = 0; i <= N; ++i) { access::rw(col_ptrs[i])    = i; }
Chris@49:   
Chris@49:   access::rw(n_nonzero) = N;
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::speye()
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   return (*this).eye(n_rows, n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::speye(const uword in_n_rows, const uword in_n_cols)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   return (*this).eye(in_n_rows, in_n_cols);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::sprandu(const uword in_rows, const uword in_cols, const double density)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_check( ( (density < double(0)) || (density > double(1)) ), "sprandu(): density must be in the [0,1] interval" );
Chris@49:   
Chris@49:   zeros(in_rows, in_cols);
Chris@49:   
Chris@49:   mem_resize( uword(density * double(in_rows) * double(in_cols) + 0.5) );
Chris@49:   
Chris@49:   if(n_nonzero == 0)
Chris@49:     {
Chris@49:     return *this;
Chris@49:     }
Chris@49:   
Chris@49:   eop_aux_randu<eT>::fill( access::rwp(values), n_nonzero );
Chris@49:   
Chris@49:   uvec indices = linspace<uvec>( 0u, in_rows*in_cols-1, n_nonzero );
Chris@49:   
Chris@49:   // perturb the indices
Chris@49:   for(uword i=1; i < n_nonzero-1; ++i)
Chris@49:     {
Chris@49:     const uword index_left  = indices[i-1];
Chris@49:     const uword index_right = indices[i+1];
Chris@49:     
Chris@49:     const uword center = (index_left + index_right) / 2;
Chris@49:     
Chris@49:     const uword delta1 = center      - index_left - 1;
Chris@49:     const uword delta2 = index_right - center     - 1;
Chris@49:     
Chris@49:     const uword min_delta = (std::min)(delta1, delta2);
Chris@49:     
Chris@49:     uword index_new = uword( double(center) + double(min_delta) * (2.0*randu()-1.0) );
Chris@49:     
Chris@49:     // paranoia, but better be safe than sorry
Chris@49:     if( (index_left < index_new) && (index_new < index_right) )
Chris@49:       {
Chris@49:       indices[i] = index_new;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   uword cur_index = 0;
Chris@49:   uword count     = 0;  
Chris@49:   
Chris@49:   for(uword lcol = 0; lcol < in_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < in_rows; ++lrow)
Chris@49:     {
Chris@49:     if(count == indices[cur_index])
Chris@49:       {
Chris@49:       access::rw(row_indices[cur_index]) = lrow;
Chris@49:       access::rw(col_ptrs[lcol + 1])++;
Chris@49:       ++cur_index;
Chris@49:       }
Chris@49:     
Chris@49:     ++count;
Chris@49:     }
Chris@49:   
Chris@49:   if(cur_index != n_nonzero)
Chris@49:     {
Chris@49:     // Fix size to correct size.
Chris@49:     mem_resize(cur_index);
Chris@49:     }
Chris@49:   
Chris@49:   // Sum column pointers.
Chris@49:   for(uword lcol = 1; lcol <= in_cols; ++lcol)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[lcol]) += col_ptrs[lcol - 1];
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: const SpMat<eT>&
Chris@49: SpMat<eT>::sprandn(const uword in_rows, const uword in_cols, const double density)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_check( ( (density < double(0)) || (density > double(1)) ), "sprandn(): density must be in the [0,1] interval" );
Chris@49:   
Chris@49:   zeros(in_rows, in_cols);
Chris@49:   
Chris@49:   mem_resize( uword(density * double(in_rows) * double(in_cols) + 0.5) );
Chris@49:   
Chris@49:   if(n_nonzero == 0)
Chris@49:     {
Chris@49:     return *this;
Chris@49:     }
Chris@49:   
Chris@49:   eop_aux_randn<eT>::fill( access::rwp(values), n_nonzero );
Chris@49:   
Chris@49:   uvec indices = linspace<uvec>( 0u, in_rows*in_cols-1, n_nonzero );
Chris@49:   
Chris@49:   // perturb the indices
Chris@49:   for(uword i=1; i < n_nonzero-1; ++i)
Chris@49:     {
Chris@49:     const uword index_left  = indices[i-1];
Chris@49:     const uword index_right = indices[i+1];
Chris@49:     
Chris@49:     const uword center = (index_left + index_right) / 2;
Chris@49:     
Chris@49:     const uword delta1 = center      - index_left - 1;
Chris@49:     const uword delta2 = index_right - center     - 1;
Chris@49:     
Chris@49:     const uword min_delta = (std::min)(delta1, delta2);
Chris@49:     
Chris@49:     uword index_new = uword( double(center) + double(min_delta) * (2.0*randu()-1.0) );
Chris@49:     
Chris@49:     // paranoia, but better be safe than sorry
Chris@49:     if( (index_left < index_new) && (index_new < index_right) )
Chris@49:       {
Chris@49:       indices[i] = index_new;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   uword cur_index = 0;
Chris@49:   uword count     = 0;  
Chris@49:   
Chris@49:   for(uword lcol = 0; lcol < in_cols; ++lcol)
Chris@49:   for(uword lrow = 0; lrow < in_rows; ++lrow)
Chris@49:     {
Chris@49:     if(count == indices[cur_index])
Chris@49:       {
Chris@49:       access::rw(row_indices[cur_index]) = lrow;
Chris@49:       access::rw(col_ptrs[lcol + 1])++;
Chris@49:       ++cur_index;
Chris@49:       }
Chris@49:     
Chris@49:     ++count;
Chris@49:     }
Chris@49:   
Chris@49:   if(cur_index != n_nonzero)
Chris@49:     {
Chris@49:     // Fix size to correct size.
Chris@49:     mem_resize(cur_index);
Chris@49:     }
Chris@49:   
Chris@49:   // Sum column pointers.
Chris@49:   for(uword lcol = 1; lcol <= in_cols; ++lcol)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[lcol]) += col_ptrs[lcol - 1];
Chris@49:     }
Chris@49:   
Chris@49:   return *this;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::reset()
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   set_size(0, 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Get the minimum or the maximum of the matrix.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::min() const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "min(): object has no elements");
Chris@49: 
Chris@49:   if (n_nonzero == 0)
Chris@49:     {
Chris@49:     return 0;
Chris@49:     }
Chris@49: 
Chris@49:   eT val = op_min::direct_min(values, n_nonzero);
Chris@49: 
Chris@49:   if ((val > 0) && (n_nonzero < n_elem)) // A sparse 0 is less.
Chris@49:     {
Chris@49:     val = 0;
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: eT
Chris@49: SpMat<eT>::min(uword& index_of_min_val) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "min(): object has no elements");
Chris@49: 
Chris@49:   eT val = 0;
Chris@49: 
Chris@49:   if (n_nonzero == 0) // There are no other elements.  It must be 0.
Chris@49:     {
Chris@49:     index_of_min_val = 0;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     uword location;
Chris@49:     val = op_min::direct_min(values, n_nonzero, location);
Chris@49: 
Chris@49:     if ((val > 0) && (n_nonzero < n_elem)) // A sparse 0 is less.
Chris@49:       {
Chris@49:       val = 0;
Chris@49: 
Chris@49:       // Give back the index to the first zero position.
Chris@49:       index_of_min_val = 0;
Chris@49:       while (get_position(index_of_min_val) == index_of_min_val) // An element exists at that position.
Chris@49:         {
Chris@49:         index_of_min_val++;
Chris@49:         }
Chris@49: 
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       index_of_min_val = get_position(location);
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: eT
Chris@49: SpMat<eT>::min(uword& row_of_min_val, uword& col_of_min_val) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "min(): object has no elements");
Chris@49: 
Chris@49:   eT val = 0;
Chris@49: 
Chris@49:   if (n_nonzero == 0) // There are no other elements.  It must be 0.
Chris@49:     {
Chris@49:     row_of_min_val = 0;
Chris@49:     col_of_min_val = 0;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     uword location;
Chris@49:     val = op_min::direct_min(values, n_nonzero, location);
Chris@49: 
Chris@49:     if ((val > 0) && (n_nonzero < n_elem)) // A sparse 0 is less.
Chris@49:       {
Chris@49:       val = 0;
Chris@49: 
Chris@49:       location = 0;
Chris@49:       while (get_position(location) == location) // An element exists at that position.
Chris@49:         {
Chris@49:         location++;
Chris@49:         }
Chris@49: 
Chris@49:       row_of_min_val = location % n_rows;
Chris@49:       col_of_min_val = location / n_rows;
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       get_position(location, row_of_min_val, col_of_min_val);
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::max() const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "max(): object has no elements");
Chris@49: 
Chris@49:   if (n_nonzero == 0)
Chris@49:     {
Chris@49:     return 0;
Chris@49:     }
Chris@49: 
Chris@49:   eT val = op_max::direct_max(values, n_nonzero);
Chris@49: 
Chris@49:   if ((val < 0) && (n_nonzero < n_elem)) // A sparse 0 is more.
Chris@49:     {
Chris@49:     return 0;
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: eT
Chris@49: SpMat<eT>::max(uword& index_of_max_val) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "max(): object has no elements");
Chris@49: 
Chris@49:   eT val = 0;
Chris@49: 
Chris@49:   if (n_nonzero == 0)
Chris@49:     {
Chris@49:     index_of_max_val = 0;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     uword location;
Chris@49:     val = op_max::direct_max(values, n_nonzero, location);
Chris@49: 
Chris@49:     if ((val < 0) && (n_nonzero < n_elem)) // A sparse 0 is more.
Chris@49:       {
Chris@49:       val = 0;
Chris@49: 
Chris@49:       location = 0;
Chris@49:       while (get_position(location) == location) // An element exists at that position.
Chris@49:         {
Chris@49:         location++;
Chris@49:         }
Chris@49: 
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       index_of_max_val = get_position(location);
Chris@49:       }
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: eT
Chris@49: SpMat<eT>::max(uword& row_of_max_val, uword& col_of_max_val) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   arma_debug_check((n_elem == 0), "max(): object has no elements");
Chris@49: 
Chris@49:   eT val = 0;
Chris@49: 
Chris@49:   if (n_nonzero == 0)
Chris@49:     {
Chris@49:     row_of_max_val = 0;
Chris@49:     col_of_max_val = 0;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     uword location;
Chris@49:     val = op_max::direct_max(values, n_nonzero, location);
Chris@49: 
Chris@49:     if ((val < 0) && (n_nonzero < n_elem)) // A sparse 0 is more.
Chris@49:       {
Chris@49:       val = 0;
Chris@49: 
Chris@49:       location = 0;
Chris@49:       while (get_position(location) == location) // An element exists at that position.
Chris@49:         {
Chris@49:         location++;
Chris@49:         }
Chris@49: 
Chris@49:       row_of_max_val = location % n_rows;
Chris@49:       col_of_max_val = location / n_rows;
Chris@49: 
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       get_position(location, row_of_max_val, col_of_max_val);
Chris@49:       }
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   return val;
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! save the matrix to a file
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::save(const std::string name, const file_type type, const bool print_status) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   bool save_okay;
Chris@49:   
Chris@49:   switch(type)
Chris@49:     {
Chris@49:     // case raw_ascii:
Chris@49:     //   save_okay = diskio::save_raw_ascii(*this, name);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case csv_ascii:
Chris@49:     //   save_okay = diskio::save_csv_ascii(*this, name);
Chris@49:     //   break;
Chris@49:     
Chris@49:     case arma_binary:
Chris@49:       save_okay = diskio::save_arma_binary(*this, name);
Chris@49:       break;
Chris@49:     
Chris@49:     case coord_ascii:
Chris@49:       save_okay = diskio::save_coord_ascii(*this, name);
Chris@49:       break;
Chris@49:     
Chris@49:     default:
Chris@49:       arma_warn(true, "SpMat::save(): unsupported file type");
Chris@49:       save_okay = false;
Chris@49:     }
Chris@49:   
Chris@49:   arma_warn( (save_okay == false), "SpMat::save(): couldn't write to ", name);
Chris@49:   
Chris@49:   return save_okay;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! save the matrix to a stream
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::save(std::ostream& os, const file_type type, const bool print_status) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   bool save_okay;
Chris@49:   
Chris@49:   switch(type)
Chris@49:     {
Chris@49:     // case raw_ascii:
Chris@49:     //   save_okay = diskio::save_raw_ascii(*this, os);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case csv_ascii:
Chris@49:     //   save_okay = diskio::save_csv_ascii(*this, os);
Chris@49:     //   break;
Chris@49:     
Chris@49:     case arma_binary:
Chris@49:       save_okay = diskio::save_arma_binary(*this, os);
Chris@49:       break;
Chris@49:     
Chris@49:     case coord_ascii:
Chris@49:       save_okay = diskio::save_coord_ascii(*this, os);
Chris@49:       break;
Chris@49:     
Chris@49:     default:
Chris@49:       arma_warn(true, "SpMat::save(): unsupported file type");
Chris@49:       save_okay = false;
Chris@49:     }
Chris@49:   
Chris@49:   arma_warn( (save_okay == false), "SpMat::save(): couldn't write to the given stream");
Chris@49:   
Chris@49:   return save_okay;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! load a matrix from a file
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::load(const std::string name, const file_type type, const bool print_status)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   bool load_okay;
Chris@49:   std::string err_msg;
Chris@49:   
Chris@49:   switch(type)
Chris@49:     {
Chris@49:     // case auto_detect:
Chris@49:     //   load_okay = diskio::load_auto_detect(*this, name, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case raw_ascii:
Chris@49:     //   load_okay = diskio::load_raw_ascii(*this, name, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case csv_ascii:
Chris@49:     //   load_okay = diskio::load_csv_ascii(*this, name, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     case arma_binary:
Chris@49:       load_okay = diskio::load_arma_binary(*this, name, err_msg);
Chris@49:       break;
Chris@49:     
Chris@49:     case coord_ascii:
Chris@49:       load_okay = diskio::load_coord_ascii(*this, name, err_msg);
Chris@49:       break;
Chris@49:     
Chris@49:     default:
Chris@49:       arma_warn(true, "SpMat::load(): unsupported file type");
Chris@49:       load_okay = false;
Chris@49:     }
Chris@49:   
Chris@49:   if(load_okay == false)
Chris@49:     {
Chris@49:     if(err_msg.length() > 0)
Chris@49:       {
Chris@49:       arma_warn(true, "SpMat::load(): ", err_msg, name);
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       arma_warn(true, "SpMat::load(): couldn't read ", name);
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   if(load_okay == false)
Chris@49:     {
Chris@49:     (*this).reset();
Chris@49:     }
Chris@49:     
Chris@49:   return load_okay;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! load a matrix from a stream
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::load(std::istream& is, const file_type type, const bool print_status)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   bool load_okay;
Chris@49:   std::string err_msg;
Chris@49:   
Chris@49:   switch(type)
Chris@49:     {
Chris@49:     // case auto_detect:
Chris@49:     //   load_okay = diskio::load_auto_detect(*this, is, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case raw_ascii:
Chris@49:     //   load_okay = diskio::load_raw_ascii(*this, is, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     // case csv_ascii:
Chris@49:     //   load_okay = diskio::load_csv_ascii(*this, is, err_msg);
Chris@49:     //   break;
Chris@49:     
Chris@49:     case arma_binary:
Chris@49:       load_okay = diskio::load_arma_binary(*this, is, err_msg);
Chris@49:       break;
Chris@49:     
Chris@49:     case coord_ascii:
Chris@49:       load_okay = diskio::load_coord_ascii(*this, is, err_msg);
Chris@49:       break;
Chris@49:     
Chris@49:     default:
Chris@49:       arma_warn(true, "SpMat::load(): unsupported file type");
Chris@49:       load_okay = false;
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   if(load_okay == false)
Chris@49:     {
Chris@49:     if(err_msg.length() > 0)
Chris@49:       {
Chris@49:       arma_warn(true, "SpMat::load(): ", err_msg, "the given stream");
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       arma_warn(true, "SpMat::load(): couldn't load from the given stream");
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   if(load_okay == false)
Chris@49:     {
Chris@49:     (*this).reset();
Chris@49:     }
Chris@49:     
Chris@49:   return load_okay;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! save the matrix to a file, without printing any error messages
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::quiet_save(const std::string name, const file_type type) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   return (*this).save(name, type, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! save the matrix to a stream, without printing any error messages
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::quiet_save(std::ostream& os, const file_type type) const
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   return (*this).save(os, type, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! load a matrix from a file, without printing any error messages
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::quiet_load(const std::string name, const file_type type)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   return (*this).load(name, type, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! load a matrix from a stream, without printing any error messages
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::quiet_load(std::istream& is, const file_type type)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   return (*this).load(is, type, false);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Initialize the matrix to the specified size.  Data is not preserved, so the matrix is assumed to be entirely sparse (empty).
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::init(uword in_rows, uword in_cols)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   // Verify that we are allowed to do this.
Chris@49:   if(vec_state > 0)
Chris@49:     {
Chris@49:     if((in_rows == 0) && (in_cols == 0))
Chris@49:       {
Chris@49:       if(vec_state == 1)
Chris@49:         {
Chris@49:         in_cols = 1;
Chris@49:         }
Chris@49:       else
Chris@49:       if(vec_state == 2)
Chris@49:         {
Chris@49:         in_rows = 1;
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       arma_debug_check
Chris@49:         (
Chris@49:         ( ((vec_state == 1) && (in_cols != 1)) || ((vec_state == 2) && (in_rows != 1)) ),
Chris@49:         "SpMat::init(): object is a row or column vector; requested size is not compatible"
Chris@49:         );
Chris@49:       }
Chris@49:     }
Chris@49: 
Chris@49:   // Ensure that n_elem can hold the result of (n_rows * n_cols)
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:       (
Chris@49:       ( (in_rows > ARMA_MAX_UHWORD) || (in_cols > ARMA_MAX_UHWORD) )
Chris@49:         ? ( (float(in_rows) * float(in_cols)) > float(ARMA_MAX_UWORD) )
Chris@49:         : false
Chris@49:       ),
Chris@49:       "SpMat::init(): requested size is too large"
Chris@49:     );
Chris@49: 
Chris@49:   // Clean out the existing memory.
Chris@49:   if (values)
Chris@49:     {
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49:     }
Chris@49: 
Chris@49:   access::rw(values)      = memory::acquire_chunked<eT>   (1);
Chris@49:   access::rw(row_indices) = memory::acquire_chunked<uword>(1);
Chris@49: 
Chris@49:   access::rw(values[0]) = 0;
Chris@49:   access::rw(row_indices[0]) = 0;
Chris@49: 
Chris@49:   memory::release(col_ptrs);
Chris@49: 
Chris@49:   // Set the new size accordingly.
Chris@49:   access::rw(n_rows)    = in_rows;
Chris@49:   access::rw(n_cols)    = in_cols;
Chris@49:   access::rw(n_elem)    = (in_rows * in_cols);
Chris@49:   access::rw(n_nonzero) = 0;
Chris@49: 
Chris@49:   // Try to allocate the column pointers, filling them with 0, except for the
Chris@49:   // last element which contains the maximum possible element (so iterators
Chris@49:   // terminate correctly).
Chris@49:   access::rw(col_ptrs) = memory::acquire<uword>(in_cols + 2);
Chris@49:   access::rw(col_ptrs[in_cols + 1]) = std::numeric_limits<uword>::max();
Chris@49:   
Chris@49:   arrayops::inplace_set(access::rwp(col_ptrs), uword(0), in_cols + 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Initialize the matrix from a string.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::init(const std::string& text)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   // Figure out the size first.
Chris@49:   uword t_n_rows = 0;
Chris@49:   uword t_n_cols = 0;
Chris@49: 
Chris@49:   bool t_n_cols_found = false;
Chris@49: 
Chris@49:   std::string token;
Chris@49: 
Chris@49:   std::string::size_type line_start = 0;
Chris@49:   std::string::size_type   line_end = 0;
Chris@49: 
Chris@49:   while (line_start < text.length())
Chris@49:     {
Chris@49: 
Chris@49:     line_end = text.find(';', line_start);
Chris@49: 
Chris@49:     if (line_end == std::string::npos)
Chris@49:       line_end = text.length() - 1;
Chris@49: 
Chris@49:     std::string::size_type line_len = line_end - line_start + 1;
Chris@49:     std::stringstream line_stream(text.substr(line_start, line_len));
Chris@49: 
Chris@49:     // Step through each column.
Chris@49:     uword line_n_cols = 0;
Chris@49: 
Chris@49:     while (line_stream >> token)
Chris@49:       {
Chris@49:       ++line_n_cols;
Chris@49:       }
Chris@49: 
Chris@49:     if (line_n_cols > 0)
Chris@49:       {
Chris@49:       if (t_n_cols_found == false)
Chris@49:         {
Chris@49:         t_n_cols = line_n_cols;
Chris@49:         t_n_cols_found = true;
Chris@49:         }
Chris@49:       else // Check it each time through, just to make sure.
Chris@49:         arma_check((line_n_cols != t_n_cols), "SpMat::init(): inconsistent number of columns in given string");
Chris@49: 
Chris@49:       ++t_n_rows;
Chris@49:       }
Chris@49: 
Chris@49:     line_start = line_end + 1;
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   set_size(t_n_rows, t_n_cols);
Chris@49: 
Chris@49:   // Second time through will pick up all the values.
Chris@49:   line_start = 0;
Chris@49:   line_end = 0;
Chris@49: 
Chris@49:   uword lrow = 0;
Chris@49: 
Chris@49:   while (line_start < text.length())
Chris@49:     {
Chris@49: 
Chris@49:     line_end = text.find(';', line_start);
Chris@49: 
Chris@49:     if (line_end == std::string::npos)
Chris@49:       line_end = text.length() - 1;
Chris@49: 
Chris@49:     std::string::size_type line_len = line_end - line_start + 1;
Chris@49:     std::stringstream line_stream(text.substr(line_start, line_len));
Chris@49: 
Chris@49:     uword lcol = 0;
Chris@49:     eT val;
Chris@49: 
Chris@49:     while (line_stream >> val)
Chris@49:       {
Chris@49:       // Only add nonzero elements.
Chris@49:       if (val != eT(0))
Chris@49:         {
Chris@49:         get_value(lrow, lcol) = val;
Chris@49:         }
Chris@49: 
Chris@49:       ++lcol;
Chris@49:       }
Chris@49: 
Chris@49:     ++lrow;
Chris@49:     line_start = line_end + 1;
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   }
Chris@49: 
Chris@49: /**
Chris@49:  * Copy from another matrix.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::init(const SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   // Ensure we are not initializing to ourselves.
Chris@49:   if (this != &x)
Chris@49:     {
Chris@49:     init(x.n_rows, x.n_cols);
Chris@49: 
Chris@49:     // values and row_indices may not be null.
Chris@49:     if (values != NULL)
Chris@49:       {
Chris@49:       memory::release(values);
Chris@49:       memory::release(row_indices);
Chris@49:       }
Chris@49: 
Chris@49:     access::rw(values)      = memory::acquire_chunked<eT>   (x.n_nonzero + 1);
Chris@49:     access::rw(row_indices) = memory::acquire_chunked<uword>(x.n_nonzero + 1);
Chris@49: 
Chris@49:     // Now copy over the elements.
Chris@49:     arrayops::copy(access::rwp(values),      x.values,      x.n_nonzero + 1);
Chris@49:     arrayops::copy(access::rwp(row_indices), x.row_indices, x.n_nonzero + 1);
Chris@49:     arrayops::copy(access::rwp(col_ptrs),    x.col_ptrs,    x.n_cols + 1);
Chris@49:     
Chris@49:     access::rw(n_nonzero) = x.n_nonzero;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::mem_resize(const uword new_n_nonzero)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   if(n_nonzero != new_n_nonzero)
Chris@49:     {
Chris@49:     if(new_n_nonzero == 0)
Chris@49:       {
Chris@49:       memory::release(values);
Chris@49:       memory::release(row_indices);
Chris@49:       
Chris@49:       access::rw(values)      = memory::acquire_chunked<eT>   (1);
Chris@49:       access::rw(row_indices) = memory::acquire_chunked<uword>(1);
Chris@49: 
Chris@49:       access::rw(values[0]) = 0;
Chris@49:       access::rw(row_indices[0]) = 0;
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       // Figure out the actual amount of memory currently allocated
Chris@49:       // NOTE: this relies on memory::acquire_chunked() being used for the 'values' and 'row_indices' arrays
Chris@49:       const uword n_alloc = memory::enlarge_to_mult_of_chunksize(n_nonzero);
Chris@49:       
Chris@49:       if(n_alloc < new_n_nonzero)
Chris@49:         {
Chris@49:         eT*    new_values      = memory::acquire_chunked<eT>   (new_n_nonzero + 1);
Chris@49:         uword* new_row_indices = memory::acquire_chunked<uword>(new_n_nonzero + 1);
Chris@49:         
Chris@49:         if(n_nonzero > 0)
Chris@49:           {
Chris@49:           // Copy old elements.
Chris@49:           uword copy_len = std::min(n_nonzero, new_n_nonzero);
Chris@49:           
Chris@49:           arrayops::copy(new_values,      values,      copy_len);
Chris@49:           arrayops::copy(new_row_indices, row_indices, copy_len);
Chris@49:           }
Chris@49:         
Chris@49:         memory::release(values);
Chris@49:         memory::release(row_indices);
Chris@49:         
Chris@49:         access::rw(values)      = new_values;
Chris@49:         access::rw(row_indices) = new_row_indices;
Chris@49:         }
Chris@49:         
Chris@49:       // Set the "fake end" of the matrix by setting the last value and row
Chris@49:       // index to 0.  This helps the iterators work correctly.
Chris@49:       access::rw(values[new_n_nonzero]) = 0;
Chris@49:       access::rw(row_indices[new_n_nonzero]) = 0;
Chris@49:       }
Chris@49:     
Chris@49:     access::rw(n_nonzero) = new_n_nonzero;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: // Steal memory from another matrix.
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::steal_mem(SpMat<eT>& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49: 
Chris@49:   if(this != &x)
Chris@49:     {
Chris@49:     // Release all the memory.
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49:     memory::release(col_ptrs);
Chris@49: 
Chris@49:     // We'll have to copy everything about the other matrix.
Chris@49:     const uword x_n_rows    = x.n_rows;
Chris@49:     const uword x_n_cols    = x.n_cols;
Chris@49:     const uword x_n_elem    = x.n_elem;
Chris@49:     const uword x_n_nonzero = x.n_nonzero;
Chris@49: 
Chris@49:     access::rw(n_rows)    = x_n_rows;
Chris@49:     access::rw(n_cols)    = x_n_cols;
Chris@49:     access::rw(n_elem)    = x_n_elem;
Chris@49:     access::rw(n_nonzero) = x_n_nonzero;
Chris@49: 
Chris@49:     access::rw(values)      = x.values;
Chris@49:     access::rw(row_indices) = x.row_indices;
Chris@49:     access::rw(col_ptrs)    = x.col_ptrs;
Chris@49: 
Chris@49:     // Set other matrix to empty.
Chris@49:     access::rw(x.n_rows)    = 0;
Chris@49:     access::rw(x.n_cols)    = 0;
Chris@49:     access::rw(x.n_elem)    = 0;
Chris@49:     access::rw(x.n_nonzero) = 0;
Chris@49: 
Chris@49:     access::rw(x.values)      = NULL;
Chris@49:     access::rw(x.row_indices) = NULL;
Chris@49:     access::rw(x.col_ptrs)    = NULL;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename T1, typename Functor>
Chris@49: arma_hot
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::init_xform(const SpBase<eT,T1>& A, const Functor& func)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   // if possible, avoid doing a copy and instead apply func to the generated elements
Chris@49:   if(SpProxy<T1>::Q_created_by_proxy == true)
Chris@49:     {
Chris@49:     (*this) = A.get_ref();
Chris@49:     
Chris@49:     const uword nnz = n_nonzero;
Chris@49:     
Chris@49:     eT* t_values = access::rwp(values);
Chris@49:     
Chris@49:     for(uword i=0; i < nnz; ++i)
Chris@49:       {
Chris@49:       t_values[i] = func(t_values[i]);
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     init_xform_mt(A.get_ref(), func);
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: template<typename eT2, typename T1, typename Functor>
Chris@49: arma_hot
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::init_xform_mt(const SpBase<eT2,T1>& A, const Functor& func)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const SpProxy<T1> P(A.get_ref());
Chris@49:   
Chris@49:   if( (P.is_alias(*this) == true) || (is_SpMat<typename SpProxy<T1>::stored_type>::value == true) )
Chris@49:     {
Chris@49:     // NOTE: unwrap_spmat will convert a submatrix to a matrix, which in effect takes care of aliasing with submatrices;
Chris@49:     // NOTE: however, when more delayed ops are implemented, more elaborate handling of aliasing will be necessary
Chris@49:     const unwrap_spmat<typename SpProxy<T1>::stored_type> tmp(P.Q);
Chris@49:     
Chris@49:     const SpMat<eT2>& x = tmp.M;
Chris@49:     
Chris@49:     if(void_ptr(this) != void_ptr(&x))
Chris@49:       {
Chris@49:       init(x.n_rows, x.n_cols);
Chris@49:       
Chris@49:       // values and row_indices may not be null.
Chris@49:       if(values != NULL)
Chris@49:         {
Chris@49:         memory::release(values);
Chris@49:         memory::release(row_indices);
Chris@49:         }
Chris@49:       
Chris@49:       access::rw(values)      = memory::acquire_chunked<eT>   (x.n_nonzero + 1);
Chris@49:       access::rw(row_indices) = memory::acquire_chunked<uword>(x.n_nonzero + 1);
Chris@49:       
Chris@49:       arrayops::copy(access::rwp(row_indices), x.row_indices, x.n_nonzero + 1);
Chris@49:       arrayops::copy(access::rwp(col_ptrs),    x.col_ptrs,    x.n_cols    + 1);
Chris@49:       
Chris@49:       access::rw(n_nonzero) = x.n_nonzero;
Chris@49:       }
Chris@49:     
Chris@49:     
Chris@49:     // initialise the elements array with a transformed version of the elements from x
Chris@49:     
Chris@49:     const uword nnz = n_nonzero;
Chris@49:     
Chris@49:     const eT2* x_values = x.values;
Chris@49:           eT*  t_values = access::rwp(values);
Chris@49:     
Chris@49:     for(uword i=0; i < nnz; ++i)
Chris@49:       {
Chris@49:       t_values[i] = func(x_values[i]);   // NOTE: func() must produce a value of type eT (ie. act as a convertor between eT2 and eT)
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     init(P.get_n_rows(), P.get_n_cols());
Chris@49:     
Chris@49:     mem_resize(P.get_n_nonzero());
Chris@49:     
Chris@49:     typename SpProxy<T1>::const_iterator_type it = P.begin();
Chris@49:     
Chris@49:     while(it != P.end())
Chris@49:       {
Chris@49:       access::rw(row_indices[it.pos()]) = it.row();
Chris@49:       access::rw(values[it.pos()]) = func(*it);   // NOTE: func() must produce a value of type eT (ie. act as a convertor between eT2 and eT)
Chris@49:       ++access::rw(col_ptrs[it.col() + 1]);
Chris@49:       ++it;
Chris@49:       }
Chris@49:     
Chris@49:     // Now sum column pointers.
Chris@49:     for(uword c = 1; c <= n_cols; ++c)
Chris@49:       {
Chris@49:       access::rw(col_ptrs[c]) += col_ptrs[c - 1];
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::iterator
Chris@49: SpMat<eT>::begin()
Chris@49:   {
Chris@49:   return iterator(*this);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_iterator
Chris@49: SpMat<eT>::begin() const
Chris@49:   {
Chris@49:   return const_iterator(*this);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::iterator
Chris@49: SpMat<eT>::end()
Chris@49:   {
Chris@49:   return iterator(*this, 0, n_cols, n_nonzero);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_iterator
Chris@49: SpMat<eT>::end() const
Chris@49:   {
Chris@49:   return const_iterator(*this, 0, n_cols, n_nonzero);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::iterator
Chris@49: SpMat<eT>::begin_col(const uword col_num)
Chris@49:   {
Chris@49:   return iterator(*this, 0, col_num);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_iterator
Chris@49: SpMat<eT>::begin_col(const uword col_num) const
Chris@49:   {
Chris@49:   return const_iterator(*this, 0, col_num);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::iterator
Chris@49: SpMat<eT>::end_col(const uword col_num)
Chris@49:   {
Chris@49:   return iterator(*this, 0, col_num + 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_iterator
Chris@49: SpMat<eT>::end_col(const uword col_num) const
Chris@49:   {
Chris@49:   return const_iterator(*this, 0, col_num + 1);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::row_iterator
Chris@49: SpMat<eT>::begin_row(const uword row_num)
Chris@49:   {
Chris@49:   return row_iterator(*this, row_num, 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_row_iterator
Chris@49: SpMat<eT>::begin_row(const uword row_num) const
Chris@49:   {
Chris@49:   return const_row_iterator(*this, row_num, 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::row_iterator
Chris@49: SpMat<eT>::end_row()
Chris@49:   {
Chris@49:   return row_iterator(*this, n_nonzero);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_row_iterator
Chris@49: SpMat<eT>::end_row() const
Chris@49:   {
Chris@49:   return const_row_iterator(*this, n_nonzero);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::row_iterator
Chris@49: SpMat<eT>::end_row(const uword row_num)
Chris@49:   {
Chris@49:   return row_iterator(*this, row_num + 1, 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: typename SpMat<eT>::const_row_iterator
Chris@49: SpMat<eT>::end_row(const uword row_num) const
Chris@49:   {
Chris@49:   return const_row_iterator(*this, row_num + 1, 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: SpMat<eT>::clear()
Chris@49:   {
Chris@49:   if (values)
Chris@49:     {
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49:     
Chris@49:     access::rw(values)      = memory::acquire_chunked<eT>   (1);
Chris@49:     access::rw(row_indices) = memory::acquire_chunked<uword>(1);
Chris@49: 
Chris@49:     access::rw(values[0]) = 0;
Chris@49:     access::rw(row_indices[0]) = 0;
Chris@49:     }
Chris@49:   
Chris@49:   memory::release(col_ptrs);
Chris@49:   
Chris@49:   access::rw(col_ptrs) = memory::acquire<uword>(n_cols + 2);
Chris@49:   access::rw(col_ptrs[n_cols + 1]) = std::numeric_limits<uword>::max();
Chris@49:   
Chris@49:   arrayops::inplace_set(col_ptrs, eT(0), n_cols + 1);
Chris@49:   
Chris@49:   access::rw(n_nonzero) = 0;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: bool
Chris@49: SpMat<eT>::empty() const
Chris@49:   {
Chris@49:   return (n_elem == 0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: uword
Chris@49: SpMat<eT>::size() const
Chris@49:   {
Chris@49:   return n_elem;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::get_value(const uword i)
Chris@49:   {
Chris@49:   // First convert to the actual location.
Chris@49:   uword lcol = i / n_rows; // Integer division.
Chris@49:   uword lrow = i % n_rows;
Chris@49: 
Chris@49:   return get_value(lrow, lcol);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::get_value(const uword i) const
Chris@49:   {
Chris@49:   // First convert to the actual location.
Chris@49:   uword lcol = i / n_rows; // Integer division.
Chris@49:   uword lrow = i % n_rows;
Chris@49: 
Chris@49:   return get_value(lrow, lcol);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: SpValProxy<SpMat<eT> >
Chris@49: SpMat<eT>::get_value(const uword in_row, const uword in_col)
Chris@49:   {
Chris@49:   const uword colptr      = col_ptrs[in_col];
Chris@49:   const uword next_colptr = col_ptrs[in_col + 1];
Chris@49: 
Chris@49:   // Step through the row indices to see if our element exists.
Chris@49:   for (uword i = colptr; i < next_colptr; ++i)
Chris@49:     {
Chris@49:     const uword row_index = row_indices[i];
Chris@49:     
Chris@49:     // First check that we have not stepped past it.
Chris@49:     if (in_row < row_index) // If we have, then it doesn't exist: return 0.
Chris@49:       {
Chris@49:       return SpValProxy<SpMat<eT> >(in_row, in_col, *this); // Proxy for a zero value.
Chris@49:       }
Chris@49: 
Chris@49:     // Now check if we are at the correct place.
Chris@49:     if (in_row == row_index) // If we are, return a reference to the value.
Chris@49:       {
Chris@49:       return SpValProxy<SpMat<eT> >(in_row, in_col, *this, &access::rw(values[i]));
Chris@49:       }
Chris@49: 
Chris@49:     }
Chris@49: 
Chris@49:   // We did not find it, so it does not exist: return 0.
Chris@49:   return SpValProxy<SpMat<eT> >(in_row, in_col, *this);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: eT
Chris@49: SpMat<eT>::get_value(const uword in_row, const uword in_col) const
Chris@49:   {
Chris@49:   const uword colptr      = col_ptrs[in_col];
Chris@49:   const uword next_colptr = col_ptrs[in_col + 1];
Chris@49:   
Chris@49:   // Step through the row indices to see if our element exists.
Chris@49:   for (uword i = colptr; i < next_colptr; ++i)
Chris@49:     {
Chris@49:     const uword row_index = row_indices[i];
Chris@49:     
Chris@49:     // First check that we have not stepped past it.
Chris@49:     if (in_row < row_index) // If we have, then it doesn't exist: return 0.
Chris@49:       {
Chris@49:       return eT(0);
Chris@49:       }
Chris@49:     
Chris@49:     // Now check if we are at the correct place.
Chris@49:     if (in_row == row_index) // If we are, return the value.
Chris@49:       {
Chris@49:       return values[i];
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   // We did not find it, so it does not exist: return 0.
Chris@49:   return eT(0);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Given the index representing which of the nonzero values this is, return its
Chris@49:  * actual location, either in row/col or just the index.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: arma_hot
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: uword
Chris@49: SpMat<eT>::get_position(const uword i) const
Chris@49:   {
Chris@49:   uword lrow, lcol;
Chris@49:   
Chris@49:   get_position(i, lrow, lcol);
Chris@49:   
Chris@49:   // Assemble the row/col into the element's location in the matrix.
Chris@49:   return (lrow + n_rows * lcol);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_hot
Chris@49: arma_inline
Chris@49: void
Chris@49: SpMat<eT>::get_position(const uword i, uword& row_of_i, uword& col_of_i) const
Chris@49:   {
Chris@49:   arma_debug_check((i >= n_nonzero), "SpMat::get_position(): index out of bounds");
Chris@49:   
Chris@49:   col_of_i = 0;
Chris@49:   while (col_ptrs[col_of_i + 1] <= i)
Chris@49:     {
Chris@49:     col_of_i++;
Chris@49:     }
Chris@49:   
Chris@49:   row_of_i = row_indices[i];
Chris@49:   
Chris@49:   return;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Add an element at the given position, and return a reference to it.  The
Chris@49:  * element will be set to 0 (unless otherwise specified).  If the element
Chris@49:  * already exists, its value will be overwritten.
Chris@49:  *
Chris@49:  * @param in_row Row of new element.
Chris@49:  * @param in_col Column of new element.
Chris@49:  * @param in_val Value to set new element to (default 0.0).
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: eT&
Chris@49: SpMat<eT>::add_element(const uword in_row, const uword in_col, const eT val)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   // We will assume the new element does not exist and begin the search for
Chris@49:   // where to insert it.  If we find that it already exists, we will then
Chris@49:   // overwrite it.
Chris@49:   uword colptr      = col_ptrs[in_col    ];
Chris@49:   uword next_colptr = col_ptrs[in_col + 1];
Chris@49:   
Chris@49:   uword pos = colptr; // The position in the matrix of this value.
Chris@49:   
Chris@49:   if (colptr != next_colptr)
Chris@49:     {
Chris@49:     // There are other elements in this column, so we must find where this
Chris@49:     // element will fit as compared to those.
Chris@49:     while (pos < next_colptr && in_row > row_indices[pos])
Chris@49:       {
Chris@49:       pos++;
Chris@49:       }
Chris@49:     
Chris@49:     // We aren't inserting into the last position, so it is still possible
Chris@49:     // that the element may exist.
Chris@49:     if (pos != next_colptr && row_indices[pos] == in_row)
Chris@49:       {
Chris@49:       // It already exists.  Then, just overwrite it.
Chris@49:       access::rw(values[pos]) = val;
Chris@49:       
Chris@49:       return access::rw(values[pos]);
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   // 
Chris@49:   // Element doesn't exist, so we have to insert it
Chris@49:   // 
Chris@49:   
Chris@49:   // We have to update the rest of the column pointers.
Chris@49:   for (uword i = in_col + 1; i < n_cols + 1; i++)
Chris@49:     {
Chris@49:     access::rw(col_ptrs[i])++; // We are only inserting one new element.
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   // Figure out the actual amount of memory currently allocated
Chris@49:   // NOTE: this relies on memory::acquire_chunked() being used for the 'values' and 'row_indices' arrays
Chris@49:   const uword n_alloc = memory::enlarge_to_mult_of_chunksize(n_nonzero + 1);
Chris@49:   
Chris@49:   // If possible, avoid time-consuming memory allocation 
Chris@49:   if(n_alloc > (n_nonzero + 1))
Chris@49:     {
Chris@49:     arrayops::copy_backwards(access::rwp(values)      + pos + 1, values      + pos, (n_nonzero - pos) + 1);
Chris@49:     arrayops::copy_backwards(access::rwp(row_indices) + pos + 1, row_indices + pos, (n_nonzero - pos) + 1);
Chris@49:     
Chris@49:     // Insert the new element.
Chris@49:     access::rw(values[pos])      = val;
Chris@49:     access::rw(row_indices[pos]) = in_row;
Chris@49:     
Chris@49:     access::rw(n_nonzero)++;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword old_n_nonzero = n_nonzero;
Chris@49:     
Chris@49:     access::rw(n_nonzero)++; // Add to count of nonzero elements.
Chris@49:     
Chris@49:     // Allocate larger memory.
Chris@49:     eT*    new_values      = memory::acquire_chunked<eT>   (n_nonzero + 1);
Chris@49:     uword* new_row_indices = memory::acquire_chunked<uword>(n_nonzero + 1);
Chris@49:     
Chris@49:     // Copy things over, before the new element.
Chris@49:     if (pos > 0)
Chris@49:       {
Chris@49:       arrayops::copy(new_values,      values,      pos);
Chris@49:       arrayops::copy(new_row_indices, row_indices, pos);
Chris@49:       }
Chris@49:     
Chris@49:     // Insert the new element.
Chris@49:     new_values[pos]      = val;
Chris@49:     new_row_indices[pos] = in_row;
Chris@49:     
Chris@49:     // Copy the rest of things over (including the extra element at the end).
Chris@49:     arrayops::copy(new_values      + pos + 1, values      + pos, (old_n_nonzero - pos) + 1);
Chris@49:     arrayops::copy(new_row_indices + pos + 1, row_indices + pos, (old_n_nonzero - pos) + 1);
Chris@49:     
Chris@49:     // Assign new pointers.
Chris@49:     memory::release(values);
Chris@49:     memory::release(row_indices);
Chris@49: 
Chris@49:     access::rw(values)      = new_values;
Chris@49:     access::rw(row_indices) = new_row_indices;
Chris@49:     }
Chris@49:   
Chris@49:   return access::rw(values[pos]);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: /**
Chris@49:  * Delete an element at the given position.
Chris@49:  *
Chris@49:  * @param in_row Row of element to be deleted.
Chris@49:  * @param in_col Column of element to be deleted.
Chris@49:  */
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: arma_hot
Chris@49: void
Chris@49: SpMat<eT>::delete_element(const uword in_row, const uword in_col)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   // We assume the element exists (although... it may not) and look for its
Chris@49:   // exact position.  If it doesn't exist... well, we don't need to do anything.
Chris@49:   uword colptr      = col_ptrs[in_col];
Chris@49:   uword next_colptr = col_ptrs[in_col + 1];
Chris@49: 
Chris@49:   if (colptr != next_colptr)
Chris@49:     {
Chris@49:     // There's at least one element in this column.
Chris@49:     // Let's see if we are one of them.
Chris@49:     for (uword pos = colptr; pos < next_colptr; pos++)
Chris@49:       {
Chris@49:       if (in_row == row_indices[pos])
Chris@49:         {
Chris@49:         const uword old_n_nonzero = n_nonzero;
Chris@49:         
Chris@49:         --access::rw(n_nonzero); // Remove one from the count of nonzero elements.
Chris@49:         
Chris@49:         // Found it.  Now remove it.
Chris@49:           
Chris@49:         // Figure out the actual amount of memory currently allocated and the actual amount that will be required
Chris@49:         // NOTE: this relies on memory::acquire_chunked() being used for the 'values' and 'row_indices' arrays
Chris@49:           
Chris@49:         const uword n_alloc     = memory::enlarge_to_mult_of_chunksize(old_n_nonzero + 1);
Chris@49:         const uword n_alloc_mod = memory::enlarge_to_mult_of_chunksize(n_nonzero + 1);
Chris@49:           
Chris@49:         // If possible, avoid time-consuming memory allocation
Chris@49:         if(n_alloc_mod == n_alloc)
Chris@49:           {
Chris@49:           if (pos < n_nonzero)  // remember, we decremented n_nonzero
Chris@49:             {
Chris@49:             arrayops::copy_forwards(access::rwp(values)      + pos, values + pos + 1, (n_nonzero - pos) + 1);
Chris@49:             arrayops::copy_forwards(access::rwp(row_indices) + pos, row_indices + pos + 1, (n_nonzero - pos) + 1);
Chris@49:             }
Chris@49:           }
Chris@49:         else
Chris@49:           {
Chris@49:           // Make new arrays.
Chris@49:           eT*    new_values      = memory::acquire_chunked<eT>   (n_nonzero + 1);
Chris@49:           uword* new_row_indices = memory::acquire_chunked<uword>(n_nonzero + 1);
Chris@49:             
Chris@49:           if (pos > 0)
Chris@49:             {
Chris@49:             arrayops::copy(new_values,      values,      pos);
Chris@49:             arrayops::copy(new_row_indices, row_indices, pos);
Chris@49:             }
Chris@49:             
Chris@49:           arrayops::copy(new_values      + pos, values      + pos + 1, (n_nonzero - pos) + 1);
Chris@49:           arrayops::copy(new_row_indices + pos, row_indices + pos + 1, (n_nonzero - pos) + 1);
Chris@49:             
Chris@49:           memory::release(values);
Chris@49:           memory::release(row_indices);
Chris@49:             
Chris@49:           access::rw(values)      = new_values;
Chris@49:           access::rw(row_indices) = new_row_indices;
Chris@49:           }
Chris@49:         
Chris@49:         // And lastly, update all the column pointers (decrement by one).
Chris@49:         for (uword i = in_col + 1; i < n_cols + 1; i++)
Chris@49:           {
Chris@49:           --access::rw(col_ptrs[i]); // We only removed one element.
Chris@49:           }
Chris@49:         
Chris@49:         return; // There is nothing left to do.
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return; // The element does not exist, so there's nothing for us to do.
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: #ifdef ARMA_EXTRA_SPMAT_MEAT
Chris@49:   #include ARMA_INCFILE_WRAP(ARMA_EXTRA_SPMAT_MEAT)
Chris@49: #endif
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}