| SpMat_meat.hpp | | SpMat_meat.hpp | |
|---|
| // Copyright (C) 2011-2013 Ryan Curtin | | // Copyright (C) 2011-2013 Ryan Curtin | |
|
| // Copyright (C) 2012-2013 Conrad Sanderson | | // Copyright (C) 2012-2014 Conrad Sanderson | |
| // Copyright (C) 2011 Matthew Amidon | | // Copyright (C) 2011 Matthew Amidon | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup SpMat | | //! \addtogroup SpMat | |
| //! @{ | | //! @{ | |
| | | | |
| /** | | /** | |
| | | | |
| skipping to change at line 198 | | skipping to change at line 198 | |
|---|
| , n_cols(0) | | , n_cols(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , n_nonzero(0) | | , n_nonzero(0) | |
| , vec_state(0) | | , vec_state(0) | |
| , values(NULL) | | , values(NULL) | |
| , row_indices(NULL) | | , row_indices(NULL) | |
| , col_ptrs(NULL) | | , col_ptrs(NULL) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| | | | |
|
| const unwrap<T1> locs_tmp( locations_expr.get_ref() ); | | const unwrap<T1> locs_tmp( locations_expr.get_ref() ); | |
| const Mat<uword>& locs = locs_tmp.M; | | const unwrap<T2> vals_tmp( vals_expr.get_ref() ); | |
| | | | |
| const unwrap<T2> vals_tmp( vals_expr.get_ref() ); | | | |
| const Mat<eT>& vals = vals_tmp.M; | | | |
| | | | |
|
| arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'value | | const Mat<uword>& locs = locs_tmp.M; | |
| s' object is not a vector" ); | | const Mat<eT>& vals = vals_tmp.M; | |
| | | | |
|
| arma_debug_check((locs.n_cols != vals.n_elem), "SpMat::SpMat(): number of | | arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'v | |
| locations is different than number of values"); | | alues' object is not a vector" ); | |
| | | arma_debug_check( (locs.n_rows != 2), "SpMat::SpMat(): location | |
| | | s matrix must have two rows" ); | |
| | | arma_debug_check( (locs.n_cols != vals.n_elem), "SpMat::SpMat(): number o | |
| | | f locations is different than number of values" ); | |
| | | | |
| // If there are no elements in the list, max() will fail. | | // If there are no elements in the list, max() will fail. | |
|
| if (locs.n_cols == 0) | | if(locs.n_cols == 0) { init(0, 0); return; } | |
| { | | | |
| init(0, 0); | | | |
| return; | | | |
| } | | | |
| | | | |
| arma_debug_check((locs.n_rows != 2), "SpMat::SpMat(): locations matrix mu | | | |
| st have two rows"); | | | |
| | | | |
|
| // Automatically determine size (and check if it's sorted). | | // Automatically determine size | |
| uvec bounds = arma::max(locs, 1); | | uvec bounds = arma::max(locs, 1); | |
| init(bounds[0] + 1, bounds[1] + 1); | | init(bounds[0] + 1, bounds[1] + 1); | |
| | | | |
|
| // Resize to correct number of elements. | | init_batch(locs, vals, sort_locations); | |
| mem_resize(vals.n_elem); | | | |
| | | | |
| // Reset column pointers to zero. | | | |
| arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1); | | | |
| | | | |
| bool actually_sorted = true; | | | |
| if(sort_locations == true) | | | |
| { | | | |
| // sort_index() uses std::sort() which may use quicksort... so we bette | | | |
| r | | | |
| // make sure it's not already sorted before taking an O(N^2) sort penal | | | |
| ty. | | | |
| for (uword i = 1; i < locs.n_cols; ++i) | | | |
| { | | | |
| 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))) | | | |
| { | | | |
| actually_sorted = false; | | | |
| break; | | | |
| } | | | |
| } | | | |
| | | | |
| if(actually_sorted == false) | | | |
| { | | | |
| // This may not be the fastest possible implementation but it maximiz | | | |
| es code reuse. | | | |
| Col<uword> abslocs(locs.n_cols); | | | |
| | | | |
| for (uword i = 0; i < locs.n_cols; ++i) | | | |
| { | | | |
| abslocs[i] = locs.at(1, i) * n_rows + locs.at(0, i); | | | |
| } | | | |
| | | | |
| // Now we will sort with sort_index(). | | | |
| uvec sorted_indices = sort_index(abslocs); // Ascending sort. | | | |
| | | | |
| // Now we add the elements in this sorted order. | | | |
| for (uword i = 0; i < sorted_indices.n_elem; ++i) | | | |
| { | | | |
| arma_debug_check((locs.at(0, sorted_indices[i]) >= n_rows), "SpMat: | | | |
| :SpMat(): invalid row index"); | | | |
| arma_debug_check((locs.at(1, sorted_indices[i]) >= n_cols), "SpMat: | | | |
| :SpMat(): invalid column index"); | | | |
| | | | |
| access::rw(values[i]) = vals[sorted_indices[i]]; | | | |
| access::rw(row_indices[i]) = locs.at(0, sorted_indices[i]); | | | |
| | | | |
| access::rw(col_ptrs[locs.at(1, sorted_indices[i]) + 1])++; | | | |
| } | | | |
| } | | | |
| } | | | |
| | | | |
| if( (sort_locations == false) || (actually_sorted == true) ) | | | |
| { | | | |
| // Now set the values and row indices correctly. | | | |
| // Increment the column pointers in each column (so they are column "co | | | |
| unts"). | | | |
| for (uword i = 0; i < vals.n_elem; ++i) | | | |
| { | | | |
| arma_debug_check((locs.at(0, i) >= n_rows), "SpMat::SpMat(): invalid | | | |
| row index"); | | | |
| arma_debug_check((locs.at(1, i) >= n_cols), "SpMat::SpMat(): invalid | | | |
| column index"); | | | |
| | | | |
| // Check ordering in debug mode. | | | |
| if(i > 0) | | | |
| { | | | |
| arma_debug_check | | | |
| ( | | | |
| ( (locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.a | | | |
| t(1, i - 1) && locs.at(0, i) < locs.at(0, i - 1)) ), | | | |
| "SpMat::SpMat(): out of order points; either pass sort_locations | | | |
| = true, or sort points in column-major ordering" | | | |
| ); | | | |
| 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"); | | | |
| } | | | |
| | | | |
| access::rw(values[i]) = vals[i]; | | | |
| access::rw(row_indices[i]) = locs.at(0, i); | | | |
| | | | |
| access::rw(col_ptrs[locs.at(1, i) + 1])++; | | | |
| } | | | |
| } | | | |
| | | | |
| // Now fix the column pointers. | | | |
| for (uword i = 0; i <= n_cols; ++i) | | | |
| { | | | |
| access::rw(col_ptrs[i + 1]) += col_ptrs[i]; | | | |
| } | | | |
| } | | } | |
| | | | |
| //! Insert a large number of values at once. | | //! Insert a large number of values at once. | |
| //! locations.row[0] should be row indices, locations.row[1] should be colu
mn indices, | | //! locations.row[0] should be row indices, locations.row[1] should be colu
mn indices, | |
| //! and values should be the corresponding values. | | //! and values should be the corresponding values. | |
| //! If sort_locations is false, then it is assumed that the locations and v
alues | | //! If sort_locations is false, then it is assumed that the locations and v
alues | |
| //! are already sorted in column-major ordering. | | //! are already sorted in column-major ordering. | |
| //! In this constructor the size is explicitly given. | | //! In this constructor the size is explicitly given. | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1, typename T2> | | template<typename T1, typename T2> | |
| | | | |
| skipping to change at line 323 | | skipping to change at line 239 | |
|---|
| , n_cols(0) | | , n_cols(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , n_nonzero(0) | | , n_nonzero(0) | |
| , vec_state(0) | | , vec_state(0) | |
| , values(NULL) | | , values(NULL) | |
| , row_indices(NULL) | | , row_indices(NULL) | |
| , col_ptrs(NULL) | | , col_ptrs(NULL) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| | | | |
|
| init(in_n_rows, in_n_cols); | | | |
| | | | |
| const unwrap<T1> locs_tmp( locations_expr.get_ref() ); | | const unwrap<T1> locs_tmp( locations_expr.get_ref() ); | |
|
| const unwrap<T2> vals_tmp( vals_expr.get_ref() ); | | const unwrap<T2> vals_tmp( vals_expr.get_ref() ); | |
| | | | |
| const Mat<uword>& locs = locs_tmp.M; | | const Mat<uword>& locs = locs_tmp.M; | |
| const Mat<eT>& vals = vals_tmp.M; | | const Mat<eT>& vals = vals_tmp.M; | |
| | | | |
|
| arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'v | | arma_debug_check( (vals.is_vec() == false), "SpMat::SpMat(): given 'v | |
| alues' object is not a vector" ); | | alues' object is not a vector" ); | |
| arma_debug_check( (locs.n_rows != 2), "SpMat::SpMat(): location | | arma_debug_check( (locs.n_rows != 2), "SpMat::SpMat(): location | |
| s matrix must have two rows" ); | | s matrix must have two rows" ); | |
| arma_debug_check( (locs.n_cols != vals.n_elem), "SpMat::SpMat(): number o
f locations is different than number of values" ); | | arma_debug_check( (locs.n_cols != vals.n_elem), "SpMat::SpMat(): number o
f locations is different than number of values" ); | |
| | | | |
|
| // Resize to correct number of elements. | | init(in_n_rows, in_n_cols); | |
| mem_resize(vals.n_elem); | | | |
| | | | |
| // Reset column pointers to zero. | | | |
| arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1); | | | |
| | | | |
| bool actually_sorted = true; | | | |
| | | | |
| if(sort_locations == true) | | | |
| { | | | |
| // sort_index() uses std::sort() which may use quicksort... so we bette | | | |
| r | | | |
| // make sure it's not already sorted before taking an O(N^2) sort penal | | | |
| ty. | | | |
| for (uword i = 1; i < locs.n_cols; ++i) | | | |
| { | | | |
| 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))) | | | |
| { | | | |
| actually_sorted = false; | | | |
| break; | | | |
| } | | | |
| } | | | |
| | | | |
| if(actually_sorted == false) | | | |
| { | | | |
| // This may not be the fastest possible implementation but it maximiz | | | |
| es code reuse. | | | |
| Col<uword> abslocs(locs.n_cols); | | | |
| | | | |
| for (uword i = 0; i < locs.n_cols; ++i) | | | |
| { | | | |
| abslocs[i] = locs.at(1, i) * n_rows + locs.at(0, i); | | | |
| } | | | |
| | | | |
| // Now we will sort with sort_index(). | | | |
| uvec sorted_indices = sort_index(abslocs); // Ascending sort. | | | |
| | | | |
| // Now we add the elements in this sorted order. | | | |
| for (uword i = 0; i < sorted_indices.n_elem; ++i) | | | |
| { | | | |
| arma_debug_check((locs.at(0, sorted_indices[i]) >= n_rows), "SpMat: | | | |
| :SpMat(): invalid row index"); | | | |
| arma_debug_check((locs.at(1, sorted_indices[i]) >= n_cols), "SpMat: | | | |
| :SpMat(): invalid column index"); | | | |
| | | | |
| access::rw(values[i]) = vals[sorted_indices[i]]; | | | |
| access::rw(row_indices[i]) = locs.at(0, sorted_indices[i]); | | | |
| | | | |
| access::rw(col_ptrs[locs.at(1, sorted_indices[i]) + 1])++; | | | |
| } | | | |
| } | | | |
| } | | | |
| | | | |
| if( (sort_locations == false) || (actually_sorted == true) ) | | | |
| { | | | |
| // Now set the values and row indices correctly. | | | |
| // Increment the column pointers in each column (so they are column "co | | | |
| unts"). | | | |
| for (uword i = 0; i < vals.n_elem; ++i) | | | |
| { | | | |
| arma_debug_check((locs.at(0, i) >= n_rows), "SpMat::SpMat(): invalid | | | |
| row index"); | | | |
| arma_debug_check((locs.at(1, i) >= n_cols), "SpMat::SpMat(): invalid | | | |
| column index"); | | | |
| | | | |
| // Check ordering in debug mode. | | | |
| if(i > 0) | | | |
| { | | | |
| arma_debug_check | | | |
| ( | | | |
| ( (locs.at(1, i) < locs.at(1, i - 1)) || (locs.at(1, i) == locs.a | | | |
| t(1, i - 1) && locs.at(0, i) < locs.at(0, i - 1)) ), | | | |
| "SpMat::SpMat(): out of order points; either pass sort_locations | | | |
| = true or sort points in column-major ordering" | | | |
| ); | | | |
| 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"); | | | |
| } | | | |
| //! If sort_locations is false, then it is assumed that the locations and v | | | |
| alues | | | |
| //! are already sorted in column-major ordering. | | | |
| | | | |
| access::rw(values[i]) = vals[i]; | | | |
| access::rw(row_indices[i]) = locs.at(0, i); | | | |
| | | | |
| access::rw(col_ptrs[locs.at(1, i) + 1])++; | | | |
| } | | | |
| } | | | |
| | | | |
|
| // Now fix the column pointers. | | init_batch(locs, vals, sort_locations); | |
| for (uword i = 0; i <= n_cols; ++i) | | | |
| { | | | |
| access::rw(col_ptrs[i + 1]) += col_ptrs[i]; | | | |
| } | | | |
| } | | } | |
| | | | |
| //! Insert a large number of values at once. | | //! Insert a large number of values at once. | |
| //! Per CSC format, rowind_expr should be row indices, | | //! Per CSC format, rowind_expr should be row indices, | |
| //! colptr_expr should column ptr indices locations, | | //! colptr_expr should column ptr indices locations, | |
| //! and values should be the corresponding values. | | //! and values should be the corresponding values. | |
| //! In this constructor the size is explicitly given. | | //! In this constructor the size is explicitly given. | |
| //! Values are assumed to be sorted, and the size | | //! Values are assumed to be sorted, and the size | |
| //! information is trusted | | //! information is trusted | |
| template<typename eT> | | template<typename eT> | |
| | | | |
| skipping to change at line 477 | | skipping to change at line 312 | |
|---|
| | | | |
| // copy supplied values into sparse matrix -- not checked for consistency | | // copy supplied values into sparse matrix -- not checked for consistency | |
| arrayops::copy(access::rwp(row_indices), rowind.memptr(), rowind.n_elem )
; | | arrayops::copy(access::rwp(row_indices), rowind.memptr(), rowind.n_elem )
; | |
| arrayops::copy(access::rwp(col_ptrs), colptr.memptr(), colptr.n_elem )
; | | arrayops::copy(access::rwp(col_ptrs), colptr.memptr(), colptr.n_elem )
; | |
| arrayops::copy(access::rwp(values), vals.memptr(), vals.n_elem )
; | | arrayops::copy(access::rwp(values), vals.memptr(), vals.n_elem )
; | |
| | | | |
| // important: set the sentinel as well | | // important: set the sentinel as well | |
| access::rw(col_ptrs[n_cols + 1]) = std::numeric_limits<uword>::max(); | | access::rw(col_ptrs[n_cols + 1]) = std::numeric_limits<uword>::max(); | |
| } | | } | |
| | | | |
|
| /** | | | |
| * Simple operators with plain values. These operate on every value in the | | | |
| * matrix, so a sparse matrix += 1 will turn all those zeroes into ones. B | | | |
| e | | | |
| * careful and make sure that's what you really want! | | | |
| */ | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator=(const eT val) | | SpMat<eT>::operator=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| // Resize to 1x1 then set that to the right value. | | if(val != eT(0)) | |
| init(1, 1); // Sets col_ptrs to 0. | | { | |
| mem_resize(1); // One element. | | // Resize to 1x1 then set that to the right value. | |
| | | init(1, 1); // Sets col_ptrs to 0. | |
| | | mem_resize(1); // One element. | |
| | | | |
|
| // Manually set element. | | // Manually set element. | |
| access::rw(values[0]) = val; | | access::rw(values[0]) = val; | |
| access::rw(row_indices[0]) = 0; | | access::rw(row_indices[0]) = 0; | |
| access::rw(col_ptrs[1]) = 1; | | access::rw(col_ptrs[1]) = 1; | |
| | | } | |
| | | else | |
| | | { | |
| | | init(0, 0); | |
| | | } | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator*=(const eT val) | | SpMat<eT>::operator*=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| if(val == eT(0)) | | if(val != eT(0)) | |
| | | { | |
| | | arrayops::inplace_mul( access::rwp(values), val, n_nonzero ); | |
| | | } | |
| | | else | |
| { | | { | |
| // Everything will be zero. | | // Everything will be zero. | |
| init(n_rows, n_cols); | | init(n_rows, n_cols); | |
|
| return *this; | | | |
| } | | } | |
| | | | |
|
| arrayops::inplace_mul( access::rwp(values), val, n_nonzero ); | | | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator/=(const eT val) | | SpMat<eT>::operator/=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 553 | | skipping to change at line 391 | |
|---|
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator+=(const SpMat<eT>& x) | | SpMat<eT>::operator+=(const SpMat<eT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| SpMat<eT> out; | | SpMat<eT> out = (*this) + x; | |
| out = (*this) + x; | | | |
| steal_mem(out); | | steal_mem(out); | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator-=(const SpMat<eT>& x) | | SpMat<eT>::operator-=(const SpMat<eT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| SpMat<eT> out; | | SpMat<eT> out = (*this) - x; | |
| out = (*this) - x; | | | |
| steal_mem(out); | | steal_mem(out); | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator*=(const SpMat<eT>& y) | | SpMat<eT>::operator*=(const SpMat<eT>& y) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| SpMat<eT> z; | | SpMat<eT> z = (*this) * y; | |
| z = (*this) * y; | | | |
| steal_mem(z); | | steal_mem(z); | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| // This is in-place element-wise matrix multiplication. | | // This is in-place element-wise matrix multiplication. | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
|
| SpMat<eT>::operator%=(const SpMat<eT>& x) | | SpMat<eT>::operator%=(const SpMat<eT>& y) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element- | | SpMat<eT> z = (*this) % y; | |
| wise multiplication"); | | | |
| | | | |
| // We can do this with two iterators rather simply. | | | |
| iterator it = begin(); | | | |
| const_iterator x_it = x.begin(); | | | |
| | | | |
| while (it != end() && x_it != x.end()) | | | |
| { | | | |
| // One of these will be further advanced than the other (or they will b | | | |
| e at the same place). | | | |
| if ((it.row() == x_it.row()) && (it.col() == x_it.col())) | | | |
| { | | | |
| // There is an element at this place in both matrices. Multiply. | | | |
| (*it) *= (*x_it); | | | |
| | | | |
| // Now move on to the next position. | | | |
| it++; | | | |
| x_it++; | | | |
| } | | | |
| | | | |
| else if ((it.col() < x_it.col()) || ((it.col() == x_it.col()) && (it.ro | | | |
| w() < x_it.row()))) | | | |
| { | | | |
| // This case is when our matrix has an element which the other matrix | | | |
| does not. | | | |
| // So we must delete this element. | | | |
| (*it) = 0; | | | |
| | | | |
| // Because we have deleted the element, we now have to manually set t | | | |
| he position... | | | |
| it.internal_pos--; | | | |
| | | | |
| // Now we can increment our iterator. | | | |
| it++; | | | |
| } | | | |
| | | | |
| else /* if our iterator is ahead of the other matrix */ | | | |
| { | | | |
| // In this case we don't need to set anything to 0; our element is al | | | |
| ready 0. | | | |
| // We can just increment the iterator of the other matrix. | | | |
| x_it++; | | | |
| } | | | |
| | | | |
| } | | | |
| | | | |
| // If we are not at the end of our matrix, then we must terminate the rem | | | |
| aining elements. | | | |
| while (it != end()) | | | |
| { | | | |
| (*it) = 0; | | | |
| | | | |
|
| // Hack to manually set the position right... | | steal_mem(z); | |
| it.internal_pos--; | | | |
| it++; // ...and then an increment. | | | |
| } | | | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| // Construct a complex matrix out of two non-complex matrices | | // Construct a complex matrix out of two non-complex matrices | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1, typename T2> | | template<typename T1, typename T2> | |
| inline | | inline | |
| SpMat<eT>::SpMat | | SpMat<eT>::SpMat | |
| ( | | ( | |
| | | | |
| skipping to change at line 1089 | | skipping to change at line 880 | |
|---|
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator+=(const SpSubview<eT>& X) | | SpMat<eT>::operator+=(const SpSubview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition | | SpMat<eT> tmp = (*this) + X; | |
| "); | | | |
| | | | |
| typename SpSubview<eT>::const_iterator it = X.begin(); | | | |
| | | | |
|
| while(it != X.end()) | | steal_mem(tmp); | |
| { | | | |
| at(it.row(), it.col()) += (*it); | | | |
| ++it; | | | |
| } | | | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator-=(const SpSubview<eT>& X) | | SpMat<eT>::operator-=(const SpSubview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtract | | SpMat<eT> tmp = (*this) - X; | |
| ion"); | | | |
| | | | |
| typename SpSubview<eT>::const_iterator it = X.begin(); | | | |
| | | | |
|
| while(it != X.end()) | | steal_mem(tmp); | |
| { | | | |
| at(it.row(), it.col()) -= (*it); | | | |
| ++it; | | | |
| } | | | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator*=(const SpSubview<eT>& y) | | SpMat<eT>::operator*=(const SpSubview<eT>& y) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| arma_debug_assert_mul_size(n_rows, n_cols, y.n_rows, y.n_cols, "matrix mu | | | |
| ltiplication"); | | | |
| | | | |
| // Cannot be done in-place (easily). | | | |
| SpMat<eT> z = (*this) * y; | | SpMat<eT> z = (*this) * y; | |
|
| | | | |
| steal_mem(z); | | steal_mem(z); | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator%=(const SpSubview<eT>& x) | | SpMat<eT>::operator%=(const SpSubview<eT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element- | | SpMat<eT> tmp = (*this) % x; | |
| wise multiplication"); | | | |
| | | | |
| iterator it = begin(); | | | |
| typename SpSubview<eT>::const_iterator xit = x.begin(); | | | |
| | | | |
|
| while((it != end()) || (xit != x.end())) | | steal_mem(tmp); | |
| { | | | |
| if((xit.row() == it.row()) && (xit.col() == it.col())) | | | |
| { | | | |
| (*it) *= (*xit); | | | |
| ++it; | | | |
| ++xit; | | | |
| } | | | |
| else | | | |
| { | | | |
| if((xit.col() > it.col()) || ((xit.col() == it.col()) && (xit.row() > | | | |
| it.row()))) | | | |
| { | | | |
| // xit is "ahead" | | | |
| (*it) = eT(0); // erase element; x has a zero here | | | |
| it.internal_pos--; // update iterator so it still works | | | |
| ++it; | | | |
| } | | | |
| else | | | |
| { | | | |
| // it is "ahead" | | | |
| ++xit; | | | |
| } | | | |
| } | | | |
| } | | | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpMat<eT>& | | const SpMat<eT>& | |
| SpMat<eT>::operator/=(const SpSubview<eT>& x) | | SpMat<eT>::operator/=(const SpSubview<eT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| skipping to change at line 3950 | | skipping to change at line 3701 | |
|---|
| arrayops::copy(access::rwp(row_indices), x.row_indices, x.n_nonzero + 1
); | | arrayops::copy(access::rwp(row_indices), x.row_indices, x.n_nonzero + 1
); | |
| arrayops::copy(access::rwp(col_ptrs), x.col_ptrs, x.n_cols + 1); | | arrayops::copy(access::rwp(col_ptrs), x.col_ptrs, x.n_cols + 1); | |
| | | | |
| access::rw(n_nonzero) = x.n_nonzero; | | access::rw(n_nonzero) = x.n_nonzero; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
|
| | | SpMat<eT>::init_batch(const Mat<uword>& locs, const Mat<eT>& vals, const bo | |
| | | ol sort_locations) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | // Resize to correct number of elements. | |
| | | mem_resize(vals.n_elem); | |
| | | | |
| | | // Reset column pointers to zero. | |
| | | arrayops::inplace_set(access::rwp(col_ptrs), uword(0), n_cols + 1); | |
| | | | |
| | | bool actually_sorted = true; | |
| | | | |
| | | if(sort_locations == true) | |
| | | { | |
| | | // sort_index() uses std::sort() which may use quicksort... so we bette | |
| | | r | |
| | | // make sure it's not already sorted before taking an O(N^2) sort penal | |
| | | ty. | |
| | | for (uword i = 1; i < locs.n_cols; ++i) | |
| | | { | |
| | | const uword* locs_i = locs.colptr(i ); | |
| | | const uword* locs_im1 = locs.colptr(i-1); | |
| | | | |
| | | if( (locs_i[1] < locs_im1[1]) || (locs_i[1] == locs_im1[1] && locs_ | |
| | | i[0] <= locs_im1[0]) ) | |
| | | { | |
| | | actually_sorted = false; | |
| | | break; | |
| | | } | |
| | | } | |
| | | | |
| | | if(actually_sorted == false) | |
| | | { | |
| | | // This may not be the fastest possible implementation but it maximiz | |
| | | es code reuse. | |
| | | Col<uword> abslocs(locs.n_cols); | |
| | | | |
| | | for (uword i = 0; i < locs.n_cols; ++i) | |
| | | { | |
| | | const uword* locs_i = locs.colptr(i); | |
| | | | |
| | | abslocs[i] = locs_i[1] * n_rows + locs_i[0]; | |
| | | } | |
| | | | |
| | | uvec sorted_indices = sort_index(abslocs); // Ascending sort. | |
| | | | |
| | | // Now we add the elements in this sorted order. | |
| | | for (uword i = 0; i < sorted_indices.n_elem; ++i) | |
| | | { | |
| | | const uword* locs_i = locs.colptr( sorted_indices[i] ); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (locs_i[0] >= n_rows) || (locs_i[1] >= n_cols) ), | |
| | | "SpMat::SpMat(): invalid row or column index" | |
| | | ); | |
| | | | |
| | | if(i > 0) | |
| | | { | |
| | | const uword* locs_im1 = locs.colptr( sorted_indices[i-1] ); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (locs_i[1] == locs_im1[1]) && (locs_i[0] == locs_im1[0]) ), | |
| | | "SpMat::SpMat(): two identical point locations in list" | |
| | | ); | |
| | | } | |
| | | | |
| | | access::rw(values[i]) = vals[ sorted_indices[i] ]; | |
| | | access::rw(row_indices[i]) = locs_i[0]; | |
| | | | |
| | | access::rw(col_ptrs[ locs_i[1] + 1 ])++; | |
| | | } | |
| | | } | |
| | | } | |
| | | | |
| | | if( (sort_locations == false) || (actually_sorted == true) ) | |
| | | { | |
| | | // Now set the values and row indices correctly. | |
| | | // Increment the column pointers in each column (so they are column "co | |
| | | unts"). | |
| | | for (uword i = 0; i < vals.n_elem; ++i) | |
| | | { | |
| | | const uword* locs_i = locs.colptr(i); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (locs_i[0] >= n_rows) || (locs_i[1] >= n_cols) ), | |
| | | "SpMat::SpMat(): invalid row or column index" | |
| | | ); | |
| | | | |
| | | if(i > 0) | |
| | | { | |
| | | const uword* locs_im1 = locs.colptr(i-1); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (locs_i[1] < locs_im1[1]) || (locs_i[1] == locs_im1[1] && loc | |
| | | s_i[0] < locs_im1[0]) ), | |
| | | "SpMat::SpMat(): out of order points; either pass sort_locations | |
| | | = true, or sort points in column-major ordering" | |
| | | ); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (locs_i[1] == locs_im1[1]) && (locs_i[0] == locs_im1[0]) ), | |
| | | "SpMat::SpMat(): two identical point locations in list" | |
| | | ); | |
| | | } | |
| | | | |
| | | access::rw(values[i]) = vals[i]; | |
| | | access::rw(row_indices[i]) = locs_i[0]; | |
| | | | |
| | | access::rw(col_ptrs[ locs_i[1] + 1 ])++; | |
| | | } | |
| | | } | |
| | | | |
| | | // Now fix the column pointers. | |
| | | for (uword i = 0; i <= n_cols; ++i) | |
| | | { | |
| | | access::rw(col_ptrs[i + 1]) += col_ptrs[i]; | |
| | | } | |
| | | | |
| | | remove_zeros(); | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | inline | |
| | | void | |
| | | SpMat<eT>::remove_zeros() | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | uword zeros_count = 0; | |
| | | | |
| | | for(uword i=0; i<n_nonzero; ++i) | |
| | | { | |
| | | if(values[i] == eT(0)) { zeros_count++; } | |
| | | } | |
| | | | |
| | | if(zeros_count == 0) | |
| | | { | |
| | | return; | |
| | | } | |
| | | | |
| | | const uword actual_n_nonzero = n_nonzero - zeros_count; | |
| | | | |
| | | SpMat<eT> out(n_rows, n_cols); | |
| | | | |
| | | out.mem_resize(actual_n_nonzero); | |
| | | | |
| | | const SpMat<eT>& x = (*this); | |
| | | | |
| | | typename SpMat<eT>::const_iterator x_it = x.begin(); | |
| | | typename SpMat<eT>::const_iterator x_end = x.end(); | |
| | | | |
| | | uword cur_val = 0; | |
| | | while(x_it != x_end) | |
| | | { | |
| | | const eT val = (*x_it); | |
| | | | |
| | | if(val != eT(0)) | |
| | | { | |
| | | access::rw(out.values[cur_val]) = val; | |
| | | access::rw(out.row_indices[cur_val]) = x_it.row(); | |
| | | | |
| | | ++access::rw(out.col_ptrs[x_it.col() + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
| | | ++x_it; | |
| | | } | |
| | | | |
| | | const uword out_n_cols = out.n_cols; | |
| | | uword* out_col_ptrs = access::rwp(out.col_ptrs); | |
| | | | |
| | | for(uword c = 1; c <= out_n_cols; ++c) | |
| | | { | |
| | | out_col_ptrs[c] += out_col_ptrs[c - 1]; | |
| | | } | |
| | | | |
| | | steal_mem(out); | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | inline | |
| | | void | |
| SpMat<eT>::mem_resize(const uword new_n_nonzero) | | SpMat<eT>::mem_resize(const uword new_n_nonzero) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(n_nonzero != new_n_nonzero) | | if(n_nonzero != new_n_nonzero) | |
| { | | { | |
| if(new_n_nonzero == 0) | | if(new_n_nonzero == 0) | |
| { | | { | |
| memory::release(values); | | memory::release(values); | |
| memory::release(row_indices); | | memory::release(row_indices); | |
| | | | |
| End of changes. 33 change blocks. |
|---|
| 345 lines changed or deleted | | 244 lines changed or added | |
|---|
|
| SpSubview_meat.hpp | | SpSubview_meat.hpp | |
|---|
| // Copyright (C) 2011-2012 Ryan Curtin | | // Copyright (C) 2011-2012 Ryan Curtin | |
|
| | | // Copyright (C) 2012-2014 Conrad Sanderson | |
| // Copyright (C) 2011 Matthew Amidon | | // Copyright (C) 2011 Matthew Amidon | |
|
| // Copyright (C) 2012 Conrad Sanderson | | | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup SpSubview | | //! \addtogroup SpSubview | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_inline | | arma_inline | |
| | | | |
| skipping to change at line 89 | | skipping to change at line 89 | |
|---|
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator+=(const eT val) | | SpSubview<eT>::operator+=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(val == eT(0)) | | if(val == eT(0)) | |
| { | | { | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
|
| const uword lstart_row = aux_row1; | | Mat<eT> tmp( (*this).n_rows, (*this).n_cols ); | |
| const uword lend_row = aux_row1 + n_rows; | | | |
| | | | |
| const uword lstart_col = aux_col1; | | | |
| const uword lend_col = aux_col1 + n_cols; | | | |
| | | | |
| const uword old_n_nonzero = m.n_nonzero; | | | |
| | | | |
|
| // iterate over our part of the sparse matrix | | tmp.fill(val); | |
| for(uword lcol = lstart_col; lcol < lend_col; ++lcol) | | | |
| for(uword lrow = lstart_row; lrow < lend_row; ++lrow) | | | |
| { | | | |
| access::rw(m).at(lrow, lcol) += val; | | | |
| } | | | |
| | | | |
| const uword new_n_nonzero = m.n_nonzero; | | | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | return (*this).operator=( (*this) + tmp ); | |
| | | | |
| return *this; | | | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator-=(const eT val) | | SpSubview<eT>::operator-=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(val == eT(0)) | | if(val == eT(0)) | |
| { | | { | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
|
| const uword lstart_row = aux_row1; | | Mat<eT> tmp( (*this).n_rows, (*this).n_cols ); | |
| const uword lend_row = aux_row1 + n_rows; | | | |
| | | | |
|
| const uword lstart_col = aux_col1; | | tmp.fill(val); | |
| const uword lend_col = aux_col1 + n_cols; | | | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | return (*this).operator=( (*this) - tmp ); | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | inline | |
| | | const SpSubview<eT>& | |
| | | SpSubview<eT>::operator*=(const eT val) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
|
| for(uword lcol = lstart_col; lcol < lend_col; ++lcol) | | if(val != eT(0)) | |
| for(uword lrow = lstart_row; lrow < lend_row; ++lrow) | | | |
| { | | { | |
|
| access::rw(m).at(lrow, lcol) -= val; | | const uword lstart_row = aux_row1; | |
| } | | const uword lend_row = aux_row1 + n_rows; | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | const uword lstart_col = aux_col1; | |
| | | const uword lend_col = aux_col1 + n_cols; | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | for(uword c = lstart_col; c < lend_col; ++c) | |
| | | { | |
| | | for(uword r = m.col_ptrs[c]; r < m.col_ptrs[c + 1]; ++r) | |
| | | { | |
| | | if(m.row_indices[r] >= lstart_row && m.row_indices[r] < lend_row) | |
| | | { | |
| | | access::rw(m.values[r]) *= val; | |
| | | } | |
| | | } | |
| | | } | |
| | | } | |
| | | else | |
| | | { | |
| | | (*this).zeros(); | |
| | | } | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
|
| SpSubview<eT>::operator*=(const eT val) | | SpSubview<eT>::operator/=(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| if(val == eT(0)) | | arma_debug_check( (val == eT(0)), "element-wise division: division by zer | |
| { | | o" ); | |
| // Turn it all into zeros. | | | |
| for(iterator it(*this); it != end(); ++it) | | | |
| { | | | |
| (*it) = eT(0); // zero out the value. | | | |
| it.internal_pos--; | | | |
| } | | | |
| | | | |
| return *this; | | | |
| } | | | |
| | | | |
| const uword lstart_row = aux_row1; | | const uword lstart_row = aux_row1; | |
| const uword lend_row = aux_row1 + n_rows; | | const uword lend_row = aux_row1 + n_rows; | |
| | | | |
| const uword lstart_col = aux_col1; | | const uword lstart_col = aux_col1; | |
| const uword lend_col = aux_col1 + n_cols; | | const uword lend_col = aux_col1 + n_cols; | |
| | | | |
| for(uword c = lstart_col; c < lend_col; ++c) | | for(uword c = lstart_col; c < lend_col; ++c) | |
| { | | { | |
| for(uword r = m.col_ptrs[c]; r < m.col_ptrs[c + 1]; ++r) | | for(uword r = m.col_ptrs[c]; r < m.col_ptrs[c + 1]; ++r) | |
| { | | { | |
| if(m.row_indices[r] >= lstart_row && m.row_indices[r] < lend_row) | | if(m.row_indices[r] >= lstart_row && m.row_indices[r] < lend_row) | |
| { | | { | |
|
| access::rw(m.values[r]) *= val; | | access::rw(m.values[r]) /= val; | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
|
| | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
|
| SpSubview<eT>::operator/=(const eT val) | | SpSubview<eT>::operator=(const Base<eT, T1>& in) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const uword lstart_col = aux_col1; | | // this is a modified version of SpSubview::operator_equ_common(const SpB | |
| const uword lend_col = aux_col1 + n_cols; | | ase) | |
| | | | |
|
| const uword lstart_row = aux_row1; | | const SpProxy< SpMat<eT> > pa((*this).m); | |
| const uword lend_row = aux_row1 + n_rows; | | | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | const unwrap<T1> b_tmp(in.get_ref()); | |
| | | const Mat<eT>& b = b_tmp.M; | |
| | | | |
|
| for(uword c = lstart_col; c < lend_col; ++c) | | arma_debug_assert_same_size(n_rows, n_cols, b.n_rows, b.n_cols, "insertio | |
| for(uword r = lstart_row; r < lend_row; ++r) | | n into sparse submatrix"); | |
| | | | |
| | | const uword pa_start_row = (*this).aux_row1; | |
| | | const uword pa_start_col = (*this).aux_col1; | |
| | | | |
| | | const uword pa_end_row = pa_start_row + (*this).n_rows - 1; | |
| | | const uword pa_end_col = pa_start_col + (*this).n_cols - 1; | |
| | | | |
| | | const uword pa_n_rows = pa.get_n_rows(); | |
| | | | |
| | | const uword b_n_elem = b.n_elem; | |
| | | const eT* b_mem = b.memptr(); | |
| | | | |
| | | uword box_count = 0; | |
| | | | |
| | | for(uword i=0; i<b_n_elem; ++i) | |
| { | | { | |
|
| access::rw(m).at(r, c) /= val; | | if(b_mem[i] != eT(0)) { ++box_count; } | |
| } | | } | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | SpMat<eT> out(pa.get_n_rows(), pa.get_n_cols()); | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | const uword alt_count = pa.get_n_nonzero() - (*this).n_nonzero + box_coun
t; | |
| | | | |
|
| return *this; | | // Resize memory to correct size. | |
| } | | out.mem_resize(alt_count); | |
| | | | |
|
| template<typename eT> | | typename SpProxy< SpMat<eT> >::const_iterator_type x_it = pa.begin(); | |
| template<typename T1> | | typename SpProxy< SpMat<eT> >::const_iterator_type x_end = pa.end(); | |
| inline | | | |
| const SpSubview<eT>& | | | |
| SpSubview<eT>::operator=(const Base<eT, T1>& x) | | | |
| { | | | |
| arma_extra_debug_sigprint(); | | | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | uword b_row = 0; | |
| | | uword b_col = 0; | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols( | | bool x_it_ok = (x_it != x_end); | |
| ), "insert into sparse submatrix"); | | bool y_it_ok = ( (b_row < b.n_rows) && (b_col < b.n_cols) ); | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | uword x_it_row = (x_it_ok) ? x_it.row() : 0; | |
| | | uword x_it_col = (x_it_ok) ? x_it.col() : 0; | |
| | | | |
|
| for(uword c = 0; c < n_cols; ++c) | | uword y_it_row = (y_it_ok) ? b_row + pa_start_row : 0; | |
| | | uword y_it_col = (y_it_ok) ? b_col + pa_start_col : 0; | |
| | | | |
| | | uword cur_val = 0; | |
| | | while(x_it_ok || y_it_ok) | |
| { | | { | |
|
| for(uword r = 0; r < n_rows; ++r) | | const bool x_inside_box = (x_it_row >= pa_start_row) && (x_it_row <= pa | |
| | | _end_row) && (x_it_col >= pa_start_col) && (x_it_col <= pa_end_col); | |
| | | const bool y_inside_box = (y_it_row >= pa_start_row) && (y_it_row <= pa | |
| | | _end_row) && (y_it_col >= pa_start_col) && (y_it_col <= pa_end_col); | |
| | | | |
| | | const eT x_val = x_inside_box ? eT(0) : ( x_it_ok ? (*x_it) : eT(0) ); | |
| | | | |
| | | const eT y_val = y_inside_box ? ( y_it_ok ? b.at(b_row,b_col) : eT(0) ) | |
| | | : eT(0); | |
| | | | |
| | | if( (x_it_row == y_it_row) && (x_it_col == y_it_col) ) | |
| { | | { | |
|
| at(r, c) = P.at(r, c); | | if( (x_val != eT(0)) || (y_val != eT(0)) ) | |
| } | | { | |
| } | | access::rw(out.values[cur_val]) = (x_val != eT(0)) ? x_val : y_val; | |
| | | access::rw(out.row_indices[cur_val]) = x_it_row; | |
| | | ++access::rw(out.col_ptrs[x_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | if(x_it_ok) | |
| | | { | |
| | | ++x_it; | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | if(x_it == x_end) { x_it_ok = false; } | |
| | | } | |
| | | | |
|
| return *this; | | if(x_it_ok) | |
| } | | { | |
| | | x_it_row = x_it.row(); | |
| | | x_it_col = x_it.col(); | |
| | | } | |
| | | else | |
| | | { | |
| | | x_it_row++; | |
| | | | |
|
| template<typename eT> | | if(x_it_row >= pa_n_rows) { x_it_row = 0; x_it_col++; } | |
| template<typename T1> | | } | |
| inline | | | |
| const SpSubview<eT>& | | | |
| SpSubview<eT>::operator+=(const Base<eT, T1>& x) | | | |
| { | | | |
| arma_extra_debug_sigprint(); | | | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | if(y_it_ok) | |
| | | { | |
| | | b_row++; | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols(
), "addition"); | | if(b_row >= b.n_rows) { b_row = 0; b_col++; } | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | if( (b_row > b.n_rows) || (b_col > b.n_cols) ) { y_it_ok = false; | |
| | | } | |
| | | } | |
| | | | |
|
| for(uword c = 0; c < n_cols; ++c) | | if(y_it_ok) | |
| { | | { | |
| for(uword r = 0; r < n_rows; ++r) | | y_it_row = b_row + pa_start_row; | |
| | | y_it_col = b_col + pa_start_col; | |
| | | } | |
| | | else | |
| | | { | |
| | | y_it_row++; | |
| | | | |
| | | if(y_it_row >= pa_n_rows) { y_it_row = 0; y_it_col++; } | |
| | | } | |
| | | } | |
| | | else | |
| { | | { | |
|
| at(r, c) += P.at(r, c); | | 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 | |
| | | { | |
| | | if(x_val != eT(0)) | |
| | | { | |
| | | access::rw(out.values[cur_val]) = x_val; | |
| | | access::rw(out.row_indices[cur_val]) = x_it_row; | |
| | | ++access::rw(out.col_ptrs[x_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
| | | if(x_it_ok) | |
| | | { | |
| | | ++x_it; | |
| | | | |
| | | if(x_it == x_end) { x_it_ok = false; } | |
| | | } | |
| | | | |
| | | if(x_it_ok) | |
| | | { | |
| | | x_it_row = x_it.row(); | |
| | | x_it_col = x_it.col(); | |
| | | } | |
| | | else | |
| | | { | |
| | | x_it_row++; | |
| | | | |
| | | if(x_it_row >= pa_n_rows) { x_it_row = 0; x_it_col++; } | |
| | | } | |
| | | } | |
| | | else | |
| | | { | |
| | | if(y_val != eT(0)) | |
| | | { | |
| | | access::rw(out.values[cur_val]) = y_val; | |
| | | access::rw(out.row_indices[cur_val]) = y_it_row; | |
| | | ++access::rw(out.col_ptrs[y_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
| | | if(y_it_ok) | |
| | | { | |
| | | b_row++; | |
| | | | |
| | | if(b_row >= b.n_rows) { b_row = 0; b_col++; } | |
| | | | |
| | | if( (b_row > b.n_rows) || (b_col > b.n_cols) ) { y_it_ok = false | |
| | | ; } | |
| | | } | |
| | | | |
| | | if(y_it_ok) | |
| | | { | |
| | | y_it_row = b_row + pa_start_row; | |
| | | y_it_col = b_col + pa_start_col; | |
| | | } | |
| | | else | |
| | | { | |
| | | y_it_row++; | |
| | | | |
| | | if(y_it_row >= pa_n_rows) { y_it_row = 0; y_it_col++; } | |
| | | } | |
| | | } | |
| } | | } | |
| } | | } | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | const uword out_n_cols = out.n_cols; | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | uword* col_ptrs = access::rwp(out.col_ptrs); | |
| | | | |
| | | // Fix column pointers to be cumulative. | |
| | | for(uword c = 1; c <= out_n_cols; ++c) | |
| | | { | |
| | | col_ptrs[c] += col_ptrs[c - 1]; | |
| | | } | |
| | | | |
| | | access::rw((*this).m).steal_mem(out); | |
| | | | |
| | | access::rw(n_nonzero) = box_count; | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
|
| SpSubview<eT>::operator-=(const Base<eT, T1>& x) | | SpSubview<eT>::operator+=(const Base<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | return (*this).operator=( (*this) + x.get_ref() ); | |
| | | } | |
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols( | | | |
| ), "subtraction"); | | | |
| | | | |
| const uword old_n_nonzero = m.n_nonzero; | | | |
| | | | |
| for(uword c = 0; c < n_cols; ++c) | | | |
| { | | | |
| for(uword r = 0; r < n_rows; ++r) | | | |
| { | | | |
| at(r, c) -= P.at(r, c); | | | |
| } | | | |
| } | | | |
| | | | |
| const uword new_n_nonzero = m.n_nonzero; | | | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | template<typename eT> | |
| | | template<typename T1> | |
| | | inline | |
| | | const SpSubview<eT>& | |
| | | SpSubview<eT>::operator-=(const Base<eT, T1>& x) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
|
| return *this; | | return (*this).operator=( (*this) - x.get_ref() ); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator*=(const Base<eT, T1>& x) | | SpSubview<eT>::operator*=(const Base<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | | |
| | | | |
| // Must be exactly the same size for this (we can't modify our own size). | | | |
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols( | | | |
| ), "matrix multiplication"); | | | |
| | | | |
| SpMat<eT> tmp(*this); | | SpMat<eT> tmp(*this); | |
|
| Mat<eT> other_tmp(x.get_ref()); | | | |
| tmp *= other_tmp; | | | |
| operator=(tmp); | | | |
| | | | |
|
| return *this; | | tmp *= x.get_ref(); | |
| | | | |
| | | return (*this).operator=(tmp); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator%=(const Base<eT, T1>& x) | | SpSubview<eT>::operator%=(const Base<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | return (*this).operator=( (*this) % x.get_ref() ); | |
| | | } | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols( | | template<typename eT> | |
| ), "element-wise multiplication"); | | template<typename T1> | |
| | | inline | |
| | | const SpSubview<eT>& | |
| | | SpSubview<eT>::operator/=(const Base<eT, T1>& x) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | return (*this).operator=( (*this) / x.get_ref() ); | |
| | | } | |
| | | | |
|
| for(iterator it(*this); it != end(); ++it) | | template<typename eT> | |
| { | | inline | |
| (*it) *= P.at(it.row(), it.col()); | | const SpSubview<eT>& | |
| if(P.at(it.row(), it.col()) == eT(0)) | | SpSubview<eT>::operator=(const SpSubview<eT>& x) | |
| { | | { | |
| it.internal_pos--; | | arma_extra_debug_sigprint(); | |
| } | | | |
| } | | | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | return (*this).operator_equ_common(x); | |
| | | } | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | template<typename eT> | |
| | | template<typename T1> | |
| | | inline | |
| | | const SpSubview<eT>& | |
| | | SpSubview<eT>::operator=(const SpBase<eT, T1>& x) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
|
| return *this; | | return (*this).operator_equ_common( x.get_ref() ); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
|
| SpSubview<eT>::operator/=(const Base<eT, T1>& x) | | SpSubview<eT>::operator_equ_common(const SpBase<eT, T1>& in) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const Proxy<T1> P(x.get_ref()); | | // algorithm: | |
| | | // instead of directly inserting values into the matrix underlying the su | |
| | | bview, | |
| | | // create a new matrix by merging the underlying matrix with the input ob | |
| | | ject, | |
| | | // and then replacing the underlying matrix with the created matrix. | |
| | | // | |
| | | // the merging process requires pretending that the input object | |
| | | // has the same size as the underlying matrix. | |
| | | // while iterating through the elements of the input object, | |
| | | // this requires adjusting the row and column locations of each element, | |
| | | // as well as providing fake zero elements. | |
| | | // in effect there is a proxy for a proxy. | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols( | | const SpProxy< SpMat<eT> > pa((*this).m ); | |
| ), "element-wise division"); | | const SpProxy< T1 > pb(in.get_ref()); | |
| | | | |
|
| const uword old_n_nonzero = m.n_nonzero; | | arma_debug_assert_same_size(n_rows, n_cols, pb.get_n_rows(), pb.get_n_col
s(), "insertion into sparse submatrix"); | |
| | | | |
|
| for(uword c = 0; c < n_cols; ++c) | | const uword pa_start_row = (*this).aux_row1; | |
| { | | const uword pa_start_col = (*this).aux_col1; | |
| for(uword r = 0; r < n_rows; ++r) | | | |
| { | | | |
| at(r, c) /= P.at(r, c); | | | |
| } | | | |
| } | | | |
| | | | |
|
| const uword new_n_nonzero = m.n_nonzero; | | const uword pa_end_row = pa_start_row + (*this).n_rows - 1; | |
| | | const uword pa_end_col = pa_start_col + (*this).n_cols - 1; | |
| | | | |
|
| access::rw(n_nonzero) += (new_n_nonzero - old_n_nonzero); | | const uword pa_n_rows = pa.get_n_rows(); | |
| | | | |
|
| return *this; | | SpMat<eT> out(pa.get_n_rows(), pa.get_n_cols()); | |
| } | | | |
| | | | |
|
| template<typename eT> | | const uword alt_count = pa.get_n_nonzero() - (*this).n_nonzero + pb.get_n | |
| inline | | _nonzero(); | |
| const SpSubview<eT>& | | | |
| SpSubview<eT>::operator=(const SpSubview<eT>& x) | | | |
| { | | | |
| arma_extra_debug_sigprint(); | | | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "insertio | | // Resize memory to correct size. | |
| n into sparse submatrix"); | | out.mem_resize(alt_count); | |
| | | | |
|
| const bool alias = ( &m == &(x.m) ); | | typename SpProxy< SpMat<eT> >::const_iterator_type x_it = pa.begin(); | |
| | | typename SpProxy< SpMat<eT> >::const_iterator_type x_end = pa.end(); | |
| | | | |
|
| if(alias == false) | | typename SpProxy<T1>::const_iterator_type y_it = pb.begin(); | |
| | | typename SpProxy<T1>::const_iterator_type y_end = pb.end(); | |
| | | | |
| | | bool x_it_ok = (x_it != x_end); | |
| | | bool y_it_ok = (y_it != y_end); | |
| | | | |
| | | uword x_it_row = (x_it_ok) ? x_it.row() : 0; | |
| | | uword x_it_col = (x_it_ok) ? x_it.col() : 0; | |
| | | | |
| | | uword y_it_row = (y_it_ok) ? y_it.row() + pa_start_row : 0; | |
| | | uword y_it_col = (y_it_ok) ? y_it.col() + pa_start_col : 0; | |
| | | | |
| | | uword cur_val = 0; | |
| | | while(x_it_ok || y_it_ok) | |
| { | | { | |
|
| const_iterator cit = x.begin(); | | const bool x_inside_box = (x_it_row >= pa_start_row) && (x_it_row <= pa | |
| iterator it = begin(); | | _end_row) && (x_it_col >= pa_start_col) && (x_it_col <= pa_end_col); | |
| | | const bool y_inside_box = (y_it_row >= pa_start_row) && (y_it_row <= pa | |
| | | _end_row) && (y_it_col >= pa_start_col) && (y_it_col <= pa_end_col); | |
| | | | |
| | | const eT x_val = x_inside_box ? eT(0) : ( x_it_ok ? (*x_it) : eT(0) ); | |
| | | | |
|
| while((cit != x.end()) || (it != end())) | | const eT y_val = y_inside_box ? ( y_it_ok ? (*y_it) : eT(0) ) : eT(0); | |
| | | | |
| | | if( (x_it_row == y_it_row) && (x_it_col == y_it_col) ) | |
| { | | { | |
|
| if((cit.row() == it.row()) && (cit.col() == it.col())) | | if( (x_val != eT(0)) || (y_val != eT(0)) ) | |
| | | { | |
| | | access::rw(out.values[cur_val]) = (x_val != eT(0)) ? x_val : y_val; | |
| | | access::rw(out.row_indices[cur_val]) = x_it_row; | |
| | | ++access::rw(out.col_ptrs[x_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
| | | if(x_it_ok) | |
| { | | { | |
|
| (*it) = (*cit); | | ++x_it; | |
| ++it; | | | |
| ++cit; | | if(x_it == x_end) { x_it_ok = false; } | |
| | | } | |
| | | | |
| | | if(x_it_ok) | |
| | | { | |
| | | x_it_row = x_it.row(); | |
| | | x_it_col = x_it.col(); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| if((cit.col() > it.col()) || ((cit.col() == it.col()) && (cit.row() | | x_it_row++; | |
| > it.row()))) | | | |
| { | | if(x_it_row >= pa_n_rows) { x_it_row = 0; x_it_col++; } | |
| // cit is "ahead" | | | |
| (*it) = eT(0); // erase element | | | |
| it.internal_pos--; // update iterator so it still works | | | |
| ++it; | | | |
| } | | | |
| else | | | |
| { | | | |
| // it is "ahead" | | | |
| at(cit.row(), cit.col()) = (*cit); | | | |
| it.internal_pos++; // update iterator so it still works | | | |
| ++cit; | | | |
| } | | | |
| } | | } | |
|
| } | | | |
| | | | |
|
| access::rw(n_nonzero) = x.n_nonzero; | | if(y_it_ok) | |
| } | | { | |
| else | | ++y_it; | |
| { | | | |
| const SpMat<eT> tmp(x); | | | |
| | | | |
|
| (*this).operator=(tmp); | | if(y_it == y_end) { y_it_ok = false; } | |
| } | | } | |
| | | | |
|
| return *this; | | if(y_it_ok) | |
| } | | { | |
| | | y_it_row = y_it.row() + pa_start_row; | |
| | | y_it_col = y_it.col() + pa_start_col; | |
| | | } | |
| | | else | |
| | | { | |
| | | y_it_row++; | |
| | | | |
|
| template<typename eT> | | if(y_it_row >= pa_n_rows) { y_it_row = 0; y_it_col++; } | |
| template<typename T1> | | } | |
| inline | | } | |
| const SpSubview<eT>& | | else | |
| SpSubview<eT>::operator=(const SpBase<eT, T1>& x) | | { | |
| { | | if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_ | |
| arma_extra_debug_sigprint(); | | it_row))) // if y is closer to the end | |
| | | { | |
| | | if(x_val != eT(0)) | |
| | | { | |
| | | access::rw(out.values[cur_val]) = x_val; | |
| | | access::rw(out.row_indices[cur_val]) = x_it_row; | |
| | | ++access::rw(out.col_ptrs[x_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
|
| const SpProxy<T1> p(x.get_ref()); | | if(x_it_ok) | |
| | | { | |
| | | ++x_it; | |
| | | | |
|
| arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols( | | if(x_it == x_end) { x_it_ok = false; } | |
| ), "insertion into sparse submatrix"); | | } | |
| | | | |
|
| if(p.is_alias(m) == false) | | if(x_it_ok) | |
| { | | { | |
| typename SpProxy<T1>::const_iterator_type cit = p.begin(); | | x_it_row = x_it.row(); | |
| iterator it(*this); | | x_it_col = x_it.col(); | |
| | | } | |
| | | else | |
| | | { | |
| | | x_it_row++; | |
| | | | |
|
| while((cit != p.end()) || (it != end())) | | if(x_it_row >= pa_n_rows) { x_it_row = 0; x_it_col++; } | |
| { | | } | |
| if(cit == it) // at the same location | | | |
| { | | | |
| (*it) = (*cit); | | | |
| ++it; | | | |
| ++cit; | | | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| if((cit.col() > it.col()) || ((cit.col() == it.col()) && (cit.row() | | if(y_val != eT(0)) | |
| > it.row()))) | | { | |
| | | access::rw(out.values[cur_val]) = y_val; | |
| | | access::rw(out.row_indices[cur_val]) = y_it_row; | |
| | | ++access::rw(out.col_ptrs[y_it_col + 1]); | |
| | | ++cur_val; | |
| | | } | |
| | | | |
| | | if(y_it_ok) | |
| { | | { | |
|
| // cit is "ahead" | | ++y_it; | |
| (*it) = eT(0); // erase element | | | |
| it.internal_pos--; // update iterator so it still works | | if(y_it == y_end) { y_it_ok = false; } | |
| ++it; | | } | |
| | | | |
| | | if(y_it_ok) | |
| | | { | |
| | | y_it_row = y_it.row() + pa_start_row; | |
| | | y_it_col = y_it.col() + pa_start_col; | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| // it is "ahead" | | y_it_row++; | |
| at(cit.row(), cit.col()) = (*cit); | | | |
| it.internal_pos++; // update iterator so it still works | | if(y_it_row >= pa_n_rows) { y_it_row = 0; y_it_col++; } | |
| ++cit; | | | |
| } | | } | |
| } | | } | |
| } | | } | |
| } | | } | |
|
| else | | | |
| { | | | |
| const SpMat<eT> tmp(p.Q); | | | |
| | | | |
|
| (*this).operator=(tmp); | | const uword out_n_cols = out.n_cols; | |
| | | | |
| | | uword* col_ptrs = access::rwp(out.col_ptrs); | |
| | | | |
| | | // Fix column pointers to be cumulative. | |
| | | for(uword c = 1; c <= out_n_cols; ++c) | |
| | | { | |
| | | col_ptrs[c] += col_ptrs[c - 1]; | |
| } | | } | |
| | | | |
|
| | | access::rw((*this).m).steal_mem(out); | |
| | | | |
| | | access::rw(n_nonzero) = pb.get_n_nonzero(); | |
| | | | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator+=(const SpBase<eT, T1>& x) | | SpSubview<eT>::operator+=(const SpBase<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const SpProxy<T1> p(x.get_ref()); | | // TODO: implement dedicated machinery | |
| | | return (*this).operator=( (*this) + x.get_ref() ); | |
| arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols( | | | |
| ), "addition"); | | | |
| | | | |
| if(p.is_alias(m) == false) | | | |
| { | | | |
| typename SpProxy<T1>::const_iterator_type cit = p.begin(); | | | |
| | | | |
| while(cit != p.end()) | | | |
| { | | | |
| at(cit.row(), cit.col()) += (*cit); | | | |
| ++cit; | | | |
| } | | | |
| } | | | |
| else | | | |
| { | | | |
| const SpMat<eT> tmp(p.Q); | | | |
| | | | |
| (*this).operator+=(tmp); | | | |
| } | | | |
| | | | |
| return *this; | | | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator-=(const SpBase<eT, T1>& x) | | SpSubview<eT>::operator-=(const SpBase<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const SpProxy<T1> p(x.get_ref()); | | // TODO: implement dedicated machinery | |
| | | return (*this).operator=( (*this) - x.get_ref() ); | |
| arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols( | | | |
| ), "subtraction"); | | | |
| | | | |
| if(p.is_alias(m) == false) | | | |
| { | | | |
| typename SpProxy<T1>::const_iterator_type cit = p.begin(); | | | |
| | | | |
| while(cit != p.end()) | | | |
| { | | | |
| at(cit.row(), cit.col()) -= (*cit); | | | |
| ++cit; | | | |
| } | | | |
| } | | | |
| else | | | |
| { | | | |
| const SpMat<eT> tmp(p.Q); | | | |
| | | | |
| (*this).operator-=(tmp); | | | |
| } | | | |
| | | | |
| return *this; | | | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator*=(const SpBase<eT, T1>& x) | | SpSubview<eT>::operator*=(const SpBase<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const SpProxy<T1> p(x.get_ref()); | | return (*this).operator=( (*this) * x.get_ref() ); | |
| | | | |
| arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols( | | | |
| ), "matrix multiplication"); | | | |
| | | | |
| // Because we have to use a temporary anyway, it doesn't make sense to | | | |
| // reimplement this here. | | | |
| return operator=((*this) * x.get_ref()); | | | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator%=(const SpBase<eT, T1>& x) | | SpSubview<eT>::operator%=(const SpBase<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| const SpProxy<T1> p(x.get_ref()); | | // TODO: implement dedicated machinery | |
| | | return (*this).operator=( (*this) % x.get_ref() ); | |
| arma_debug_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols( | | | |
| ), "element-wise multiplication"); | | | |
| | | | |
| if(p.is_alias(m) == false) | | | |
| { | | | |
| typename SpProxy<T1>::const_iterator_type cit = p.begin(); | | | |
| iterator it(*this); | | | |
| | | | |
| while((it != end()) || (cit != p.end())) | | | |
| { | | | |
| if((cit.row() == it.row()) && (cit.col() == it.col())) | | | |
| { | | | |
| (*it) *= (*cit); | | | |
| ++it; | | | |
| ++cit; | | | |
| } | | | |
| else | | | |
| { | | | |
| if((cit.col() > it.col()) || ((cit.col() == it.col()) && (cit.row() | | | |
| > it.row()))) | | | |
| { | | | |
| // cit is "ahead" | | | |
| (*it) = eT(0); // erase element -- x has a zero here | | | |
| it.internal_pos--; // update iterator so it still works | | | |
| ++it; | | | |
| } | | | |
| else | | | |
| { | | | |
| // it is "ahead" | | | |
| ++cit; | | | |
| } | | | |
| } | | | |
| } | | | |
| } | | | |
| else | | | |
| { | | | |
| const SpMat<eT> tmp(p.Q); | | | |
| | | | |
| (*this).operator%=(tmp); | | | |
| } | | | |
| | | | |
| return *this; | | | |
| } | | } | |
| | | | |
| //! If you are using this function, you are probably misguided. | | //! If you are using this function, you are probably misguided. | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const SpSubview<eT>& | | const SpSubview<eT>& | |
| SpSubview<eT>::operator/=(const SpBase<eT, T1>& x) | | SpSubview<eT>::operator/=(const SpBase<eT, T1>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| skipping to change at line 660 | | skipping to change at line 737 | |
|---|
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
| SpSubview<eT>::fill(const eT val) | | SpSubview<eT>::fill(const eT val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(val != eT(0)) | | if(val != eT(0)) | |
| { | | { | |
|
| // TODO: implement a faster version; the code below is slow | | Mat<eT> tmp( (*this).n_rows, (*this).n_cols ); | |
| | | | |
| const uword lstart_row = aux_row1; | | | |
| const uword lend_row = aux_row1 + n_rows; | | | |
| | | | |
| const uword lstart_col = aux_col1; | | | |
| const uword lend_col = aux_col1 + n_cols; | | | |
| | | | |
|
| const uword orig_nonzero = m.n_nonzero; | | tmp.fill(val); | |
| | | | |
| // iterate over our part of the sparse matrix | | | |
| for(uword lcol = lstart_col; lcol < lend_col; ++lcol) | | | |
| for(uword lrow = lstart_row; lrow < lend_row; ++lrow) | | | |
| { | | | |
| access::rw(m).at(lrow, lcol) = val; | | | |
| } | | | |
| | | | |
| access::rw(n_nonzero) += (m.n_nonzero - orig_nonzero); | | | |
| | | | |
|
| | | (*this).operator=(tmp); | |
| } | | } | |
| else | | else | |
| { | | { | |
| (*this).zeros(); | | (*this).zeros(); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
| SpSubview<eT>::zeros() | | SpSubview<eT>::zeros() | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| // we can be a little smarter here | | const SpMat<eT> tmp( (*this).n_rows, (*this).n_cols ); | |
| iterator it(*this); | | | |
| | | | |
|
| while(it != end()) | | (*this).operator=(tmp); | |
| { | | | |
| (*it) = eT(0); | | | |
| it.internal_pos--; // hack to update iterator without requiring a new o | | | |
| ne | | | |
| ++it; | | | |
| } | | | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
| SpSubview<eT>::ones() | | SpSubview<eT>::ones() | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| (*this).fill(eT(1)); | | (*this).fill(eT(1)); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
| SpSubview<eT>::eye() | | SpSubview<eT>::eye() | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| // clear other things | | SpMat<eT> tmp; | |
| (*this).zeros(); | | | |
| | | | |
| const uword orig_nonzero = m.n_nonzero; | | | |
| | | | |
| // now the diagonal ones | | | |
| const uword end_index = std::min(n_rows, n_cols); | | | |
| | | | |
|
| for(uword ind = 0; ind < end_index; ++ind) | | tmp.eye( (*this).n_rows, (*this).n_cols ); | |
| { | | | |
| access::rw(m).at(ind + aux_row1, ind + aux_col1) = eT(1); | | | |
| } | | | |
| | | | |
|
| access::rw(n_nonzero) += (m.n_nonzero - orig_nonzero); | | (*this).operator=(tmp); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| SpValProxy<SpSubview<eT> > | | SpValProxy<SpSubview<eT> > | |
| SpSubview<eT>::operator[](const uword i) | | SpSubview<eT>::operator[](const uword i) | |
| { | | { | |
| const uword lrow = i % n_rows; | | const uword lrow = i % n_rows; | |
| const uword lcol = i / n_rows; | | const uword lcol = i / n_rows; | |
| | | | |
| End of changes. 101 change blocks. |
|---|
| 365 lines changed or deleted | | 413 lines changed or added | |
|---|
|
| arrayops_meat.hpp | | arrayops_meat.hpp | |
|---|
| | | | |
| skipping to change at line 142 | | skipping to change at line 142 | |
|---|
| ; | | ; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::fill_zeros(eT* dest, const uword n_elem) | | arrayops::fill_zeros(eT* dest, const uword n_elem) | |
| { | | { | |
|
| typedef typename get_pod_type<eT>::result pod_type; | | arrayops::inplace_set(dest, eT(0), n_elem); | |
| | | | |
| if( (n_elem >= 16) && (std::numeric_limits<eT>::is_integer || (std::numer | | | |
| ic_limits<pod_type>::is_iec559 && is_real<pod_type>::value)) ) | | | |
| { | | | |
| std::memset(dest, 0, sizeof(eT)*n_elem); | | | |
| } | | | |
| else | | | |
| { | | | |
| arrayops::inplace_set(dest, eT(0), n_elem); | | | |
| } | | | |
| } | | } | |
| | | | |
| template<typename out_eT, typename in_eT> | | template<typename out_eT, typename in_eT> | |
| arma_hot | | arma_hot | |
| arma_inline | | arma_inline | |
| void | | void | |
| arrayops::convert_cx_scalar | | arrayops::convert_cx_scalar | |
| ( | | ( | |
| out_eT& out, | | out_eT& out, | |
| const in_eT& in, | | const in_eT& in, | |
| | | | |
| skipping to change at line 410 | | skipping to change at line 401 | |
|---|
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_plus_base(eT* dest, const eT* src, const uword n_elem) | | arrayops::inplace_plus_base(eT* dest, const eT* src, const uword n_elem) | |
| { | | { | |
|
| uword i,j; | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | { | |
|
| const eT tmp_i = src[i]; | | for(uword i=0; i<n_elem; ++i) | |
| const eT tmp_j = src[j]; | | { | |
| | | dest[i] += src[i]; | |
| dest[i] += tmp_i; | | } | |
| dest[j] += tmp_j; | | | |
| } | | } | |
|
| | | #else | |
| if(i < n_elem) | | | |
| { | | { | |
|
| dest[i] += src[i]; | | uword i,j; | |
| | | | |
| | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| | | { | |
| | | const eT tmp_i = src[i]; | |
| | | const eT tmp_j = src[j]; | |
| | | | |
| | | dest[i] += tmp_i; | |
| | | dest[j] += tmp_j; | |
| | | } | |
| | | | |
| | | if(i < n_elem) | |
| | | { | |
| | | dest[i] += src[i]; | |
| | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_minus_base(eT* dest, const eT* src, const uword n_elem) | | arrayops::inplace_minus_base(eT* dest, const eT* src, const uword n_elem) | |
| { | | { | |
|
| uword i,j; | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | { | |
|
| const eT tmp_i = src[i]; | | for(uword i=0; i<n_elem; ++i) | |
| const eT tmp_j = src[j]; | | { | |
| | | dest[i] -= src[i]; | |
| dest[i] -= tmp_i; | | } | |
| dest[j] -= tmp_j; | | | |
| } | | } | |
|
| | | #else | |
| if(i < n_elem) | | | |
| { | | { | |
|
| dest[i] -= src[i]; | | uword i,j; | |
| | | | |
| | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| | | { | |
| | | const eT tmp_i = src[i]; | |
| | | const eT tmp_j = src[j]; | |
| | | | |
| | | dest[i] -= tmp_i; | |
| | | dest[j] -= tmp_j; | |
| | | } | |
| | | | |
| | | if(i < n_elem) | |
| | | { | |
| | | dest[i] -= src[i]; | |
| | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_mul_base(eT* dest, const eT* src, const uword n_elem) | | arrayops::inplace_mul_base(eT* dest, const eT* src, const uword n_elem) | |
| { | | { | |
|
| uword i,j; | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | { | |
|
| const eT tmp_i = src[i]; | | for(uword i=0; i<n_elem; ++i) | |
| const eT tmp_j = src[j]; | | { | |
| | | dest[i] *= src[i]; | |
| dest[i] *= tmp_i; | | } | |
| dest[j] *= tmp_j; | | | |
| } | | } | |
|
| | | #else | |
| if(i < n_elem) | | | |
| { | | { | |
|
| dest[i] *= src[i]; | | uword i,j; | |
| | | | |
| | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| | | { | |
| | | const eT tmp_i = src[i]; | |
| | | const eT tmp_j = src[j]; | |
| | | | |
| | | dest[i] *= tmp_i; | |
| | | dest[j] *= tmp_j; | |
| | | } | |
| | | | |
| | | if(i < n_elem) | |
| | | { | |
| | | dest[i] *= src[i]; | |
| | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_div_base(eT* dest, const eT* src, const uword n_elem) | | arrayops::inplace_div_base(eT* dest, const eT* src, const uword n_elem) | |
| { | | { | |
|
| uword i,j; | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | { | |
|
| const eT tmp_i = src[i]; | | for(uword i=0; i<n_elem; ++i) | |
| const eT tmp_j = src[j]; | | { | |
| | | dest[i] /= src[i]; | |
| dest[i] /= tmp_i; | | } | |
| dest[j] /= tmp_j; | | | |
| } | | } | |
|
| | | #else | |
| if(i < n_elem) | | | |
| { | | { | |
|
| dest[i] /= src[i]; | | uword i,j; | |
| | | | |
| | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| | | { | |
| | | const eT tmp_i = src[i]; | |
| | | const eT tmp_j = src[j]; | |
| | | | |
| | | dest[i] /= tmp_i; | |
| | | dest[j] /= tmp_j; | |
| | | } | |
| | | | |
| | | if(i < n_elem) | |
| | | { | |
| | | dest[i] /= src[i]; | |
| | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_set(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_set(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
|
| | | typedef typename get_pod_type<eT>::result pod_type; | |
| | | | |
| if( (n_elem <= 16) && (is_cx<eT>::no) ) | | if( (n_elem <= 16) && (is_cx<eT>::no) ) | |
| { | | { | |
| arrayops::inplace_set_small(dest, val, n_elem); | | arrayops::inplace_set_small(dest, val, n_elem); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| if(memory::is_aligned(dest)) | | if( (val == eT(0)) && (std::numeric_limits<eT>::is_integer || (std::num
eric_limits<pod_type>::is_iec559 && is_real<pod_type>::value)) ) | |
| { | | { | |
|
| memory::mark_as_aligned(dest); | | std::memset(dest, 0, sizeof(eT)*n_elem); | |
| | | } | |
| uword i,j; | | else | |
| | | { | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | if(memory::is_aligned(dest)) | |
| { | | { | |
|
| dest[i] = val; | | memory::mark_as_aligned(dest); | |
| dest[j] = val; | | | |
| } | | | |
| | | | |
|
| if(i < n_elem) | | arrayops::inplace_set_base(dest, val, n_elem); | |
| | | } | |
| | | else | |
| { | | { | |
|
| dest[i] = val; | | arrayops::inplace_set_base(dest, val, n_elem); | |
| } | | } | |
| } | | } | |
|
| else | | } | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | arma_hot | |
| | | inline | |
| | | void | |
| | | arrayops::inplace_set_base(eT* dest, const eT val, const uword n_elem) | |
| | | { | |
| | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | { | |
| | | for(uword i=0; i<n_elem; ++i) | |
| { | | { | |
|
| uword i,j; | | dest[i] = val; | |
| | | } | |
| | | } | |
| | | #else | |
| | | { | |
| | | uword i,j; | |
| | | | |
|
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| dest[i] = val; | | dest[i] = val; | |
| dest[j] = val; | | dest[j] = val; | |
| } | | } | |
| | | | |
|
| if(i < n_elem) | | if(i < n_elem) | |
| { | | { | |
| dest[i] = val; | | dest[i] = val; | |
| } | | | |
| } | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_set_small(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_set_small(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
| switch(n_elem) | | switch(n_elem) | |
| { | | { | |
| | | | |
| skipping to change at line 593 | | skipping to change at line 646 | |
|---|
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
| arrayops::inplace_plus(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_plus(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
| if(memory::is_aligned(dest)) | | if(memory::is_aligned(dest)) | |
| { | | { | |
| memory::mark_as_aligned(dest); | | memory::mark_as_aligned(dest); | |
| | | | |
|
| uword i,j; | | arrayops::inplace_plus_base(dest, val, n_elem); | |
| | | } | |
| | | else | |
| | | { | |
| | | arrayops::inplace_plus_base(dest, val, n_elem); | |
| | | } | |
| | | } | |
| | | | |
|
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | template<typename eT> | |
| { | | arma_hot | |
| dest[i] += val; | | inline | |
| dest[j] += val; | | void | |
| } | | arrayops::inplace_minus(eT* dest, const eT val, const uword n_elem) | |
| | | { | |
| | | if(memory::is_aligned(dest)) | |
| | | { | |
| | | memory::mark_as_aligned(dest); | |
| | | | |
|
| if(i < n_elem) | | arrayops::inplace_minus_base(dest, val, n_elem); | |
| | | } | |
| | | else | |
| | | { | |
| | | arrayops::inplace_minus_base(dest, val, n_elem); | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | arma_hot | |
| | | inline | |
| | | void | |
| | | arrayops::inplace_mul(eT* dest, const eT val, const uword n_elem) | |
| | | { | |
| | | if(memory::is_aligned(dest)) | |
| | | { | |
| | | memory::mark_as_aligned(dest); | |
| | | | |
| | | arrayops::inplace_mul_base(dest, val, n_elem); | |
| | | } | |
| | | else | |
| | | { | |
| | | arrayops::inplace_mul_base(dest, val, n_elem); | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | arma_hot | |
| | | inline | |
| | | void | |
| | | arrayops::inplace_div(eT* dest, const eT val, const uword n_elem) | |
| | | { | |
| | | if(memory::is_aligned(dest)) | |
| | | { | |
| | | memory::mark_as_aligned(dest); | |
| | | | |
| | | arrayops::inplace_div_base(dest, val, n_elem); | |
| | | } | |
| | | else | |
| | | { | |
| | | arrayops::inplace_div_base(dest, val, n_elem); | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | arma_hot | |
| | | inline | |
| | | void | |
| | | arrayops::inplace_plus_base(eT* dest, const eT val, const uword n_elem) | |
| | | { | |
| | | #if defined(ARMA_SIMPLE_LOOPS) | |
| | | { | |
| | | for(uword i=0; i<n_elem; ++i) | |
| { | | { | |
| dest[i] += val; | | dest[i] += val; | |
| } | | } | |
| } | | } | |
|
| else | | #else | |
| { | | { | |
| uword i,j; | | uword i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| dest[i] += val; | | dest[i] += val; | |
| dest[j] += val; | | dest[j] += val; | |
| } | | } | |
| | | | |
| if(i < n_elem) | | if(i < n_elem) | |
| { | | { | |
| dest[i] += val; | | dest[i] += val; | |
| } | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
|
| arrayops::inplace_minus(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_minus_base(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
|
| if(memory::is_aligned(dest)) | | #if defined(ARMA_SIMPLE_LOOPS) | |
| { | | { | |
|
| memory::mark_as_aligned(dest); | | for(uword i=0; i<n_elem; ++i) | |
| | | | |
| uword i,j; | | | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | | |
| dest[i] -= val; | | | |
| dest[j] -= val; | | | |
| } | | | |
| | | | |
| if(i < n_elem) | | | |
| { | | { | |
| dest[i] -= val; | | dest[i] -= val; | |
| } | | } | |
| } | | } | |
|
| else | | #else | |
| { | | { | |
| uword i,j; | | uword i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| dest[i] -= val; | | dest[i] -= val; | |
| dest[j] -= val; | | dest[j] -= val; | |
| } | | } | |
| | | | |
| if(i < n_elem) | | if(i < n_elem) | |
| { | | { | |
| dest[i] -= val; | | dest[i] -= val; | |
| } | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
|
| arrayops::inplace_mul(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_mul_base(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
|
| if(memory::is_aligned(dest)) | | #if defined(ARMA_SIMPLE_LOOPS) | |
| { | | { | |
|
| memory::mark_as_aligned(dest); | | for(uword i=0; i<n_elem; ++i) | |
| | | | |
| uword i,j; | | | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | | |
| dest[i] *= val; | | | |
| dest[j] *= val; | | | |
| } | | | |
| | | | |
| if(i < n_elem) | | | |
| { | | { | |
| dest[i] *= val; | | dest[i] *= val; | |
| } | | } | |
| } | | } | |
|
| else | | #else | |
| { | | { | |
| uword i,j; | | uword i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| dest[i] *= val; | | dest[i] *= val; | |
| dest[j] *= val; | | dest[j] *= val; | |
| } | | } | |
| | | | |
| if(i < n_elem) | | if(i < n_elem) | |
| { | | { | |
| dest[i] *= val; | | dest[i] *= val; | |
| } | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| inline | | inline | |
| void | | void | |
|
| arrayops::inplace_div(eT* dest, const eT val, const uword n_elem) | | arrayops::inplace_div_base(eT* dest, const eT val, const uword n_elem) | |
| { | | { | |
|
| if(memory::is_aligned(dest)) | | #if defined(ARMA_SIMPLE_LOOPS) | |
| { | | { | |
|
| memory::mark_as_aligned(dest); | | for(uword i=0; i<n_elem; ++i) | |
| | | | |
| uword i,j; | | | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | | |
| { | | | |
| dest[i] /= val; | | | |
| dest[j] /= val; | | | |
| } | | | |
| | | | |
| if(i < n_elem) | | | |
| { | | { | |
| dest[i] /= val; | | dest[i] /= val; | |
| } | | } | |
| } | | } | |
|
| else | | #else | |
| { | | { | |
| uword i,j; | | uword i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| dest[i] /= val; | | dest[i] /= val; | |
| dest[j] /= val; | | dest[j] /= val; | |
| } | | } | |
| | | | |
| if(i < n_elem) | | if(i < n_elem) | |
| { | | { | |
| dest[i] /= val; | | dest[i] /= val; | |
| } | | } | |
| } | | } | |
|
| | | #endif | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_hot | | arma_hot | |
| arma_pure | | arma_pure | |
| inline | | inline | |
| eT | | eT | |
| arrayops::accumulate(const eT* src, const uword n_elem) | | arrayops::accumulate(const eT* src, const uword n_elem) | |
| { | | { | |
| uword i,j; | | uword i,j; | |
| | | | |
| End of changes. 52 change blocks. |
|---|
| 127 lines changed or deleted | | 215 lines changed or added | |
|---|
|
| field_bones.hpp | | field_bones.hpp | |
|---|
|
| // Copyright (C) 2008-2013 Conrad Sanderson | | // Copyright (C) 2008-2014 Conrad Sanderson | |
| // Copyright (C) 2008-2013 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2014 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2009-2010 Ian Cullinan | | // Copyright (C) 2009-2010 Ian Cullinan | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup field | | //! \addtogroup field | |
| //! @{ | | //! @{ | |
| | | | |
| struct field_prealloc_n_elem | | struct field_prealloc_n_elem | |
| | | | |
| skipping to change at line 29 | | skipping to change at line 29 | |
|---|
| | | | |
| template<typename oT> | | template<typename oT> | |
| class field | | class field | |
| { | | { | |
| public: | | public: | |
| | | | |
| typedef oT object_type; | | typedef oT object_type; | |
| | | | |
| const uword n_rows; //!< number of rows in the field (read-only) | | const uword n_rows; //!< number of rows in the field (read-only) | |
| const uword n_cols; //!< number of columns in the field (read-only) | | const uword n_cols; //!< number of columns in the field (read-only) | |
|
| | | const uword n_slices; //!< number of slices in the field (read-only) | |
| const uword n_elem; //!< number of elements in the field (read-only) | | const uword n_elem; //!< number of elements in the field (read-only) | |
| | | | |
| private: | | private: | |
| | | | |
| arma_aligned oT** mem; //!< pointer t
o memory used by the object | | arma_aligned oT** mem; //!< pointer t
o memory used by the object | |
| arma_aligned oT* mem_local[ field_prealloc_n_elem::val ]; //!< Internal
memory, to avoid calling the 'new' operator for small amounts of memory | | arma_aligned oT* mem_local[ field_prealloc_n_elem::val ]; //!< Internal
memory, to avoid calling the 'new' operator for small amounts of memory | |
| | | | |
| public: | | public: | |
| | | | |
| inline ~field(); | | inline ~field(); | |
| inline field(); | | inline field(); | |
| | | | |
| inline field(const field& x); | | inline field(const field& x); | |
| inline const field& operator=(const field& x); | | inline const field& operator=(const field& x); | |
| | | | |
| inline field(const subview_field<oT>& x); | | inline field(const subview_field<oT>& x); | |
| inline const field& operator=(const subview_field<oT>& x); | | inline const field& operator=(const subview_field<oT>& x); | |
| | | | |
| inline explicit field(const uword n_elem_in); | | inline explicit field(const uword n_elem_in); | |
| inline field(const uword n_rows_in, const uword n_cols_in); | | inline field(const uword n_rows_in, const uword n_cols_in); | |
|
| | | inline field(const uword n_rows_in, const uword n_cols_in, const
uword n_slices_in); | |
| | | | |
| inline void set_size(const uword n_obj_in); | | inline void set_size(const uword n_obj_in); | |
| inline void set_size(const uword n_rows_in, const uword n_cols_in); | | inline void set_size(const uword n_rows_in, const uword n_cols_in); | |
|
| | | inline void set_size(const uword n_rows_in, const uword n_cols_in, const
uword n_slices_in); | |
| | | | |
| template<typename oT2> | | template<typename oT2> | |
| inline void copy_size(const field<oT2>& x); | | inline void copy_size(const field<oT2>& x); | |
| | | | |
| arma_inline oT& operator[](const uword i); | | arma_inline oT& operator[](const uword i); | |
| arma_inline const oT& operator[](const uword i) const; | | arma_inline const oT& operator[](const uword i) const; | |
| | | | |
| arma_inline oT& at(const uword i); | | arma_inline oT& at(const uword i); | |
| arma_inline const oT& at(const uword i) const; | | arma_inline const oT& at(const uword i) const; | |
| | | | |
| arma_inline oT& operator()(const uword i); | | arma_inline oT& operator()(const uword i); | |
| arma_inline const oT& operator()(const uword i) const; | | arma_inline const oT& operator()(const uword i) const; | |
| | | | |
| arma_inline oT& at(const uword row, const uword col); | | arma_inline oT& at(const uword row, const uword col); | |
| arma_inline const oT& at(const uword row, const uword col) const; | | arma_inline const oT& at(const uword row, const uword col) const; | |
| | | | |
|
| | | arma_inline oT& at(const uword row, const uword col, const | |
| | | uword slice); | |
| | | arma_inline const oT& at(const uword row, const uword col, const | |
| | | uword slice) const; | |
| | | | |
| arma_inline oT& operator()(const uword row, const uword col); | | arma_inline oT& operator()(const uword row, const uword col); | |
| arma_inline const oT& operator()(const uword row, const uword col) const; | | arma_inline const oT& operator()(const uword row, const uword col) const; | |
| | | | |
|
| | | arma_inline oT& operator()(const uword row, const uword col, const | |
| | | uword slice); | |
| | | arma_inline const oT& operator()(const uword row, const uword col, const | |
| | | uword slice) const; | |
| | | | |
| inline field_injector<field> operator<<(const oT& val); | | inline field_injector<field> operator<<(const oT& val); | |
| inline field_injector<field> operator<<(const injector_end_of_row<>& x); | | inline field_injector<field> operator<<(const injector_end_of_row<>& x); | |
| | | | |
| inline subview_field<oT> row(const uword row_num); | | inline subview_field<oT> row(const uword row_num); | |
| inline const subview_field<oT> row(const uword row_num) const; | | inline const subview_field<oT> row(const uword row_num) const; | |
| | | | |
| inline subview_field<oT> col(const uword col_num); | | inline subview_field<oT> col(const uword col_num); | |
| inline const subview_field<oT> col(const uword col_num) const; | | inline const subview_field<oT> col(const uword col_num) const; | |
| | | | |
|
| | | inline subview_field<oT> slice(const uword slice_num); | |
| | | inline const subview_field<oT> slice(const uword slice_num) const; | |
| | | | |
| inline subview_field<oT> rows(const uword in_row1, const uword in_r
ow2); | | inline subview_field<oT> rows(const uword in_row1, const uword in_r
ow2); | |
| inline const subview_field<oT> rows(const uword in_row1, const uword in_r
ow2) const; | | inline const subview_field<oT> rows(const uword in_row1, const uword in_r
ow2) const; | |
| | | | |
| inline subview_field<oT> cols(const uword in_col1, const uword in_c
ol2); | | inline subview_field<oT> cols(const uword in_col1, const uword in_c
ol2); | |
| inline const subview_field<oT> cols(const uword in_col1, const uword in_c
ol2) const; | | inline const subview_field<oT> cols(const uword in_col1, const uword in_c
ol2) const; | |
| | | | |
|
| | | inline subview_field<oT> slices(const uword in_slice1, const uword | |
| | | in_slice2); | |
| | | inline const subview_field<oT> slices(const uword in_slice1, const uword | |
| | | in_slice2) const; | |
| | | | |
| inline subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const uword in_row2, const uword in_col2); | | inline subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const uword in_row2, const uword in_col2); | |
| inline const subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const uword in_row2, const uword in_col2) const; | | inline const subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const uword in_row2, const uword in_col2) const; | |
| | | | |
|
| | | inline subview_field<oT> subfield(const uword in_row1, const uword | |
| | | in_col1, const uword in_slice1, const uword in_row2, const uword in_col2, c | |
| | | onst uword in_slice2); | |
| | | inline const subview_field<oT> subfield(const uword in_row1, const uword | |
| | | in_col1, const uword in_slice1, const uword in_row2, const uword in_col2, c | |
| | | onst uword in_slice2) const; | |
| | | | |
| inline subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const SizeMat& s); | | inline subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const SizeMat& s); | |
| inline const subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const SizeMat& s) const; | | inline const subview_field<oT> subfield(const uword in_row1, const uword
in_col1, const SizeMat& s) const; | |
| | | | |
|
| inline subview_field<oT> subfield (const span& row_span, const spa | | inline subview_field<oT> subfield(const uword in_row1, const uword | |
| n& col_span); | | in_col1, const uword in_slice1, const SizeCube& s); | |
| inline const subview_field<oT> subfield (const span& row_span, const spa | | inline const subview_field<oT> subfield(const uword in_row1, const uword | |
| n& col_span) const; | | in_col1, const uword in_slice1, const SizeCube& s) const; | |
| | | | |
| | | inline subview_field<oT> subfield(const span& row_span, const span& | |
| | | col_span); | |
| | | inline const subview_field<oT> subfield(const span& row_span, const span& | |
| | | col_span) const; | |
| | | | |
| | | inline subview_field<oT> subfield(const span& row_span, const span& | |
| | | col_span, const span& slice_span); | |
| | | inline const subview_field<oT> subfield(const span& row_span, const span& | |
| | | col_span, const span& slice_span) const; | |
| | | | |
| inline subview_field<oT> operator()(const span& row_span, const spa
n& col_span); | | inline subview_field<oT> operator()(const span& row_span, const spa
n& col_span); | |
| inline const subview_field<oT> operator()(const span& row_span, const spa
n& col_span) const; | | inline const subview_field<oT> operator()(const span& row_span, const spa
n& col_span) const; | |
| | | | |
|
| | | inline subview_field<oT> operator()(const span& row_span, const spa | |
| | | n& col_span, const span& slice_span); | |
| | | inline const subview_field<oT> operator()(const span& row_span, const spa | |
| | | n& col_span, const span& slice_span) const; | |
| | | | |
| inline subview_field<oT> operator()(const uword in_row1, const uwor
d in_col1, const SizeMat& s); | | inline subview_field<oT> operator()(const uword in_row1, const uwor
d in_col1, const SizeMat& s); | |
| inline const subview_field<oT> operator()(const uword in_row1, const uwor
d in_col1, const SizeMat& s) const; | | inline const subview_field<oT> operator()(const uword in_row1, const uwor
d in_col1, const SizeMat& s) const; | |
| | | | |
|
| | | inline subview_field<oT> operator()(const uword in_row1, const uwor | |
| | | d in_col1, const uword in_slice1, const SizeCube& s); | |
| | | inline const subview_field<oT> operator()(const uword in_row1, const uwor | |
| | | d in_col1, const uword in_slice1, const SizeCube& s) const; | |
| | | | |
| inline void print(const std::string extra_text = "") const; | | inline void print(const std::string extra_text = "") const; | |
| inline void print(std::ostream& user_stream, const std::string extra_text
= "") const; | | inline void print(std::ostream& user_stream, const std::string extra_text
= "") const; | |
| | | | |
| inline void fill(const oT& x); | | inline void fill(const oT& x); | |
| | | | |
| inline void reset(); | | inline void reset(); | |
| inline void reset_objects(); | | inline void reset_objects(); | |
| | | | |
| arma_inline bool is_empty() const; | | arma_inline bool is_empty() const; | |
| | | | |
| arma_inline arma_warn_unused bool in_range(const uword i) const; | | arma_inline arma_warn_unused bool in_range(const uword i) const; | |
| arma_inline arma_warn_unused bool in_range(const span& x) const; | | arma_inline arma_warn_unused bool in_range(const span& x) const; | |
| | | | |
| arma_inline arma_warn_unused bool in_range(const uword in_row, const uw
ord in_col) const; | | arma_inline arma_warn_unused bool in_range(const uword in_row, const uw
ord in_col) const; | |
| arma_inline arma_warn_unused bool in_range(const span& row_span, const uw
ord in_col) const; | | arma_inline arma_warn_unused bool in_range(const span& row_span, const uw
ord in_col) const; | |
| arma_inline arma_warn_unused bool in_range(const uword in_row, const sp
an& col_span) const; | | arma_inline arma_warn_unused bool in_range(const uword in_row, const sp
an& col_span) const; | |
| arma_inline arma_warn_unused bool in_range(const span& row_span, const sp
an& col_span) const; | | arma_inline arma_warn_unused bool in_range(const span& row_span, const sp
an& col_span) const; | |
| | | | |
|
| arma_inline arma_warn_unused bool in_range(const uword in_row, const uwor | | arma_inline arma_warn_unused bool in_range(const uword in_row, const uw | |
| d in_col, const SizeMat& s) const; | | ord in_col, const SizeMat& s) const; | |
| | | | |
| | | arma_inline arma_warn_unused bool in_range(const uword in_row, const uw | |
| | | ord in_col, const uword in_slice) const; | |
| | | arma_inline arma_warn_unused bool in_range(const span& row_span, const sp | |
| | | an& col_span, const span& slice_span) const; | |
| | | | |
| | | arma_inline arma_warn_unused bool in_range(const uword in_row, const uw | |
| | | ord in_col, const uword in_slice, const SizeCube& s) const; | |
| | | | |
| inline bool save(const std::string name, const file_type type = arma_bi
nary, const bool print_status = true) const; | | inline bool save(const std::string name, const file_type type = arma_bi
nary, const bool print_status = true) const; | |
| inline bool save( std::ostream& os, const file_type type = arma_bi
nary, const bool print_status = true) const; | | inline bool save( std::ostream& os, const file_type type = arma_bi
nary, const bool print_status = true) const; | |
| | | | |
| inline bool load(const std::string name, const file_type type = auto_de
tect, const bool print_status = true); | | inline bool load(const std::string name, const file_type type = auto_de
tect, const bool print_status = true); | |
| inline bool load( std::istream& is, const file_type type = auto_de
tect, const bool print_status = true); | | inline bool load( std::istream& is, const file_type type = auto_de
tect, const bool print_status = true); | |
| | | | |
| inline bool quiet_save(const std::string name, const file_type type = a
rma_binary) const; | | inline bool quiet_save(const std::string name, const file_type type = a
rma_binary) const; | |
| inline bool quiet_save( std::ostream& os, const file_type type = a
rma_binary) const; | | inline bool quiet_save( std::ostream& os, const file_type type = a
rma_binary) const; | |
| | | | |
| | | | |
| skipping to change at line 197 | | skipping to change at line 232 | |
|---|
| inline const_iterator cend() const; | | inline const_iterator cend() const; | |
| | | | |
| inline void clear(); | | inline void clear(); | |
| inline bool empty() const; | | inline bool empty() const; | |
| inline uword size() const; | | inline uword size() const; | |
| | | | |
| private: | | private: | |
| | | | |
| inline void init(const field<oT>& x); | | inline void init(const field<oT>& x); | |
| inline void init(const uword n_rows_in, const uword n_cols_in); | | inline void init(const uword n_rows_in, const uword n_cols_in); | |
|
| | | inline void init(const uword n_rows_in, const uword n_cols_in, const uwor
d n_slices_in); | |
| | | | |
| inline void delete_objects(); | | inline void delete_objects(); | |
| inline void create_objects(); | | inline void create_objects(); | |
| | | | |
| friend class field_aux; | | friend class field_aux; | |
| friend class subview_field<oT>; | | friend class subview_field<oT>; | |
| | | | |
| public: | | public: | |
| | | | |
| #ifdef ARMA_EXTRA_FIELD_PROTO | | #ifdef ARMA_EXTRA_FIELD_PROTO | |
| | | | |
| End of changes. 14 change blocks. |
|---|
| 8 lines changed or deleted | | 65 lines changed or added | |
|---|
|
| field_meat.hpp | | field_meat.hpp | |
|---|
|
| // Copyright (C) 2008-2013 Conrad Sanderson | | // Copyright (C) 2008-2014 Conrad Sanderson | |
| // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2014 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2009-2010 Ian Cullinan | | // Copyright (C) 2009-2010 Ian Cullinan | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup field | | //! \addtogroup field | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| | | | |
| skipping to change at line 37 | | skipping to change at line 37 | |
|---|
| // try to expose buggy user code that accesses deleted objects | | // try to expose buggy user code that accesses deleted objects | |
| mem = 0; | | mem = 0; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field<oT>::field() | | field<oT>::field() | |
| : n_rows(0) | | : n_rows(0) | |
| , n_cols(0) | | , n_cols(0) | |
|
| | | , n_slices(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , mem(0) | | , mem(0) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| } | | } | |
| | | | |
| //! construct a field from a given field | | //! construct a field from a given field | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field<oT>::field(const field& x) | | field<oT>::field(const field& x) | |
| : n_rows(0) | | : n_rows(0) | |
| , n_cols(0) | | , n_cols(0) | |
|
| | | , n_slices(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , mem(0) | | , mem(0) | |
| { | | { | |
| arma_extra_debug_sigprint(arma_boost::format("this = %x x = %x") % this
% &x); | | arma_extra_debug_sigprint(arma_boost::format("this = %x x = %x") % this
% &x); | |
| | | | |
| init(x); | | init(x); | |
| } | | } | |
| | | | |
| //! construct a field from a given field | | //! construct a field from a given field | |
| template<typename oT> | | template<typename oT> | |
| | | | |
| skipping to change at line 75 | | skipping to change at line 77 | |
|---|
| init(x); | | init(x); | |
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| //! construct a field from subview_field (e.g. construct a field from a del
ayed subfield operation) | | //! construct a field from subview_field (e.g. construct a field from a del
ayed subfield operation) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field<oT>::field(const subview_field<oT>& X) | | field<oT>::field(const subview_field<oT>& X) | |
| : n_rows(0) | | : n_rows(0) | |
| , n_cols(0) | | , n_cols(0) | |
|
| | | , n_slices(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , mem(0) | | , mem(0) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| | | | |
| this->operator=(X); | | this->operator=(X); | |
| } | | } | |
| | | | |
| //! construct a field from subview_field (e.g. construct a field from a del
ayed subfield operation) | | //! construct a field from subview_field (e.g. construct a field from a del
ayed subfield operation) | |
| template<typename oT> | | template<typename oT> | |
| | | | |
| skipping to change at line 102 | | skipping to change at line 105 | |
|---|
| return *this; | | return *this; | |
| } | | } | |
| | | | |
| //! construct the field with the specified number of elements, | | //! construct the field with the specified number of elements, | |
| //! assuming a column-major layout | | //! assuming a column-major layout | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field<oT>::field(const uword n_elem_in) | | field<oT>::field(const uword n_elem_in) | |
| : n_rows(0) | | : n_rows(0) | |
| , n_cols(0) | | , n_cols(0) | |
|
| | | , n_slices(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , mem(0) | | , mem(0) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| | | | |
| init(n_elem_in, 1); | | init(n_elem_in, 1); | |
| } | | } | |
| | | | |
| //! construct the field with the specified dimensions | | //! construct the field with the specified dimensions | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field<oT>::field(const uword n_rows_in, const uword n_cols_in) | | field<oT>::field(const uword n_rows_in, const uword n_cols_in) | |
| : n_rows(0) | | : n_rows(0) | |
| , n_cols(0) | | , n_cols(0) | |
|
| | | , n_slices(0) | |
| , n_elem(0) | | , n_elem(0) | |
| , mem(0) | | , mem(0) | |
| { | | { | |
| arma_extra_debug_sigprint_this(this); | | arma_extra_debug_sigprint_this(this); | |
| | | | |
| init(n_rows_in, n_cols_in); | | init(n_rows_in, n_cols_in); | |
| } | | } | |
| | | | |
|
| | | //! construct the field with the specified dimensions | |
| | | template<typename oT> | |
| | | inline | |
| | | field<oT>::field(const uword n_rows_in, const uword n_cols_in, const uword | |
| | | n_slices_in) | |
| | | : n_rows(0) | |
| | | , n_cols(0) | |
| | | , n_slices(0) | |
| | | , n_elem(0) | |
| | | , mem(0) | |
| | | { | |
| | | arma_extra_debug_sigprint_this(this); | |
| | | | |
| | | init(n_rows_in, n_cols_in, n_slices_in); | |
| | | } | |
| | | | |
| //! change the field to have the specified number of elements, | | //! change the field to have the specified number of elements, | |
| //! assuming a column-major layout (data is not preserved) | | //! assuming a column-major layout (data is not preserved) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| field<oT>::set_size(const uword n_elem_in) | | field<oT>::set_size(const uword n_elem_in) | |
| { | | { | |
| arma_extra_debug_sigprint(arma_boost::format("n_elem_in = %d") % n_elem_i
n); | | arma_extra_debug_sigprint(arma_boost::format("n_elem_in = %d") % n_elem_i
n); | |
| | | | |
| init(n_elem_in, 1); | | init(n_elem_in, 1); | |
| | | | |
| skipping to change at line 149 | | skipping to change at line 169 | |
|---|
| void | | void | |
| field<oT>::set_size(const uword n_rows_in, const uword n_cols_in) | | field<oT>::set_size(const uword n_rows_in, const uword n_cols_in) | |
| { | | { | |
| arma_extra_debug_sigprint(arma_boost::format("n_rows_in = %d, n_cols_in =
%d") % n_rows_in % n_cols_in); | | arma_extra_debug_sigprint(arma_boost::format("n_rows_in = %d, n_cols_in =
%d") % n_rows_in % n_cols_in); | |
| | | | |
| init(n_rows_in, n_cols_in); | | init(n_rows_in, n_cols_in); | |
| } | | } | |
| | | | |
| //! change the field to have the specified dimensions (data is not preserve
d) | | //! change the field to have the specified dimensions (data is not preserve
d) | |
| template<typename oT> | | template<typename oT> | |
|
| | | inline | |
| | | void | |
| | | field<oT>::set_size(const uword n_rows_in, const uword n_cols_in, const uwo | |
| | | rd n_slices_in) | |
| | | { | |
| | | arma_extra_debug_sigprint(arma_boost::format("n_rows_in = %d, n_cols_in = | |
| | | %d, n_slices_in = %d") % n_rows_in % n_cols_in % n_slices_in); | |
| | | | |
| | | init(n_rows_in, n_cols_in, n_slices_in); | |
| | | } | |
| | | | |
| | | //! change the field to have the specified dimensions (data is not preserve | |
| | | d) | |
| | | template<typename oT> | |
| template<typename oT2> | | template<typename oT2> | |
| inline | | inline | |
| void | | void | |
| field<oT>::copy_size(const field<oT2>& x) | | field<oT>::copy_size(const field<oT2>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| init(x.n_rows, x.n_cols); | | init(x.n_rows, x.n_cols, x.n_slices); | |
| } | | } | |
| | | | |
| //! linear element accessor (treats the field as a vector); no bounds check | | //! linear element accessor (treats the field as a vector); no bounds check | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
| field<oT>::operator[] (const uword i) | | field<oT>::operator[] (const uword i) | |
| { | | { | |
| return (*mem[i]); | | return (*mem[i]); | |
| } | | } | |
| | | | |
| skipping to change at line 235 | | skipping to change at line 266 | |
|---|
| //! element accessor; bounds checking not done when ARMA_NO_DEBUG is define
d | | //! element accessor; bounds checking not done when ARMA_NO_DEBUG is define
d | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| const oT& | | const oT& | |
| field<oT>::operator() (const uword in_row, const uword in_col) const | | field<oT>::operator() (const uword in_row, const uword in_col) const | |
| { | | { | |
| arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "field::ope
rator(): index out of bounds"); | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "field::ope
rator(): index out of bounds"); | |
| return (*mem[in_row + in_col*n_rows]); | | return (*mem[in_row + in_col*n_rows]); | |
| } | | } | |
| | | | |
|
| | | //! element accessor; bounds checking not done when ARMA_NO_DEBUG is define | |
| | | d | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | oT& | |
| | | field<oT>::operator() (const uword in_row, const uword in_col, const uword | |
| | | in_slice) | |
| | | { | |
| | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols) || (in_slice | |
| | | >= n_slices)), "field::operator(): index out of bounds"); | |
| | | return (*mem[in_row + in_col*n_rows + in_slice*(n_rows*n_cols)]); | |
| | | } | |
| | | | |
| | | //! element accessor; bounds checking not done when ARMA_NO_DEBUG is define | |
| | | d | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | const oT& | |
| | | field<oT>::operator() (const uword in_row, const uword in_col, const uword | |
| | | in_slice) const | |
| | | { | |
| | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols) || (in_slice | |
| | | >= n_slices)), "field::operator(): index out of bounds"); | |
| | | return (*mem[in_row + in_col*n_rows + in_slice*(n_rows*n_cols)]); | |
| | | } | |
| | | | |
| //! element accessor; no bounds check | | //! element accessor; no bounds check | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
| field<oT>::at(const uword in_row, const uword in_col) | | field<oT>::at(const uword in_row, const uword in_col) | |
| { | | { | |
| return (*mem[in_row + in_col*n_rows]); | | return (*mem[in_row + in_col*n_rows]); | |
| } | | } | |
| | | | |
| //! element accessor; no bounds check | | //! element accessor; no bounds check | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| const oT& | | const oT& | |
| field<oT>::at(const uword in_row, const uword in_col) const | | field<oT>::at(const uword in_row, const uword in_col) const | |
| { | | { | |
| return (*mem[in_row + in_col*n_rows]); | | return (*mem[in_row + in_col*n_rows]); | |
| } | | } | |
| | | | |
|
| | | //! element accessor; no bounds check | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | oT& | |
| | | field<oT>::at(const uword in_row, const uword in_col, const uword in_slice) | |
| | | { | |
| | | return (*mem[in_row + in_col*n_rows + in_slice*(n_rows*n_cols)]); | |
| | | } | |
| | | | |
| | | //! element accessor; no bounds check | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | const oT& | |
| | | field<oT>::at(const uword in_row, const uword in_col, const uword in_slice) | |
| | | const | |
| | | { | |
| | | return (*mem[in_row + in_col*n_rows + in_slice*(n_rows*n_cols)]); | |
| | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| field_injector< field<oT> > | | field_injector< field<oT> > | |
| field<oT>::operator<<(const oT& val) | | field<oT>::operator<<(const oT& val) | |
| { | | { | |
| return field_injector< field<oT> >(*this, val); | | return field_injector< field<oT> >(*this, val); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| | | | |
| skipping to change at line 277 | | skipping to change at line 346 | |
|---|
| } | | } | |
| | | | |
| //! creation of subview_field (row of a field) | | //! creation of subview_field (row of a field) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::row(const uword row_num) | | field<oT>::row(const uword row_num) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::row(): field must be 2D" ); | |
| | | | |
| arma_debug_check( (row_num >= n_rows), "field::row(): row out of bounds"
); | | arma_debug_check( (row_num >= n_rows), "field::row(): row out of bounds"
); | |
| | | | |
| return subview_field<oT>(*this, row_num, 0, 1, n_cols); | | return subview_field<oT>(*this, row_num, 0, 1, n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (row of a field) | | //! creation of subview_field (row of a field) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::row(const uword row_num) const | | field<oT>::row(const uword row_num) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::row(): field must be 2D" ); | |
| | | | |
| arma_debug_check( (row_num >= n_rows), "field::row(): row out of bounds"
); | | arma_debug_check( (row_num >= n_rows), "field::row(): row out of bounds"
); | |
| | | | |
| return subview_field<oT>(*this, row_num, 0, 1, n_cols); | | return subview_field<oT>(*this, row_num, 0, 1, n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (column of a field) | | //! creation of subview_field (column of a field) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::col(const uword col_num) | | field<oT>::col(const uword col_num) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::col(): field must be 2D" ); | |
| | | | |
| arma_debug_check( (col_num >= n_cols), "field::col(): out of bounds"); | | arma_debug_check( (col_num >= n_cols), "field::col(): out of bounds"); | |
| | | | |
| return subview_field<oT>(*this, 0, col_num, n_rows, 1); | | return subview_field<oT>(*this, 0, col_num, n_rows, 1); | |
| } | | } | |
| | | | |
| //! creation of subview_field (column of a field) | | //! creation of subview_field (column of a field) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::col(const uword col_num) const | | field<oT>::col(const uword col_num) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::col(): field must be 2D" ); | |
| | | | |
| arma_debug_check( (col_num >= n_cols), "field::col(): out of bounds"); | | arma_debug_check( (col_num >= n_cols), "field::col(): out of bounds"); | |
| | | | |
| return subview_field<oT>(*this, 0, col_num, n_rows, 1); | | return subview_field<oT>(*this, 0, col_num, n_rows, 1); | |
| } | | } | |
| | | | |
|
| | | //! creation of subview_field (slice of a field) | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::slice(const uword slice_num) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check( (slice_num >= n_slices), "field::slice(): out of bounds | |
| | | "); | |
| | | | |
| | | return subview_field<oT>(*this, 0, 0, slice_num, n_rows, n_cols, 1); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (slice of a field) | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::slice(const uword slice_num) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check( (slice_num >= n_slices), "field::slice(): out of bounds | |
| | | "); | |
| | | | |
| | | return subview_field<oT>(*this, 0, 0, slice_num, n_rows, n_cols, 1); | |
| | | } | |
| | | | |
| //! creation of subview_field (subfield comprised of specified rows) | | //! creation of subview_field (subfield comprised of specified rows) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::rows(const uword in_row1, const uword in_row2) | | field<oT>::rows(const uword in_row1, const uword in_row2) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::rows(): field must be 2D" ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ( (in_row1 > in_row2) || (in_row2 >= n_rows) ), | | ( (in_row1 > in_row2) || (in_row2 >= n_rows) ), | |
| "field::rows(): indicies out of bounds or incorrectly used" | | "field::rows(): indicies out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_rows = in_row2 - in_row1 + 1; | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| | | | |
| return subview_field<oT>(*this, in_row1, 0, sub_n_rows, n_cols); | | return subview_field<oT>(*this, in_row1, 0, sub_n_rows, n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield comprised of specified rows) | | //! creation of subview_field (subfield comprised of specified rows) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::rows(const uword in_row1, const uword in_row2) const | | field<oT>::rows(const uword in_row1, const uword in_row2) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::rows(): field must be 2D" ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ( (in_row1 > in_row2) || (in_row2 >= n_rows) ), | | ( (in_row1 > in_row2) || (in_row2 >= n_rows) ), | |
| "field::rows(): indicies out of bounds or incorrectly used" | | "field::rows(): indicies out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_rows = in_row2 - in_row1 + 1; | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| | | | |
| return subview_field<oT>(*this, in_row1, 0, sub_n_rows, n_cols); | | return subview_field<oT>(*this, in_row1, 0, sub_n_rows, n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield comprised of specified columns) | | //! creation of subview_field (subfield comprised of specified columns) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::cols(const uword in_col1, const uword in_col2) | | field<oT>::cols(const uword in_col1, const uword in_col2) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::cols(): field must be 2D" ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ( (in_col1 > in_col2) || (in_col2 >= n_cols) ), | | ( (in_col1 > in_col2) || (in_col2 >= n_cols) ), | |
| "field::cols(): indicies out of bounds or incorrectly used" | | "field::cols(): indicies out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_cols = in_col2 - in_col1 + 1; | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | | |
| return subview_field<oT>(*this, 0, in_col1, n_rows, sub_n_cols); | | return subview_field<oT>(*this, 0, in_col1, n_rows, sub_n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield comprised of specified columns) | | //! creation of subview_field (subfield comprised of specified columns) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::cols(const uword in_col1, const uword in_col2) const | | field<oT>::cols(const uword in_col1, const uword in_col2) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::cols(): field must be 2D" ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ( (in_col1 > in_col2) || (in_col2 >= n_cols) ), | | ( (in_col1 > in_col2) || (in_col2 >= n_cols) ), | |
| "field::cols(): indicies out of bounds or incorrectly used" | | "field::cols(): indicies out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_cols = in_col2 - in_col1 + 1; | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | | |
| return subview_field<oT>(*this, 0, in_col1, n_rows, sub_n_cols); | | return subview_field<oT>(*this, 0, in_col1, n_rows, sub_n_cols); | |
| } | | } | |
| | | | |
|
| | | //! creation of subview_field (subfield comprised of specified slices) | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::slices(const uword in_slice1, const uword in_slice2) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (in_slice1 > in_slice2) || (in_slice2 >= n_slices) ), | |
| | | "field::slices(): indicies out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | const uword sub_n_slices = in_slice2 - in_slice1 + 1; | |
| | | | |
| | | return subview_field<oT>(*this, 0, 0, in_slice1, n_rows, n_cols, sub_n_sl | |
| | | ices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield comprised of specified slices) | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::slices(const uword in_slice1, const uword in_slice2) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( (in_slice1 > in_slice2) || (in_slice2 >= n_slices) ), | |
| | | "field::slices(): indicies out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | const uword sub_n_slices = in_slice2 - in_slice1 + 1; | |
| | | | |
| | | return subview_field<oT>(*this, 0, 0, in_slice1, n_rows, n_cols, sub_n_sl | |
| | | ices); | |
| | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i
n_row2, const uword in_col2) | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i
n_row2, const uword in_col2) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| (in_row1 > in_row2) || (in_col1 > in_col2) || (in_row2 >= n_rows) || (i
n_col2 >= n_cols), | | (in_row1 > in_row2) || (in_col1 > in_col2) || (in_row2 >= n_rows) || (i
n_col2 >= n_cols), | |
| "field::subfield(): indices out of bounds or incorrectly used" | | "field::subfield(): indices out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_rows = in_row2 - in_row1 + 1; | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| const uword sub_n_cols = in_col2 - in_col1 + 1; | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | | |
| return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | | return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i
n_row2, const uword in_col2) const | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i
n_row2, const uword in_col2) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| (in_row1 > in_row2) || (in_col1 > in_col2) || (in_row2 >= n_rows) || (i
n_col2 >= n_cols), | | (in_row1 > in_row2) || (in_col1 > in_col2) || (in_row2 >= n_rows) || (i
n_col2 >= n_cols), | |
| "field::subfield(): indices out of bounds or incorrectly used" | | "field::subfield(): indices out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| const uword sub_n_rows = in_row2 - in_row1 + 1; | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| const uword sub_n_cols = in_col2 - in_col1 + 1; | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | | |
| return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | | return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
|
| | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i | |
| | | n_slice1, const uword in_row2, const uword in_col2, const uword in_slice2) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | (in_row1 > in_row2) || (in_col1 > in_col2) || (in_slice1 > in_slice2) | | |
| | | | (in_row2 >= n_rows) || (in_col2 >= n_cols) || (in_slice2 >= n_slices), | |
| | | "field::subfield(): indices out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | const uword sub_n_slices = in_slice2 - in_slice1 + 1; | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, sub_n_rows, | |
| | | sub_n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i | |
| | | n_slice1, const uword in_row2, const uword in_col2, const uword in_slice2) | |
| | | const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | (in_row1 > in_row2) || (in_col1 > in_col2) || (in_slice1 > in_slice2) | | |
| | | | (in_row2 >= n_rows) || (in_col2 >= n_cols) || (in_slice2 >= n_slices), | |
| | | "field::subfield(): indices out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | const uword sub_n_rows = in_row2 - in_row1 + 1; | |
| | | const uword sub_n_cols = in_col2 - in_col1 + 1; | |
| | | const uword sub_n_slices = in_slice2 - in_slice1 + 1; | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, sub_n_rows, | |
| | | sub_n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| field<oT>::subfield(const uword in_row1, const uword in_col1, const SizeMat
& s) | | field<oT>::subfield(const uword in_row1, const uword in_col1, const SizeMat
& s) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| const uword l_n_rows = n_rows; | | const uword l_n_rows = n_rows; | |
| const uword l_n_cols = n_cols; | | const uword l_n_cols = n_cols; | |
| | | | |
| const uword s_n_rows = s.n_rows; | | const uword s_n_rows = s.n_rows; | |
| const uword s_n_cols = s.n_cols; | | const uword s_n_cols = s.n_cols; | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || ((in_row1 + s_n_rows
) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_cols)), | | ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || ((in_row1 + s_n_rows
) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_cols)), | |
| "field::subfield(): indices or size out of bounds" | | "field::subfield(): indices or size out of bounds" | |
| | | | |
| skipping to change at line 468 | | skipping to change at line 665 | |
|---|
| } | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::subfield(const uword in_row1, const uword in_col1, const SizeMat
& s) const | | field<oT>::subfield(const uword in_row1, const uword in_col1, const SizeMat
& s) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| const uword l_n_rows = n_rows; | | const uword l_n_rows = n_rows; | |
| const uword l_n_cols = n_cols; | | const uword l_n_cols = n_cols; | |
| | | | |
| const uword s_n_rows = s.n_rows; | | const uword s_n_rows = s.n_rows; | |
| const uword s_n_cols = s.n_cols; | | const uword s_n_cols = s.n_cols; | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || ((in_row1 + s_n_rows
) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_cols)), | | ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || ((in_row1 + s_n_rows
) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_cols)), | |
| "field::subfield(): indices or size out of bounds" | | "field::subfield(): indices or size out of bounds" | |
| ); | | ); | |
| | | | |
| return subview_field<oT>(*this, in_row1, in_col1, s_n_rows, s_n_cols); | | return subview_field<oT>(*this, in_row1, in_col1, s_n_rows, s_n_cols); | |
| } | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
|
| | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i | |
| | | n_slice1, const SizeCube& s) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const uword l_n_rows = n_rows; | |
| | | const uword l_n_cols = n_cols; | |
| | | const uword l_n_slices = n_slices; | |
| | | | |
| | | const uword s_n_rows = s.n_rows; | |
| | | const uword s_n_cols = s.n_cols; | |
| | | const uword sub_n_slices = s.n_slices; | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || (in_slice1 >= l_n_sl | |
| | | ices) || ((in_row1 + s_n_rows) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_c | |
| | | ols) || ((in_slice1 + sub_n_slices) > l_n_slices)), | |
| | | "field::subfield(): indices or size out of bounds" | |
| | | ); | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, s_n_rows, s_ | |
| | | n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::subfield(const uword in_row1, const uword in_col1, const uword i | |
| | | n_slice1, const SizeCube& s) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const uword l_n_rows = n_rows; | |
| | | const uword l_n_cols = n_cols; | |
| | | const uword l_n_slices = n_slices; | |
| | | | |
| | | const uword s_n_rows = s.n_rows; | |
| | | const uword s_n_cols = s.n_cols; | |
| | | const uword sub_n_slices = s.n_slices; | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ((in_row1 >= l_n_rows) || (in_col1 >= l_n_cols) || (in_slice1 >= l_n_sl | |
| | | ices) || ((in_row1 + s_n_rows) > l_n_rows) || ((in_col1 + s_n_cols) > l_n_c | |
| | | ols) || ((in_slice1 + sub_n_slices) > l_n_slices)), | |
| | | "field::subfield(): indices or size out of bounds" | |
| | | ); | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, s_n_rows, s_ | |
| | | n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| field<oT>::subfield(const span& row_span, const span& col_span) | | field<oT>::subfield(const span& row_span, const span& col_span) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| const bool row_all = row_span.whole; | | const bool row_all = row_span.whole; | |
| const bool col_all = col_span.whole; | | const bool col_all = col_span.whole; | |
| | | | |
| const uword local_n_rows = n_rows; | | const uword local_n_rows = n_rows; | |
| const uword local_n_cols = n_cols; | | const uword local_n_cols = n_cols; | |
| | | | |
| const uword in_row1 = row_all ? 0 : row_span.a; | | const uword in_row1 = row_all ? 0 : row_span.a; | |
| const uword in_row2 = row_span.b; | | const uword in_row2 = row_span.b; | |
| const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | | const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | |
| | | | |
| | | | |
| skipping to change at line 525 | | skipping to change at line 776 | |
|---|
| } | | } | |
| | | | |
| //! creation of subview_field (subfield with arbitrary dimensions) | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::subfield(const span& row_span, const span& col_span) const | | field<oT>::subfield(const span& row_span, const span& col_span) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | arma_debug_check( (n_slices >= 2), "field::subfield(): field must be 2D" | |
| | | ); | |
| | | | |
| const bool row_all = row_span.whole; | | const bool row_all = row_span.whole; | |
| const bool col_all = col_span.whole; | | const bool col_all = col_span.whole; | |
| | | | |
| const uword local_n_rows = n_rows; | | const uword local_n_rows = n_rows; | |
| const uword local_n_cols = n_cols; | | const uword local_n_cols = n_cols; | |
| | | | |
| const uword in_row1 = row_all ? 0 : row_span.a; | | const uword in_row1 = row_all ? 0 : row_span.a; | |
| const uword in_row2 = row_span.b; | | const uword in_row2 = row_span.b; | |
| const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | | const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | |
| | | | |
| | | | |
| skipping to change at line 551 | | skipping to change at line 804 | |
|---|
| ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows))
) | | ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows))
) | |
| || | | || | |
| ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols))
) | | ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols))
) | |
| , | | , | |
| "field::subfield(): indices out of bounds or incorrectly used" | | "field::subfield(): indices out of bounds or incorrectly used" | |
| ); | | ); | |
| | | | |
| return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | | return subview_field<oT>(*this, in_row1, in_col1, sub_n_rows, sub_n_cols)
; | |
| } | | } | |
| | | | |
|
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::subfield(const span& row_span, const span& col_span, const span& | |
| | | slice_span) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const bool row_all = row_span.whole; | |
| | | const bool col_all = col_span.whole; | |
| | | const bool slice_all = slice_span.whole; | |
| | | | |
| | | const uword local_n_rows = n_rows; | |
| | | const uword local_n_cols = n_cols; | |
| | | const uword local_n_slices = n_slices; | |
| | | | |
| | | const uword in_row1 = row_all ? 0 : row_span.a; | |
| | | const uword in_row2 = row_span.b; | |
| | | const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | |
| | | | |
| | | const uword in_col1 = col_all ? 0 : col_span.a; | |
| | | const uword in_col2 = col_span.b; | |
| | | const uword sub_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1; | |
| | | | |
| | | const uword in_slice1 = slice_all ? 0 : slice_span.a; | |
| | | const uword in_slice2 = slice_span.b; | |
| | | const uword sub_n_slices = slice_all ? local_n_slices : in_slice2 - in_sl | |
| | | ice1 + 1; | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) | |
| | | ) | |
| | | || | |
| | | ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) | |
| | | ) | |
| | | || | |
| | | ( slice_all ? false : ((in_slice1 > in_slice2) || (in_slice2 >= local_n | |
| | | _slices)) ) | |
| | | , | |
| | | "field::subfield(): indices out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, sub_n_rows, | |
| | | sub_n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| | | //! creation of subview_field (subfield with arbitrary dimensions) | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::subfield(const span& row_span, const span& col_span, const span& | |
| | | slice_span) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const bool row_all = row_span.whole; | |
| | | const bool col_all = col_span.whole; | |
| | | const bool slice_all = slice_span.whole; | |
| | | | |
| | | const uword local_n_rows = n_rows; | |
| | | const uword local_n_cols = n_cols; | |
| | | const uword local_n_slices = n_slices; | |
| | | | |
| | | const uword in_row1 = row_all ? 0 : row_span.a; | |
| | | const uword in_row2 = row_span.b; | |
| | | const uword sub_n_rows = row_all ? local_n_rows : in_row2 - in_row1 + 1; | |
| | | | |
| | | const uword in_col1 = col_all ? 0 : col_span.a; | |
| | | const uword in_col2 = col_span.b; | |
| | | const uword sub_n_cols = col_all ? local_n_cols : in_col2 - in_col1 + 1; | |
| | | | |
| | | const uword in_slice1 = slice_all ? 0 : slice_span.a; | |
| | | const uword in_slice2 = slice_span.b; | |
| | | const uword sub_n_slices = slice_all ? local_n_slices : in_slice2 - in_sl | |
| | | ice1 + 1; | |
| | | | |
| | | arma_debug_check | |
| | | ( | |
| | | ( row_all ? false : ((in_row1 > in_row2) || (in_row2 >= local_n_rows)) | |
| | | ) | |
| | | || | |
| | | ( col_all ? false : ((in_col1 > in_col2) || (in_col2 >= local_n_cols)) | |
| | | ) | |
| | | || | |
| | | ( slice_all ? false : ((in_slice1 > in_slice2) || (in_slice2 >= local_n | |
| | | _slices)) ) | |
| | | , | |
| | | "field::subfield(): indices out of bounds or incorrectly used" | |
| | | ); | |
| | | | |
| | | return subview_field<oT>(*this, in_row1, in_col1, in_slice1, sub_n_rows, | |
| | | sub_n_cols, sub_n_slices); | |
| | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| subview_field<oT> | | subview_field<oT> | |
| field<oT>::operator()(const span& row_span, const span& col_span) | | field<oT>::operator()(const span& row_span, const span& col_span) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| return (*this).subfield(row_span, col_span); | | return (*this).subfield(row_span, col_span); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| const subview_field<oT> | | const subview_field<oT> | |
| field<oT>::operator()(const span& row_span, const span& col_span) const | | field<oT>::operator()(const span& row_span, const span& col_span) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| return (*this).subfield(row_span, col_span); | | return (*this).subfield(row_span, col_span); | |
| } | | } | |
| | | | |
|
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::operator()(const span& row_span, const span& col_span, const spa | |
| | | n& slice_span) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(row_span, col_span, slice_span); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::operator()(const span& row_span, const span& col_span, const spa | |
| | | n& slice_span) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(row_span, col_span, slice_span); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::operator()(const uword in_row1, const uword in_col1, const SizeM | |
| | | at& s) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(in_row1, in_col1, s); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::operator()(const uword in_row1, const uword in_col1, const SizeM | |
| | | at& s) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(in_row1, in_col1, s); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | inline | |
| | | subview_field<oT> | |
| | | field<oT>::operator()(const uword in_row1, const uword in_col1, const uword | |
| | | in_slice1, const SizeCube& s) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(in_row1, in_col1, in_slice1, s); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | inline | |
| | | const subview_field<oT> | |
| | | field<oT>::operator()(const uword in_row1, const uword in_col1, const uword | |
| | | in_slice1, const SizeCube& s) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | return (*this).subfield(in_row1, in_col1, in_slice1, s); | |
| | | } | |
| | | | |
| //! print contents of the field (to the cout stream), | | //! print contents of the field (to the cout stream), | |
| //! optionally preceding with a user specified line of text. | | //! optionally preceding with a user specified line of text. | |
| //! the field class preserves the stream's flags | | //! the field class preserves the stream's flags | |
| //! but the associated operator<< function for type oT | | //! but the associated operator<< function for type oT | |
| //! may still modify the stream's parameters. | | //! may still modify the stream's parameters. | |
| //! NOTE: this function assumes that type oT can be printed, | | //! NOTE: this function assumes that type oT can be printed, | |
| //! i.e. the function "std::ostream& operator<< (std::ostream&, const oT&)" | | //! i.e. the function "std::ostream& operator<< (std::ostream&, const oT&)" | |
| //! has been defined. | | //! has been defined. | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| | | | |
| skipping to change at line 651 | | skipping to change at line 1048 | |
|---|
| } | | } | |
| | | | |
| //! reset the field to an empty state (i.e. the field will have no objects) | | //! reset the field to an empty state (i.e. the field will have no objects) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| field<oT>::reset() | | field<oT>::reset() | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| init(0,0); | | init(0,0,0); | |
| } | | } | |
| | | | |
| //! reset each object | | //! reset each object | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| field<oT>::reset_objects() | | field<oT>::reset_objects() | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 798 | | skipping to change at line 1195 | |
|---|
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| else | | else | |
| { | | { | |
| return true; | | return true; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
|
| | | arma_inline | |
| | | arma_warn_unused | |
| | | bool | |
| | | field<oT>::in_range(const uword in_row, const uword in_col, const uword in_ | |
| | | slice) const | |
| | | { | |
| | | return ( (in_row < n_rows) && (in_col < n_cols) && (in_slice < n_slices) | |
| | | ); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | arma_warn_unused | |
| | | bool | |
| | | field<oT>::in_range(const span& row_span, const span& col_span, const span& | |
| | | slice_span) const | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const uword in_row1 = row_span.a; | |
| | | const uword in_row2 = row_span.b; | |
| | | | |
| | | const uword in_col1 = col_span.a; | |
| | | const uword in_col2 = col_span.b; | |
| | | | |
| | | const uword in_slice1 = slice_span.a; | |
| | | const uword in_slice2 = slice_span.b; | |
| | | | |
| | | const bool rows_ok = row_span.whole ? true : ( (in_row1 <= in_row2 | |
| | | ) && (in_row2 < n_rows ) ); | |
| | | const bool cols_ok = col_span.whole ? true : ( (in_col1 <= in_col2 | |
| | | ) && (in_col2 < n_cols ) ); | |
| | | const bool slices_ok = slice_span.whole ? true : ( (in_slice1 <= in_slice | |
| | | 2) && (in_slice2 < n_slices) ); | |
| | | | |
| | | return ( (rows_ok == true) && (cols_ok == true) && (slices_ok == true) ); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | arma_warn_unused | |
| | | bool | |
| | | field<oT>::in_range(const uword in_row, const uword in_col, const uword in_ | |
| | | slice, const SizeCube& s) const | |
| | | { | |
| | | const uword l_n_rows = n_rows; | |
| | | const uword l_n_cols = n_cols; | |
| | | const uword l_n_slices = n_slices; | |
| | | | |
| | | if( (in_row >= l_n_rows) || (in_col >= l_n_cols) || (in_slice >= l_n_slic | |
| | | es) || ((in_row + s.n_rows) > l_n_rows) || ((in_col + s.n_cols) > l_n_cols) | |
| | | || ((in_slice + s.n_slices) > l_n_slices) ) | |
| | | { | |
| | | return false; | |
| | | } | |
| | | else | |
| | | { | |
| | | return true; | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| inline | | inline | |
| bool | | bool | |
| field<oT>::save(const std::string name, const file_type type, const bool pr
int_status) const | | field<oT>::save(const std::string name, const file_type type, const bool pr
int_status) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| std::string err_msg; | | std::string err_msg; | |
| const bool save_okay = field_aux::save(*this, name, type, err_msg); | | const bool save_okay = field_aux::save(*this, name, type, err_msg); | |
| | | | |
| if( (print_status == true) && (save_okay == false) ) | | if( (print_status == true) && (save_okay == false) ) | |
| | | | |
| skipping to change at line 957 | | skipping to change at line 1407 | |
|---|
| //! construct a field from a given field | | //! construct a field from a given field | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| field<oT>::init(const field<oT>& x) | | field<oT>::init(const field<oT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(this != &x) | | if(this != &x) | |
| { | | { | |
|
| const uword x_n_rows = x.n_rows; | | const uword x_n_rows = x.n_rows; | |
| const uword x_n_cols = x.n_cols; | | const uword x_n_cols = x.n_cols; | |
| | | const uword x_n_slices = x.n_slices; | |
| | | | |
|
| init(x_n_rows, x_n_cols); | | init(x_n_rows, x_n_cols, x_n_slices); | |
| | | | |
| field& t = *this; | | field& t = *this; | |
| | | | |
|
| for(uword ucol=0; ucol < x_n_cols; ++ucol) | | if(x_n_slices == 1) | |
| for(uword urow=0; urow < x_n_rows; ++urow) | | | |
| { | | { | |
|
| t.at(urow,ucol) = x.at(urow,ucol); | | for(uword ucol=0; ucol < x_n_cols; ++ucol) | |
| | | for(uword urow=0; urow < x_n_rows; ++urow) | |
| | | { | |
| | | t.at(urow,ucol) = x.at(urow,ucol); | |
| | | } | |
| | | } | |
| | | else | |
| | | { | |
| | | for(uword uslice=0; uslice < x_n_slices; ++uslice) | |
| | | for(uword ucol=0; ucol < x_n_cols; ++ucol ) | |
| | | for(uword urow=0; urow < x_n_rows; ++urow ) | |
| | | { | |
| | | t.at(urow,ucol,uslice) = x.at(urow,ucol,uslice); | |
| | | } | |
| } | | } | |
| } | | } | |
|
| | | } | |
| | | | |
|
| | | template<typename oT> | |
| | | inline | |
| | | void | |
| | | field<oT>::init(const uword n_rows_in, const uword n_cols_in) | |
| | | { | |
| | | (*this).init(n_rows_in, n_cols_in, 1); | |
| } | | } | |
| | | | |
| //! internal field construction; if the requested size is small enough, mem
ory from the stack is used. otherwise memory is allocated via 'new' | | //! internal field construction; if the requested size is small enough, mem
ory from the stack is used. otherwise memory is allocated via 'new' | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
|
| field<oT>::init(const uword n_rows_in, const uword n_cols_in) | | field<oT>::init(const uword n_rows_in, const uword n_cols_in, const uword n
_slices_in) | |
| { | | { | |
|
| arma_extra_debug_sigprint( arma_boost::format("n_rows_in = %d, n_cols_in
= %d") % n_rows_in % n_cols_in ); | | arma_extra_debug_sigprint( arma_boost::format("n_rows_in = %d, n_cols_in
= %d, n_slices_in = %d") % n_rows_in % n_cols_in % n_slices_in ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ( | | ( | |
|
| ( (n_rows_in > ARMA_MAX_UHWORD) || (n_cols_in > ARMA_MAX_UHWORD) ) | | ( (n_rows_in > 0x0FFF) || (n_cols_in > 0x0FFF) || (n_slices_in > 0xFF | |
| ? ( (float(n_rows_in) * float(n_cols_in)) > float(ARMA_MAX_UWORD) ) | | ) ) | |
| | | ? ( (float(n_rows_in) * float(n_cols_in) * float(n_slices_in)) > fl | |
| | | oat(ARMA_MAX_UWORD) ) | |
| : false | | : false | |
| ), | | ), | |
| "field::init(): requested size is too large; suggest to enable ARMA_64B
IT_WORD" | | "field::init(): requested size is too large; suggest to enable ARMA_64B
IT_WORD" | |
| ); | | ); | |
| | | | |
|
| const uword n_elem_new = n_rows_in * n_cols_in; | | const uword n_elem_new = n_rows_in * n_cols_in * n_slices_in; | |
| | | | |
| if(n_elem == n_elem_new) | | if(n_elem == n_elem_new) | |
| { | | { | |
| // delete_objects(); | | // delete_objects(); | |
| // create_objects(); | | // create_objects(); | |
|
| access::rw(n_rows) = n_rows_in; | | access::rw(n_rows) = n_rows_in; | |
| access::rw(n_cols) = n_cols_in; | | access::rw(n_cols) = n_cols_in; | |
| | | access::rw(n_slices) = n_slices_in; | |
| } | | } | |
| else | | else | |
| { | | { | |
| delete_objects(); | | delete_objects(); | |
| | | | |
| if(n_elem > sizeof(mem_local)/sizeof(oT*) ) | | if(n_elem > sizeof(mem_local)/sizeof(oT*) ) | |
| { | | { | |
| delete [] mem; | | delete [] mem; | |
| } | | } | |
| | | | |
| | | | |
| skipping to change at line 1023 | | skipping to change at line 1494 | |
|---|
| else | | else | |
| { | | { | |
| mem = new(std::nothrow) oT* [n_elem_new]; | | mem = new(std::nothrow) oT* [n_elem_new]; | |
| arma_check_bad_alloc( (mem == 0), "field::init(): out of memory" ); | | arma_check_bad_alloc( (mem == 0), "field::init(): out of memory" ); | |
| } | | } | |
| | | | |
| access::rw(n_elem) = n_elem_new; | | access::rw(n_elem) = n_elem_new; | |
| | | | |
| if(n_elem_new == 0) | | if(n_elem_new == 0) | |
| { | | { | |
|
| access::rw(n_rows) = 0; | | access::rw(n_rows) = 0; | |
| access::rw(n_cols) = 0; | | access::rw(n_cols) = 0; | |
| | | access::rw(n_slices) = 0; | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| access::rw(n_rows) = n_rows_in; | | access::rw(n_rows) = n_rows_in; | |
| access::rw(n_cols) = n_cols_in; | | access::rw(n_cols) = n_cols_in; | |
| | | access::rw(n_slices) = n_slices_in; | |
| } | | } | |
| | | | |
| create_objects(); | | create_objects(); | |
| | | | |
| } | | } | |
| | | | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| | | | |
| End of changes. 46 change blocks. |
|---|
| 21 lines changed or deleted | | 559 lines changed or added | |
|---|
|
| sp_auxlib_meat.hpp | | sp_auxlib_meat.hpp | |
|---|
|
| // Copyright (C) 2013 Ryan Curtin | | // Copyright (C) 2013-2014 Ryan Curtin | |
| // Copyright (C) 2013 Conrad Sanderson | | // Copyright (C) 2013-2014 Conrad Sanderson | |
| | | // Copyright (C) 2013-2014 NICTA | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup sp_auxlib | | //! \addtogroup sp_auxlib | |
| //! @{ | | //! @{ | |
| | | | |
|
| | | inline | |
| | | sp_auxlib::form_type | |
| | | sp_auxlib::interpret_form_str(const char* form_str) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | // the order of the 3 if statements below is important | |
| | | if( form_str == NULL ) { return form_none; } | |
| | | if( form_str[0] == char(0) ) { return form_none; } | |
| | | if( form_str[1] == char(0) ) { return form_none; } | |
| | | | |
| | | const char c1 = form_str[0]; | |
| | | const char c2 = form_str[1]; | |
| | | | |
| | | if(c1 == 'l') | |
| | | { | |
| | | if(c2 == 'm') { return form_lm; } | |
| | | if(c2 == 'r') { return form_lr; } | |
| | | if(c2 == 'i') { return form_li; } | |
| | | } | |
| | | else | |
| | | if(c1 == 's') | |
| | | { | |
| | | if(c2 == 'm') { return form_sm; } | |
| | | if(c2 == 'r') { return form_sr; } | |
| | | if(c2 == 'i') { return form_si; } | |
| | | } | |
| | | | |
| | | return form_none; | |
| | | } | |
| | | | |
| //! immediate eigendecomposition of symmetric real sparse object | | //! immediate eigendecomposition of symmetric real sparse object | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
| bool | | bool | |
|
| sp_auxlib::eigs_sym(Col<eT>& eigval, Mat<eT>& eigvec, const SpBase<eT, T1>&
X, const uword n_eigvals) | | sp_auxlib::eigs_sym(Col<eT>& eigval, Mat<eT>& eigvec, const SpBase<eT, T1>&
X, const uword n_eigvals, const char* form_str) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_ARPACK) | | #if defined(ARMA_USE_ARPACK) | |
| { | | { | |
|
| | | const form_type form_val = sp_auxlib::interpret_form_str(form_str); | |
| | | | |
| | | arma_debug_check( (form_val != form_lm) && (form_val != form_sm), "eigs | |
| | | _sym(): unknown form specified" ); | |
| | | | |
| | | char which_sm[3] = "SM"; | |
| | | char which_lm[3] = "LM"; | |
| | | char* which = (form_val == form_sm) ? which_sm : which_lm; // sel | |
| | | ect which eigenvalues we want: smallest magnitude or largest magnitude | |
| | | | |
| // Make a sparse proxy object. | | // Make a sparse proxy object. | |
| SpProxy<T1> p(X.get_ref()); | | SpProxy<T1> p(X.get_ref()); | |
| | | | |
| // Make sure it's square. | | // Make sure it's square. | |
| arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_sym(): give
n sparse matrix is not square"); | | arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_sym(): give
n sparse matrix is not square"); | |
| | | | |
| // Make sure we aren't asking for every eigenvalue. | | // Make sure we aren't asking for every eigenvalue. | |
| arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_sym(): n_eigvals
must be less than the number of rows in the matrix"); | | arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_sym(): n_eigvals
must be less than the number of rows in the matrix"); | |
| | | | |
| // If the matrix is empty, the case is trivial. | | // If the matrix is empty, the case is trivial. | |
|
| if (p.get_n_cols() == 0) // We already know n_cols == n_rows. | | if(p.get_n_cols() == 0) // We already know n_cols == n_rows. | |
| { | | { | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| // Set up variables that get used for neupd(). | | // Set up variables that get used for neupd(). | |
| blas_int n, ncv, ldv, lworkl, info; | | blas_int n, ncv, ldv, lworkl, info; | |
| eT tol; | | eT tol; | |
| podarray<eT> resid, v, workd, workl; | | podarray<eT> resid, v, workd, workl; | |
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<eT> rwork; // Not used in this case. | | podarray<eT> rwork; // Not used in this case. | |
| | | | |
|
| run_aupd(n_eigvals, p, true /* sym, not gen */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, which, p, true /* sym, not gen */, n, tol, resid, n
cv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
|
| // The process has converged, and now we need to recover the actual eig
envectors using seupd(). | | // The process has converged, and now we need to recover the actual eig
envectors using seupd() | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
|
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| | | | |
| char howmny = 'A'; | | char howmny = 'A'; | |
|
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue proble | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni | | m. | |
| tude. | | | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| | | | |
| // seupd() will output directly into the eigval and eigvec objects. | | // seupd() will output directly into the eigval and eigvec objects. | |
| eigval.set_size(n_eigvals); | | eigval.set_size(n_eigvals); | |
| eigvec.set_size(n, n_eigvals); | | eigvec.set_size(n, n_eigvals); | |
| | | | |
| arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.
memptr(), &ldz, (eT*) NULL, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | | arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.
memptr(), &ldz, (eT*) NULL, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | |
| | | | |
| // Check for errors. | | // Check for errors. | |
| | | | |
| skipping to change at line 94 | | skipping to change at line 135 | |
|---|
| arma_stop("eigs_sym(): use of ARPACK needs to be enabled"); | | arma_stop("eigs_sym(): use of ARPACK needs to be enabled"); | |
| return false; | | return false; | |
| } | | } | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| //! immediate eigendecomposition of non-symmetric real sparse object | | //! immediate eigendecomposition of non-symmetric real sparse object | |
| template<typename T, typename T1> | | template<typename T, typename T1> | |
| inline | | inline | |
| bool | | bool | |
|
| sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >&
eigvec, const SpBase<T, T1>& X, const uword n_eigvals) | | sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >&
eigvec, const SpBase<T, T1>& X, const uword n_eigvals, const char* form_st
r) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_ARPACK) | | #if defined(ARMA_USE_ARPACK) | |
| { | | { | |
|
| | | const form_type form_val = sp_auxlib::interpret_form_str(form_str); | |
| | | | |
| | | arma_debug_check( (form_val == form_none), "eigs_gen(): unknown form sp | |
| | | ecified" ); | |
| | | | |
| | | char which_lm[3] = "LM"; | |
| | | char which_sm[3] = "SM"; | |
| | | char which_lr[3] = "LR"; | |
| | | char which_sr[3] = "SR"; | |
| | | char which_li[3] = "LI"; | |
| | | char which_si[3] = "SI"; | |
| | | | |
| | | char* which; | |
| | | | |
| | | switch(form_val) | |
| | | { | |
| | | case form_lm: which = which_lm; break; | |
| | | case form_sm: which = which_sm; break; | |
| | | case form_lr: which = which_lr; break; | |
| | | case form_sr: which = which_sr; break; | |
| | | case form_li: which = which_li; break; | |
| | | case form_si: which = which_si; break; | |
| | | | |
| | | default: which = which_lm; | |
| | | } | |
| | | | |
| // Make a sparse proxy object. | | // Make a sparse proxy object. | |
| SpProxy<T1> p(X.get_ref()); | | SpProxy<T1> p(X.get_ref()); | |
| | | | |
| // Make sure it's square. | | // Make sure it's square. | |
| arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_gen(): give
n sparse matrix is not square"); | | arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_gen(): give
n sparse matrix is not square"); | |
| | | | |
| // Make sure we aren't asking for every eigenvalue. | | // Make sure we aren't asking for every eigenvalue. | |
| arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_gen(): n_eigvals
must be less than the number of rows in the matrix"); | | arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_gen(): n_eigvals
must be less than the number of rows in the matrix"); | |
| | | | |
| // If the matrix is empty, the case is trivial. | | // If the matrix is empty, the case is trivial. | |
|
| if (p.get_n_cols() == 0) // We already know n_cols == n_rows. | | if(p.get_n_cols() == 0) // We already know n_cols == n_rows. | |
| { | | { | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| // Set up variables that get used for neupd(). | | // Set up variables that get used for neupd(). | |
| blas_int n, ncv, ldv, lworkl, info; | | blas_int n, ncv, ldv, lworkl, info; | |
| T tol; | | T tol; | |
| podarray<T> resid, v, workd, workl; | | podarray<T> resid, v, workd, workl; | |
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<T> rwork; // Not used in the real case. | | podarray<T> rwork; // Not used in the real case. | |
| | | | |
|
| run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, which, p, false /* gen, not sym */, n, tol, resid,
ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
| // The process has converged, and now we need to recover the actual eig
envectors using neupd(). | | // The process has converged, and now we need to recover the actual eig
envectors using neupd(). | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
|
| | | | |
| char howmny = 'A'; | | char howmny = 'A'; | |
|
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue proble | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni | | m. | |
| tude. | | | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| podarray<T> dr(nev + 1); // Real array of dimension NEV + 1. | | podarray<T> dr(nev + 1); // Real array of dimension NEV + 1. | |
| podarray<T> di(nev + 1); // Real array of dimension NEV + 1. | | podarray<T> di(nev + 1); // Real array of dimension NEV + 1. | |
| podarray<T> z(n * (nev + 1)); // Real N by NEV array if HOWMNY = 'A'. | | podarray<T> z(n * (nev + 1)); // Real N by NEV array if HOWMNY = 'A'. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| podarray<T> workev(3 * ncv); | | podarray<T> workev(3 * ncv); | |
| | | | |
| arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr()
, z.memptr(), &ldz, (T*) NULL, (T*) NULL, workev.memptr(), &bmat, &n, which
, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipnt
r.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info)
; | | arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr()
, z.memptr(), &ldz, (T*) NULL, (T*) NULL, workev.memptr(), &bmat, &n, which
, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipnt
r.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info)
; | |
| | | | |
| // Check for errors. | | // Check for errors. | |
| | | | |
| skipping to change at line 169 | | skipping to change at line 236 | |
|---|
| for (uword i = 0; i < n_eigvals; ++i) | | for (uword i = 0; i < n_eigvals; ++i) | |
| { | | { | |
| eigval[i] = std::complex<T>(dr[i], di[i]); | | eigval[i] = std::complex<T>(dr[i], di[i]); | |
| } | | } | |
| | | | |
| // Now recover the eigenvectors. | | // Now recover the eigenvectors. | |
| for (uword i = 0; i < n_eigvals; ++i) | | for (uword i = 0; i < n_eigvals; ++i) | |
| { | | { | |
| // ARPACK ?neupd lays things out kinda odd in memory; so does LAPACK | | // ARPACK ?neupd lays things out kinda odd in memory; so does LAPACK | |
| // ?geev (see auxlib::eig_gen()). | | // ?geev (see auxlib::eig_gen()). | |
|
| if ((i < n_eigvals - 1) && (eigval[i] == std::conj(eigval[i + 1]))) | | if((i < n_eigvals - 1) && (eigval[i] == std::conj(eigval[i + 1]))) | |
| { | | { | |
| for (uword j = 0; j < n; ++j) | | for (uword j = 0; j < n; ++j) | |
| { | | { | |
| eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1)
+ j]); | | eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1)
+ j]); | |
| eigvec.at(j, i + 1) = std::complex<T>(z[n * i + j], -z[n * (i + 1
) + j]); | | eigvec.at(j, i + 1) = std::complex<T>(z[n * i + j], -z[n * (i + 1
) + j]); | |
| } | | } | |
| ++i; // Skip the next one. | | ++i; // Skip the next one. | |
| } | | } | |
| else | | else | |
|
| if ((i == n_eigvals - 1) && (std::complex<T>(eigval[i]).imag() != 0.0
)) | | if((i == n_eigvals - 1) && (std::complex<T>(eigval[i]).imag() != 0.0)
) | |
| { | | { | |
| // We don't have the matched conjugate eigenvalue. | | // We don't have the matched conjugate eigenvalue. | |
| for (uword j = 0; j < n; ++j) | | for (uword j = 0; j < n; ++j) | |
| { | | { | |
| eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1) + j
]); | | eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1) + j
]); | |
| } | | } | |
| } | | } | |
| else | | else | |
| { | | { | |
| // The eigenvector is entirely real. | | // The eigenvector is entirely real. | |
| | | | |
| skipping to change at line 209 | | skipping to change at line 276 | |
|---|
| #else | | #else | |
| { | | { | |
| arma_ignore(eigval); | | arma_ignore(eigval); | |
| arma_ignore(eigvec); | | arma_ignore(eigvec); | |
| arma_ignore(X); | | arma_ignore(X); | |
| arma_ignore(n_eigvals); | | arma_ignore(n_eigvals); | |
| arma_stop("eigs_gen(): use of ARPACK needs to be enabled"); | | arma_stop("eigs_gen(): use of ARPACK needs to be enabled"); | |
| return false; | | return false; | |
| } | | } | |
| #endif | | #endif | |
|
| | | | |
| } | | } | |
| | | | |
| //! immediate eigendecomposition of non-symmetric complex sparse object | | //! immediate eigendecomposition of non-symmetric complex sparse object | |
| template<typename T, typename T1> | | template<typename T, typename T1> | |
| inline | | inline | |
| bool | | bool | |
|
| sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >&
eigvec, const SpBase< std::complex<T>, T1>& X, const uword n_eigvals) | | sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >&
eigvec, const SpBase< std::complex<T>, T1>& X, const uword n_eigvals, cons
t char* form_str) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_ARPACK) | | #if defined(ARMA_USE_ARPACK) | |
| { | | { | |
|
| | | const form_type form_val = sp_auxlib::interpret_form_str(form_str); | |
| | | | |
| | | arma_debug_check( (form_val == form_none), "eigs_gen(): unknown form sp | |
| | | ecified" ); | |
| | | | |
| | | char which_lm[3] = "LM"; | |
| | | char which_sm[3] = "SM"; | |
| | | char which_lr[3] = "LR"; | |
| | | char which_sr[3] = "SR"; | |
| | | char which_li[3] = "LI"; | |
| | | char which_si[3] = "SI"; | |
| | | | |
| | | char* which; | |
| | | | |
| | | switch(form_val) | |
| | | { | |
| | | case form_lm: which = which_lm; break; | |
| | | case form_sm: which = which_sm; break; | |
| | | case form_lr: which = which_lr; break; | |
| | | case form_sr: which = which_sr; break; | |
| | | case form_li: which = which_li; break; | |
| | | case form_si: which = which_si; break; | |
| | | | |
| | | default: which = which_lm; | |
| | | } | |
| | | | |
| // Make a sparse proxy object. | | // Make a sparse proxy object. | |
| SpProxy<T1> p(X.get_ref()); | | SpProxy<T1> p(X.get_ref()); | |
| | | | |
| // Make sure it's square. | | // Make sure it's square. | |
| arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_gen(): give
n sparse matrix is not square"); | | arma_debug_check( (p.get_n_rows() != p.get_n_cols()), "eigs_gen(): give
n sparse matrix is not square"); | |
| | | | |
| // Make sure we aren't asking for every eigenvalue. | | // Make sure we aren't asking for every eigenvalue. | |
| arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_gen(): n_eigvals
must be less than the number of rows in the matrix"); | | arma_debug_check( (n_eigvals >= p.get_n_rows()), "eigs_gen(): n_eigvals
must be less than the number of rows in the matrix"); | |
| | | | |
| // If the matrix is empty, the case is trivial. | | // If the matrix is empty, the case is trivial. | |
|
| if (p.get_n_cols() == 0) // We already know n_cols == n_rows. | | if(p.get_n_cols() == 0) // We already know n_cols == n_rows. | |
| { | | { | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| // Set up variables that get used for neupd(). | | // Set up variables that get used for neupd(). | |
| blas_int n, ncv, ldv, lworkl, info; | | blas_int n, ncv, ldv, lworkl, info; | |
| T tol; | | T tol; | |
| podarray< std::complex<T> > resid, v, workd, workl; | | podarray< std::complex<T> > resid, v, workd, workl; | |
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<T> rwork; | | podarray<T> rwork; | |
| | | | |
|
| run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, which, p, false /* gen, not sym */, n, tol, resid,
ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
| // The process has converged, and now we need to recover the actual eig
envectors using neupd(). | | // The process has converged, and now we need to recover the actual eig
envectors using neupd(). | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
|
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| | | | |
| char howmny = 'A'; | | char howmny = 'A'; | |
|
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue proble | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni | | m. | |
| tude. | | | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| podarray<std::complex<T> > d(nev + 1); // Real array of dimension NEV +
1. | | podarray<std::complex<T> > d(nev + 1); // Real array of dimension NEV +
1. | |
| podarray<std::complex<T> > z(n * nev); // Real N by NEV array if HOWMNY
= 'A'. | | podarray<std::complex<T> > z(n * nev); // Real N by NEV array if HOWMNY
= 'A'. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| podarray<std::complex<T> > workev(2 * ncv); | | podarray<std::complex<T> > workev(2 * ncv); | |
| | | | |
| // Prepare the outputs; neupd() will write directly to them. | | // Prepare the outputs; neupd() will write directly to them. | |
| eigval.set_size(n_eigvals); | | eigval.set_size(n_eigvals); | |
| eigvec.set_size(n, n_eigvals); | | eigvec.set_size(n, n_eigvals); | |
| std::complex<T> sigma; | | std::complex<T> sigma; | |
| | | | |
| skipping to change at line 294 | | skipping to change at line 386 | |
|---|
| #else | | #else | |
| { | | { | |
| arma_ignore(eigval); | | arma_ignore(eigval); | |
| arma_ignore(eigvec); | | arma_ignore(eigvec); | |
| arma_ignore(X); | | arma_ignore(X); | |
| arma_ignore(n_eigvals); | | arma_ignore(n_eigvals); | |
| arma_stop("eigs_gen(): use of ARPACK needs to be enabled"); | | arma_stop("eigs_gen(): use of ARPACK needs to be enabled"); | |
| return false; | | return false; | |
| } | | } | |
| #endif | | #endif | |
|
| | | | |
| } | | } | |
| | | | |
| template<typename eT, typename T, typename T1> | | template<typename eT, typename T, typename T1> | |
| inline | | inline | |
| void | | void | |
| sp_auxlib::run_aupd | | sp_auxlib::run_aupd | |
| ( | | ( | |
|
| const uword n_eigvals, const SpProxy<T1>& p, const bool sym, | | const uword n_eigvals, char* which, const SpProxy<T1>& p, const bool sym, | |
| blas_int& n, eT& tol, | | blas_int& n, eT& tol, | |
| podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv, | | podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv, | |
| podarray<blas_int>& iparam, podarray<blas_int>& ipntr, | | podarray<blas_int>& iparam, podarray<blas_int>& ipntr, | |
| podarray<T>& workd, podarray<T>& workl, blas_int& lworkl, podarray<eT>& r
work, | | podarray<T>& workd, podarray<T>& workl, blas_int& lworkl, podarray<eT>& r
work, | |
| blas_int& info | | blas_int& info | |
| ) | | ) | |
| { | | { | |
| #if defined(ARMA_USE_ARPACK) | | #if defined(ARMA_USE_ARPACK) | |
| { | | { | |
| // ARPACK provides a "reverse communication interface" which is an | | // ARPACK provides a "reverse communication interface" which is an | |
| // entertainingly archaic FORTRAN software engineering technique that | | // entertainingly archaic FORTRAN software engineering technique that | |
| // basically means that we call saupd()/naupd() and it tells us with so
me | | // basically means that we call saupd()/naupd() and it tells us with so
me | |
| // return code what we need to do next (usually a matrix-vector product
) and | | // return code what we need to do next (usually a matrix-vector product
) and | |
| // then call it again. So this results in some type of iterative proce
ss | | // then call it again. So this results in some type of iterative proce
ss | |
| // where we call saupd()/naupd() many times. | | // where we call saupd()/naupd() many times. | |
| blas_int ido = 0; // This must be 0 for the first call. | | blas_int ido = 0; // This must be 0 for the first call. | |
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue problem. | |
| n = p.get_n_rows(); // The size of the matrix. | | n = p.get_n_rows(); // The size of the matrix. | |
|
| char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | | | |
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| | | | |
| tol = std::numeric_limits<eT>::epsilon(); // Machine epsilon as toleran
ce. | | tol = std::numeric_limits<eT>::epsilon(); // Machine epsilon as toleran
ce. | |
| resid.set_size(n); | | resid.set_size(n); | |
| | | | |
| // "NCV must satisfy the two inequalities 2 <= NCV-NEV and NCV <= N". | | // "NCV must satisfy the two inequalities 2 <= NCV-NEV and NCV <= N". | |
| // "It is recommended that NCV >= 2 * NEV". | | // "It is recommended that NCV >= 2 * NEV". | |
| ncv = (2 * nev < n) ? 2 * nev : ((nev + 2 < n) ? nev + 2 : n); | | ncv = (2 * nev < n) ? 2 * nev : ((nev + 2 < n) ? nev + 2 : n); | |
| v.set_size(n * ncv); // Array N by NCV (output). | | v.set_size(n * ncv); // Array N by NCV (output). | |
| rwork.set_size(ncv); // Work array of size NCV for complex calls. | | rwork.set_size(ncv); // Work array of size NCV for complex calls. | |
| | | | |
| skipping to change at line 343 | | skipping to change at line 433 | |
|---|
| | | | |
| // IPARAM: integer array of length 11. | | // IPARAM: integer array of length 11. | |
| iparam.zeros(11); | | iparam.zeros(11); | |
| iparam(0) = 1; // Exact shifts (not provided by us). | | iparam(0) = 1; // Exact shifts (not provided by us). | |
| iparam(2) = 1000; // Maximum iterations; all the examples use 300, but
they were written in the ancient times. | | iparam(2) = 1000; // Maximum iterations; all the examples use 300, but
they were written in the ancient times. | |
| iparam(6) = 1; // Mode 1: A * x = lambda * x. | | iparam(6) = 1; // Mode 1: A * x = lambda * x. | |
| | | | |
| // IPNTR: integer array of length 14 (output). | | // IPNTR: integer array of length 14 (output). | |
| ipntr.set_size(14); | | ipntr.set_size(14); | |
| | | | |
|
| // Real work array used in the basic Arnoldi iteration for reverse | | // Real work array used in the basic Arnoldi iteration for reverse comm | |
| // communication. | | unication. | |
| workd.set_size(3 * n); | | workd.set_size(3 * n); | |
| | | | |
| // lworkl must be at least 3 * NCV^2 + 6 * NCV. | | // lworkl must be at least 3 * NCV^2 + 6 * NCV. | |
| lworkl = 3 * (ncv * ncv) + 6 * ncv; | | lworkl = 3 * (ncv * ncv) + 6 * ncv; | |
| | | | |
| // Real work array of length lworkl. | | // Real work array of length lworkl. | |
| workl.set_size(lworkl); | | workl.set_size(lworkl); | |
| | | | |
| info = 0; // Set to 0 initially to use random initial vector. | | info = 0; // Set to 0 initially to use random initial vector. | |
| | | | |
| // All the parameters have been set or created. Time to loop a lot. | | // All the parameters have been set or created. Time to loop a lot. | |
| while (ido != 99) | | while (ido != 99) | |
| { | | { | |
| // Call saupd() or naupd() with the current parameters. | | // Call saupd() or naupd() with the current parameters. | |
|
| if (sym) | | if(sym) | |
| arpack::saupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | | arpack::saupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | |
| else | | else | |
| arpack::naupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, rwork.memptr(), &info); | | arpack::naupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, rwork.memptr(), &info); | |
| | | | |
| // What do we do now? | | // What do we do now? | |
| switch (ido) | | switch (ido) | |
| { | | { | |
| case -1: | | case -1: | |
| case 1: | | case 1: | |
| { | | { | |
| | | | |
| End of changes. 29 change blocks. |
|---|
| 32 lines changed or deleted | | 126 lines changed or added | |
|---|
|
| subview_field_meat.hpp | | subview_field_meat.hpp | |
|---|
|
| // Copyright (C) 2008-2012 Conrad Sanderson | | // Copyright (C) 2008-2014 Conrad Sanderson | |
| // Copyright (C) 2008-2012 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2014 NICTA (www.nicta.com.au) | |
| // | | // | |
| // This Source Code Form is subject to the terms of the Mozilla Public | | // This Source Code Form is subject to the terms of the Mozilla Public | |
| // License, v. 2.0. If a copy of the MPL was not distributed with this | | // License, v. 2.0. If a copy of the MPL was not distributed with this | |
| // file, You can obtain one at http://mozilla.org/MPL/2.0/. | | // file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| | | | |
| //! \addtogroup subview_field | | //! \addtogroup subview_field | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| | | | |
| skipping to change at line 31 | | skipping to change at line 31 | |
|---|
| ( | | ( | |
| const field<oT>& in_f, | | const field<oT>& in_f, | |
| const uword in_row1, | | const uword in_row1, | |
| const uword in_col1, | | const uword in_col1, | |
| const uword in_n_rows, | | const uword in_n_rows, | |
| const uword in_n_cols | | const uword in_n_cols | |
| ) | | ) | |
| : f(in_f) | | : f(in_f) | |
| , aux_row1(in_row1) | | , aux_row1(in_row1) | |
| , aux_col1(in_col1) | | , aux_col1(in_col1) | |
|
| | | , aux_slice1(0) | |
| , n_rows(in_n_rows) | | , n_rows(in_n_rows) | |
| , n_cols(in_n_cols) | | , n_cols(in_n_cols) | |
|
| | | , n_slices(1) | |
| , n_elem(in_n_rows*in_n_cols) | | , n_elem(in_n_rows*in_n_cols) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
|
| | | arma_inline | |
| | | subview_field<oT>::subview_field | |
| | | ( | |
| | | const field<oT>& in_f, | |
| | | const uword in_row1, | |
| | | const uword in_col1, | |
| | | const uword in_slice1, | |
| | | const uword in_n_rows, | |
| | | const uword in_n_cols, | |
| | | const uword in_n_slices | |
| | | ) | |
| | | : f(in_f) | |
| | | , aux_row1(in_row1) | |
| | | , aux_col1(in_col1) | |
| | | , aux_slice1(in_slice1) | |
| | | , n_rows(in_n_rows) | |
| | | , n_cols(in_n_cols) | |
| | | , n_slices(in_n_slices) | |
| | | , n_elem(in_n_rows*in_n_cols*in_n_slices) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| subview_field<oT>::operator= (const field<oT>& x) | | subview_field<oT>::operator= (const field<oT>& x) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| subview_field<oT>& t = *this; | | subview_field<oT>& t = *this; | |
| | | | |
|
| arma_debug_check( (t.n_rows != x.n_rows) || (t.n_cols != x.n_cols), "inco
mpatible field dimensions"); | | arma_debug_check( (t.n_rows != x.n_rows) || (t.n_cols != x.n_cols) || (t.
n_slices != x.n_slices), "incompatible field dimensions"); | |
| | | | |
|
| for(uword col=0; col<t.n_cols; ++col) | | if(t.n_slices == 1) | |
| { | | { | |
|
| for(uword row=0; row<t.n_rows; ++row) | | for(uword col=0; col < t.n_cols; ++col) | |
| | | for(uword row=0; row < t.n_rows; ++row) | |
| { | | { | |
| t.at(row,col) = x.at(row,col); | | t.at(row,col) = x.at(row,col); | |
| } | | } | |
| } | | } | |
|
| | | else | |
| | | { | |
| | | for(uword slice=0; slice < t.n_slices; ++slice) | |
| | | for(uword col=0; col < t.n_cols; ++col ) | |
| | | for(uword row=0; row < t.n_rows; ++row ) | |
| | | { | |
| | | t.at(row,col,slice) = x.at(row,col,slice); | |
| | | } | |
| | | } | |
| } | | } | |
| | | | |
| //! x.subfield(...) = y.subfield(...) | | //! x.subfield(...) = y.subfield(...) | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| subview_field<oT>::operator= (const subview_field<oT>& x_in) | | subview_field<oT>::operator= (const subview_field<oT>& x_in) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const bool overlap = check_overlap(x_in); | | const bool overlap = check_overlap(x_in); | |
| | | | |
| field<oT>* tmp_field = overlap ? new field<oT>(x_in.f) :
0; | | field<oT>* tmp_field = overlap ? new field<oT>(x_in.f) :
0; | |
|
| const subview_field<oT>* tmp_subview = overlap ? new subview_field<oT>(*t
mp_field, x_in.aux_row1, x_in.aux_col1, x_in.n_rows, x_in.n_cols) : 0; | | const subview_field<oT>* tmp_subview = overlap ? new subview_field<oT>(*t
mp_field, x_in.aux_row1, x_in.aux_col1, x_in.aux_slice1, x_in.n_rows, x_in.
n_cols, x_in.n_slices) : 0; | |
| const subview_field<oT>& x = overlap ? (*tmp_subview) : x_in; | | const subview_field<oT>& x = overlap ? (*tmp_subview) : x_in; | |
| | | | |
| subview_field<oT>& t = *this; | | subview_field<oT>& t = *this; | |
| | | | |
|
| arma_debug_check( (t.n_rows != x.n_rows) || (t.n_cols != x.n_cols), "inco
mpatible field dimensions"); | | arma_debug_check( (t.n_rows != x.n_rows) || (t.n_cols != x.n_cols) || (t.
n_slices != x.n_slices), "incompatible field dimensions"); | |
| | | | |
|
| for(uword col=0; col<t.n_cols; ++col) | | if(t.n_slices == 1) | |
| { | | { | |
|
| for(uword row=0; row<t.n_rows; ++row) | | for(uword col=0; col < t.n_cols; ++col) | |
| | | for(uword row=0; row < t.n_rows; ++row) | |
| { | | { | |
| t.at(row,col) = x.at(row,col); | | t.at(row,col) = x.at(row,col); | |
| } | | } | |
| } | | } | |
|
| | | else | |
| | | { | |
| | | for(uword slice=0; slice < t.n_slices; ++slice) | |
| | | for(uword col=0; col < t.n_cols; ++col ) | |
| | | for(uword row=0; row < t.n_rows; ++row ) | |
| | | { | |
| | | t.at(row,col,slice) = x.at(row,col,slice); | |
| | | } | |
| | | } | |
| | | | |
| if(overlap) | | if(overlap) | |
| { | | { | |
| delete tmp_subview; | | delete tmp_subview; | |
| delete tmp_field; | | delete tmp_field; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
| subview_field<oT>::operator[](const uword i) | | subview_field<oT>::operator[](const uword i) | |
| { | | { | |
|
| const uword in_col = i / n_rows; | | uword index; | |
| const uword in_row = i % n_rows; | | | |
| | | | |
|
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | if(n_slices == 1) | |
| | | { | |
| | | const uword in_col = i / n_rows; | |
| | | const uword in_row = i % n_rows; | |
| | | | |
| | | index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | } | |
| | | else | |
| | | { | |
| | | const uword n_elem_slice = n_rows*n_cols; | |
| | | | |
| | | const uword in_slice = i / n_elem_slice; | |
| | | const uword offset = in_slice * n_elem_slice; | |
| | | const uword j = i - offset; | |
| | | | |
| | | const uword in_col = j / n_rows; | |
| | | const uword in_row = j % n_rows; | |
| | | | |
| | | index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col + aux_col | |
| | | 1)*f.n_rows + aux_row1 + in_row; | |
| | | } | |
| | | | |
| return *((const_cast< field<oT>& >(f)).mem[index]); | | return *((const_cast< field<oT>& >(f)).mem[index]); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| const oT& | | const oT& | |
| subview_field<oT>::operator[](const uword i) const | | subview_field<oT>::operator[](const uword i) const | |
| { | | { | |
|
| const uword in_col = i / n_rows; | | uword index; | |
| const uword in_row = i % n_rows; | | | |
| | | | |
|
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | if(n_slices == 1) | |
| | | { | |
| | | const uword in_col = i / n_rows; | |
| | | const uword in_row = i % n_rows; | |
| | | | |
| | | index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | } | |
| | | else | |
| | | { | |
| | | const uword n_elem_slice = n_rows*n_cols; | |
| | | | |
| | | const uword in_slice = i / n_elem_slice; | |
| | | const uword offset = in_slice * n_elem_slice; | |
| | | const uword j = i - offset; | |
| | | | |
| | | const uword in_col = j / n_rows; | |
| | | const uword in_row = j % n_rows; | |
| | | | |
| | | index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col + aux_col | |
| | | 1)*f.n_rows + aux_row1 + in_row; | |
| | | } | |
| | | | |
| return *(f.mem[index]); | | return *(f.mem[index]); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
| subview_field<oT>::operator()(const uword i) | | subview_field<oT>::operator()(const uword i) | |
| { | | { | |
| arma_debug_check( (i >= n_elem), "subview_field::operator(): index out of
bounds"); | | arma_debug_check( (i >= n_elem), "subview_field::operator(): index out of
bounds"); | |
| | | | |
|
| const uword in_col = i / n_rows; | | return operator[](i); | |
| const uword in_row = i % n_rows; | | | |
| | | | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | | |
| | | | |
| return *((const_cast< field<oT>& >(f)).mem[index]); | | | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| const oT& | | const oT& | |
| subview_field<oT>::operator()(const uword i) const | | subview_field<oT>::operator()(const uword i) const | |
| { | | { | |
| arma_debug_check( (i >= n_elem), "subview_field::operator(): index out of
bounds"); | | arma_debug_check( (i >= n_elem), "subview_field::operator(): index out of
bounds"); | |
| | | | |
|
| const uword in_col = i / n_rows; | | return operator[](i); | |
| const uword in_row = i % n_rows; | | | |
| | | | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | | |
| | | | |
| return *(f.mem[index]); | | | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
| subview_field<oT>::operator()(const uword in_row, const uword in_col) | | subview_field<oT>::operator()(const uword in_row, const uword in_col) | |
| { | | { | |
| arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "subview_fi
eld::operator(): index out of bounds"); | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "subview_fi
eld::operator(): index out of bounds"); | |
| | | | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| skipping to change at line 174 | | skipping to change at line 246 | |
|---|
| arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "subview_fi
eld::operator(): index out of bounds"); | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols)), "subview_fi
eld::operator(): index out of bounds"); | |
| | | | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| return *(f.mem[index]); | | return *(f.mem[index]); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| oT& | | oT& | |
|
| | | subview_field<oT>::operator()(const uword in_row, const uword in_col, const | |
| | | uword in_slice) | |
| | | { | |
| | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols) || (in_slice | |
| | | >= n_slices)), "subview_field::operator(): index out of bounds"); | |
| | | | |
| | | const uword index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col | |
| | | + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| | | return *((const_cast< field<oT>& >(f)).mem[index]); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | const oT& | |
| | | subview_field<oT>::operator()(const uword in_row, const uword in_col, const | |
| | | uword in_slice) const | |
| | | { | |
| | | arma_debug_check( ((in_row >= n_rows) || (in_col >= n_cols) || (in_slice | |
| | | >= n_slices)), "subview_field::operator(): index out of bounds"); | |
| | | | |
| | | const uword index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col | |
| | | + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| | | return *(f.mem[index]); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | oT& | |
| subview_field<oT>::at(const uword in_row, const uword in_col) | | subview_field<oT>::at(const uword in_row, const uword in_col) | |
| { | | { | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| return *((const_cast< field<oT>& >(f)).mem[index]); | | return *((const_cast< field<oT>& >(f)).mem[index]); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| arma_inline | | arma_inline | |
| const oT& | | const oT& | |
| subview_field<oT>::at(const uword in_row, const uword in_col) const | | subview_field<oT>::at(const uword in_row, const uword in_col) const | |
| { | | { | |
|
| //arma_extra_debug_sigprint(); | | | |
| | | | |
| const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | | const uword index = (in_col + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| return *(f.mem[index]); | | return *(f.mem[index]); | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
|
| | | arma_inline | |
| | | oT& | |
| | | subview_field<oT>::at(const uword in_row, const uword in_col, const uword i | |
| | | n_slice) | |
| | | { | |
| | | const uword index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col | |
| | | + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| | | return *((const_cast< field<oT>& >(f)).mem[index]); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| | | arma_inline | |
| | | const oT& | |
| | | subview_field<oT>::at(const uword in_row, const uword in_col, const uword i | |
| | | n_slice) const | |
| | | { | |
| | | const uword index = (in_slice + aux_slice1)*(f.n_rows*f.n_cols) + (in_col | |
| | | + aux_col1)*f.n_rows + aux_row1 + in_row; | |
| | | | |
| | | return *(f.mem[index]); | |
| | | } | |
| | | | |
| | | template<typename oT> | |
| inline | | inline | |
| bool | | bool | |
| subview_field<oT>::check_overlap(const subview_field<oT>& x) const | | subview_field<oT>::check_overlap(const subview_field<oT>& x) const | |
| { | | { | |
| const subview_field<oT>& t = *this; | | const subview_field<oT>& t = *this; | |
| | | | |
| if(&t.f != &x.f) | | if(&t.f != &x.f) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| else | | else | |
| { | | { | |
| if( (t.n_elem == 0) || (x.n_elem == 0) ) | | if( (t.n_elem == 0) || (x.n_elem == 0) ) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| const uword t_row_start = t.aux_row1; | | const uword t_row_start = t.aux_row1; | |
| const uword t_row_end_p1 = t_row_start + t.n_rows; | | const uword t_row_end_p1 = t_row_start + t.n_rows; | |
| | | | |
| | | const uword t_col_start = t.aux_col1; | |
| | | const uword t_col_end_p1 = t_col_start + t.n_cols; | |
| | | | |
|
| const uword t_col_start = t.aux_col1; | | const uword t_slice_start = t.aux_slice1; | |
| const uword t_col_end_p1 = t_col_start + t.n_cols; | | const uword t_slice_end_p1 = t_slice_start + t.n_slices; | |
| | | | |
|
| const uword x_row_start = x.aux_row1; | | const uword x_row_start = x.aux_row1; | |
| const uword x_row_end_p1 = x_row_start + x.n_rows; | | const uword x_row_end_p1 = x_row_start + x.n_rows; | |
| | | | |
|
| const uword x_col_start = x.aux_col1; | | const uword x_col_start = x.aux_col1; | |
| const uword x_col_end_p1 = x_col_start + x.n_cols; | | const uword x_col_end_p1 = x_col_start + x.n_cols; | |
| | | | |
|
| const bool outside_rows = ( (x_row_start >= t_row_end_p1) || (t_row_s | | const uword x_slice_start = x.aux_slice1; | |
| tart >= x_row_end_p1) ); | | const uword x_slice_end_p1 = x_slice_start + x.n_slices; | |
| const bool outside_cols = ( (x_col_start >= t_col_end_p1) || (t_col_s | | | |
| tart >= x_col_end_p1) ); | | | |
| | | | |
|
| return ( (outside_rows == false) && (outside_cols == false) ); | | const bool outside_rows = ( (x_row_start >= t_row_end_p1 ) || (t | |
| | | _row_start >= x_row_end_p1 ) ); | |
| | | const bool outside_cols = ( (x_col_start >= t_col_end_p1 ) || (t | |
| | | _col_start >= x_col_end_p1 ) ); | |
| | | const bool outside_slices = ( (x_slice_start >= t_slice_end_p1) || (t | |
| | | _slice_start >= x_slice_end_p1) ); | |
| | | | |
| | | return ( (outside_rows == false) && (outside_cols == false) && (outsi | |
| | | de_slices == false) ); | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename oT> | | template<typename oT> | |
| inline | | inline | |
| void | | void | |
| subview_field<oT>::print(const std::string extra_text) const | | subview_field<oT>::print(const std::string extra_text) const | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| skipping to change at line 286 | | skipping to change at line 407 | |
|---|
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| // | | // | |
| const bool alias = (&actual_out == &in.f); | | const bool alias = (&actual_out == &in.f); | |
| | | | |
| field<oT>* tmp = (alias) ? new field<oT> : 0; | | field<oT>* tmp = (alias) ? new field<oT> : 0; | |
| field<oT>& out = (alias) ? (*tmp) : actual_out; | | field<oT>& out = (alias) ? (*tmp) : actual_out; | |
| | | | |
| // | | // | |
| | | | |
|
| const uword n_rows = in.n_rows; | | const uword n_rows = in.n_rows; | |
| const uword n_cols = in.n_cols; | | const uword n_cols = in.n_cols; | |
| | | const uword n_slices = in.n_slices; | |
| | | | |
|
| out.set_size(n_rows, n_cols); | | out.set_size(n_rows, n_cols, n_slices); | |
| | | | |
|
| arma_extra_debug_print(arma_boost::format("out.n_rows = %d out.n_cols =
%d in.m.n_rows = %d in.m.n_cols = %d") % out.n_rows % out.n_cols % in.
f.n_rows % in.f.n_cols ); | | arma_extra_debug_print(arma_boost::format("out.n_rows = %d out.n_cols =
%d out.n_slices = %d in.m.n_rows = %d in.m.n_cols = %d in.m.n_slice
s = %d") % out.n_rows % out.n_cols % out.n_slices % in.f.n_rows % in.f.n_co
ls % in.f.n_slices); | |
| | | | |
|
| for(uword col = 0; col<n_cols; ++col) | | if(n_slices == 1) | |
| { | | { | |
|
| for(uword row = 0; row<n_rows; ++row) | | for(uword col = 0; col < n_cols; ++col) | |
| | | for(uword row = 0; row < n_rows; ++row) | |
| { | | { | |
| out.at(row,col) = in.at(row,col); | | out.at(row,col) = in.at(row,col); | |
| } | | } | |
| } | | } | |
|
| | | else | |
| | | { | |
| | | for(uword slice = 0; slice < n_slices; ++slice) | |
| | | for(uword col = 0; col < n_cols; ++col ) | |
| | | for(uword row = 0; row < n_rows; ++row ) | |
| | | { | |
| | | out.at(row,col,slice) = in.at(row,col,slice); | |
| | | } | |
| | | } | |
| | | | |
| if(alias) | | if(alias) | |
| { | | { | |
| actual_out = out; | | actual_out = out; | |
| delete tmp; | | delete tmp; | |
| } | | } | |
| | | | |
| } | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
| End of changes. 34 change blocks. |
|---|
| 48 lines changed or deleted | | 194 lines changed or added | |
|---|
|