annotate armadillo-3.900.4/include/armadillo_bits/spglue_minus_meat.hpp @ 84:55a047986812 tip

Update library URI so as not to be document-local
author Chris Cannam
date Wed, 22 Apr 2020 14:21:57 +0100
parents 1ec0e2823891
children
rev   line source
Chris@49 1 // Copyright (C) 2012 Ryan Curtin
Chris@49 2 // Copyright (C) 2012 Conrad Sanderson
Chris@49 3 //
Chris@49 4 // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49 5 // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49 6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49 7
Chris@49 8
Chris@49 9 //! \addtogroup spglue_minus
Chris@49 10 //! @{
Chris@49 11
Chris@49 12
Chris@49 13
Chris@49 14 template<typename T1, typename T2>
Chris@49 15 inline
Chris@49 16 void
Chris@49 17 spglue_minus::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus>& X)
Chris@49 18 {
Chris@49 19 arma_extra_debug_sigprint();
Chris@49 20
Chris@49 21 typedef typename T1::elem_type eT;
Chris@49 22
Chris@49 23 const SpProxy<T1> pa(X.A);
Chris@49 24 const SpProxy<T2> pb(X.B);
Chris@49 25
Chris@49 26 const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
Chris@49 27
Chris@49 28 if(is_alias == false)
Chris@49 29 {
Chris@49 30 spglue_minus::apply_noalias(out, pa, pb);
Chris@49 31 }
Chris@49 32 else
Chris@49 33 {
Chris@49 34 SpMat<eT> tmp;
Chris@49 35 spglue_minus::apply_noalias(tmp, pa, pb);
Chris@49 36
Chris@49 37 out.steal_mem(tmp);
Chris@49 38 }
Chris@49 39 }
Chris@49 40
Chris@49 41
Chris@49 42
Chris@49 43 template<typename eT, typename T1, typename T2>
Chris@49 44 arma_hot
Chris@49 45 inline
Chris@49 46 void
Chris@49 47 spglue_minus::apply_noalias(SpMat<eT>& result, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
Chris@49 48 {
Chris@49 49 arma_extra_debug_sigprint();
Chris@49 50
Chris@49 51 arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "subtraction");
Chris@49 52
Chris@49 53 if( (pa.get_n_nonzero() != 0) && (pb.get_n_nonzero() != 0) )
Chris@49 54 {
Chris@49 55 result.set_size(pa.get_n_rows(), pa.get_n_cols());
Chris@49 56
Chris@49 57 // Resize memory to correct size.
Chris@49 58 result.mem_resize(n_unique(pa, pb, op_n_unique_sub()));
Chris@49 59
Chris@49 60 // Now iterate across both matrices.
Chris@49 61 typename SpProxy<T1>::const_iterator_type x_it = pa.begin();
Chris@49 62 typename SpProxy<T2>::const_iterator_type y_it = pb.begin();
Chris@49 63
Chris@49 64 typename SpProxy<T1>::const_iterator_type x_end = pa.end();
Chris@49 65 typename SpProxy<T2>::const_iterator_type y_end = pb.end();
Chris@49 66
Chris@49 67 uword cur_val = 0;
Chris@49 68 while((x_it != x_end) || (y_it != y_end))
Chris@49 69 {
Chris@49 70 if(x_it == y_it)
Chris@49 71 {
Chris@49 72 const eT val = (*x_it) - (*y_it);
Chris@49 73
Chris@49 74 if (val != eT(0))
Chris@49 75 {
Chris@49 76 access::rw(result.values[cur_val]) = val;
Chris@49 77 access::rw(result.row_indices[cur_val]) = x_it.row();
Chris@49 78 ++access::rw(result.col_ptrs[x_it.col() + 1]);
Chris@49 79 ++cur_val;
Chris@49 80 }
Chris@49 81
Chris@49 82 ++x_it;
Chris@49 83 ++y_it;
Chris@49 84 }
Chris@49 85 else
Chris@49 86 {
Chris@49 87 const uword x_it_row = x_it.row();
Chris@49 88 const uword x_it_col = x_it.col();
Chris@49 89
Chris@49 90 const uword y_it_row = y_it.row();
Chris@49 91 const uword y_it_col = y_it.col();
Chris@49 92
Chris@49 93 if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_it_row))) // if y is closer to the end
Chris@49 94 {
Chris@49 95 access::rw(result.values[cur_val]) = (*x_it);
Chris@49 96 access::rw(result.row_indices[cur_val]) = x_it_row;
Chris@49 97 ++access::rw(result.col_ptrs[x_it_col + 1]);
Chris@49 98 ++cur_val;
Chris@49 99 ++x_it;
Chris@49 100 }
Chris@49 101 else
Chris@49 102 {
Chris@49 103 access::rw(result.values[cur_val]) = -(*y_it);
Chris@49 104 access::rw(result.row_indices[cur_val]) = y_it_row;
Chris@49 105 ++access::rw(result.col_ptrs[y_it_col + 1]);
Chris@49 106 ++cur_val;
Chris@49 107 ++y_it;
Chris@49 108 }
Chris@49 109 }
Chris@49 110 }
Chris@49 111
Chris@49 112 // Fix column pointers to be cumulative.
Chris@49 113 for(uword c = 1; c <= result.n_cols; ++c)
Chris@49 114 {
Chris@49 115 access::rw(result.col_ptrs[c]) += result.col_ptrs[c - 1];
Chris@49 116 }
Chris@49 117 }
Chris@49 118 else
Chris@49 119 {
Chris@49 120 if(pa.get_n_nonzero() == 0)
Chris@49 121 {
Chris@49 122 result = pb.Q;
Chris@49 123 result *= eT(-1);
Chris@49 124
Chris@49 125 return;
Chris@49 126 }
Chris@49 127
Chris@49 128 if(pb.get_n_nonzero() == 0)
Chris@49 129 {
Chris@49 130 result = pa.Q;
Chris@49 131 return;
Chris@49 132 }
Chris@49 133 }
Chris@49 134 }
Chris@49 135
Chris@49 136
Chris@49 137
Chris@49 138 //
Chris@49 139 //
Chris@49 140 // spglue_minus2: scalar*(A - B)
Chris@49 141
Chris@49 142
Chris@49 143
Chris@49 144 template<typename T1, typename T2>
Chris@49 145 inline
Chris@49 146 void
Chris@49 147 spglue_minus2::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus2>& X)
Chris@49 148 {
Chris@49 149 arma_extra_debug_sigprint();
Chris@49 150
Chris@49 151 typedef typename T1::elem_type eT;
Chris@49 152
Chris@49 153 const SpProxy<T1> pa(X.A);
Chris@49 154 const SpProxy<T2> pb(X.B);
Chris@49 155
Chris@49 156 const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
Chris@49 157
Chris@49 158 if(is_alias == false)
Chris@49 159 {
Chris@49 160 spglue_minus::apply_noalias(out, pa, pb);
Chris@49 161 }
Chris@49 162 else
Chris@49 163 {
Chris@49 164 SpMat<eT> tmp;
Chris@49 165 spglue_minus::apply_noalias(tmp, pa, pb);
Chris@49 166
Chris@49 167 out.steal_mem(tmp);
Chris@49 168 }
Chris@49 169
Chris@49 170 out *= X.aux;
Chris@49 171 }
Chris@49 172
Chris@49 173
Chris@49 174
Chris@49 175 //! @}