1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup spglue_min
20 //! @{
21
22
23
24 template<typename T1, typename T2>
25 inline
26 void
apply(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_min> & X)27 spglue_min::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_min>& X)
28 {
29 arma_extra_debug_sigprint();
30
31 typedef typename T1::elem_type eT;
32
33 const SpProxy<T1> pa(X.A);
34 const SpProxy<T2> pb(X.B);
35
36 const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
37
38 if(is_alias == false)
39 {
40 spglue_min::apply_noalias(out, pa, pb);
41 }
42 else
43 {
44 SpMat<eT> tmp;
45
46 spglue_min::apply_noalias(tmp, pa, pb);
47
48 out.steal_mem(tmp);
49 }
50 }
51
52
53
54 template<typename eT, typename T1, typename T2>
55 inline
56 void
apply_noalias(SpMat<eT> & out,const SpProxy<T1> & pa,const SpProxy<T2> & pb)57 spglue_min::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
58 {
59 arma_extra_debug_sigprint();
60
61 arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "element-wise minimum");
62
63 const uword max_n_nonzero = pa.get_n_nonzero() + pb.get_n_nonzero();
64
65 // Resize memory to upper bound
66 out.reserve(pa.get_n_rows(), pa.get_n_cols(), max_n_nonzero);
67
68 // Now iterate across both matrices.
69 typename SpProxy<T1>::const_iterator_type x_it = pa.begin();
70 typename SpProxy<T1>::const_iterator_type x_end = pa.end();
71
72 typename SpProxy<T2>::const_iterator_type y_it = pb.begin();
73 typename SpProxy<T2>::const_iterator_type y_end = pb.end();
74
75 uword count = 0;
76
77 while( (x_it != x_end) || (y_it != y_end) )
78 {
79 eT out_val;
80
81 const uword x_it_col = x_it.col();
82 const uword x_it_row = x_it.row();
83
84 const uword y_it_col = y_it.col();
85 const uword y_it_row = y_it.row();
86
87 bool use_y_loc = false;
88
89 if(x_it == y_it)
90 {
91 out_val = elem_min(eT(*x_it), eT(*y_it));
92
93 ++x_it;
94 ++y_it;
95 }
96 else
97 {
98 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
99 {
100 out_val = elem_min(eT(*x_it), eT(0));
101
102 ++x_it;
103 }
104 else
105 {
106 out_val = elem_min(eT(*y_it), eT(0));
107
108 ++y_it;
109
110 use_y_loc = true;
111 }
112 }
113
114 if(out_val != eT(0))
115 {
116 access::rw(out.values[count]) = out_val;
117
118 const uword out_row = (use_y_loc == false) ? x_it_row : y_it_row;
119 const uword out_col = (use_y_loc == false) ? x_it_col : y_it_col;
120
121 access::rw(out.row_indices[count]) = out_row;
122 access::rw(out.col_ptrs[out_col + 1])++;
123 ++count;
124 }
125
126 arma_check( (count > max_n_nonzero), "internal error: spglue_min::apply_noalias(): count > max_n_nonzero" );
127 }
128
129 const uword out_n_cols = out.n_cols;
130
131 uword* col_ptrs = access::rwp(out.col_ptrs);
132
133 // Fix column pointers to be cumulative.
134 for(uword c = 1; c <= out_n_cols; ++c)
135 {
136 col_ptrs[c] += col_ptrs[c - 1];
137 }
138
139 if(count < max_n_nonzero)
140 {
141 if(count <= (max_n_nonzero/2))
142 {
143 out.mem_resize(count);
144 }
145 else
146 {
147 // quick resize without reallocating memory and copying data
148 access::rw( out.n_nonzero) = count;
149 access::rw( out.values[count]) = eT(0);
150 access::rw(out.row_indices[count]) = uword(0);
151 }
152 }
153 }
154
155
156
157 template<typename eT>
158 inline
159 void
apply_noalias(SpMat<eT> & out,const SpMat<eT> & A,const SpMat<eT> & B)160 spglue_min::apply_noalias(SpMat<eT>& out, const SpMat<eT>& A, const SpMat<eT>& B)
161 {
162 arma_extra_debug_sigprint();
163
164 const SpProxy< SpMat<eT> > pa(A);
165 const SpProxy< SpMat<eT> > pb(B);
166
167 spglue_min::apply_noalias(out, pa, pb);
168 }
169
170
171
172 template<typename eT, typename T1, typename T2>
173 inline
174 void
dense_sparse_min(Mat<eT> & out,const Base<eT,T1> & X,const SpBase<eT,T2> & Y)175 spglue_min::dense_sparse_min(Mat<eT>& out, const Base<eT,T1>& X, const SpBase<eT,T2>& Y)
176 {
177 arma_extra_debug_sigprint();
178
179 // NOTE: this function assumes there is no aliasing between matrix 'out' and X
180
181 const Proxy<T1> pa(X.get_ref());
182 const SpProxy<T2> pb(Y.get_ref());
183
184 const uword n_rows = pa.get_n_rows();
185 const uword n_cols = pa.get_n_cols();
186
187 arma_debug_assert_same_size( n_rows, n_cols, pb.get_n_rows(), pb.get_n_cols(), "element-wise minimum" );
188
189 out.set_size(n_rows, n_cols);
190
191 for(uword c=0; c < n_cols; ++c)
192 for(uword r=0; r < n_rows; ++r)
193 {
194 out.at(r,c) = elem_min(pa.at(r,c), pb.at(r,c));
195 }
196 }
197
198
199
200 // min of non-complex elements
201 template<typename eT>
202 inline
203 typename enable_if2<is_cx<eT>::no, eT>::result
elem_min(const eT & a,const eT & b)204 spglue_min::elem_min(const eT& a, const eT& b)
205 {
206 return (std::min)(a, b);
207 }
208
209
210
211 // min of complex elements
212 template<typename eT>
213 inline
214 typename enable_if2<is_cx<eT>::yes, eT>::result
elem_min(const eT & a,const eT & b)215 spglue_min::elem_min(const eT& a, const eT& b)
216 {
217 return (std::abs(a) < std::abs(b)) ? a : b;
218 }
219
220
221
222 //! @}
223