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
20 //! \addtogroup op_orth_null
21 //! @{
22
23
24
25 template<typename T1>
26 inline
27 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_orth> & expr)28 op_orth::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_orth>& expr)
29 {
30 arma_extra_debug_sigprint();
31
32 typedef typename T1::pod_type T;
33
34 const T tol = access::tmp_real(expr.aux);
35
36 const bool status = op_orth::apply_direct(out, expr.m, tol);
37
38 if(status == false)
39 {
40 out.soft_reset();
41 arma_stop_runtime_error("orth(): svd failed");
42 }
43 }
44
45
46
47 template<typename T1>
48 inline
49 bool
apply_direct(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & expr,typename T1::pod_type tol)50 op_orth::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr, typename T1::pod_type tol)
51 {
52 arma_extra_debug_sigprint();
53
54 typedef typename T1::elem_type eT;
55 typedef typename T1::pod_type T;
56
57 arma_debug_check((tol < T(0)), "orth(): tolerance must be >= 0");
58
59 Mat<eT> A(expr.get_ref());
60
61 Mat<eT> U;
62 Col< T> s;
63 Mat<eT> V;
64
65 const bool status = auxlib::svd_dc(U, s, V, A);
66
67 V.reset();
68
69 if(status == false) { return false; }
70
71 if(s.is_empty()) { out.reset(); return true; }
72
73 const uword s_n_elem = s.n_elem;
74 const T* s_mem = s.memptr();
75
76 // set tolerance to default if it hasn't been specified
77 if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * s_mem[0] * std::numeric_limits<T>::epsilon(); }
78
79 uword count = 0;
80
81 for(uword i=0; i < s_n_elem; ++i) { count += (s_mem[i] > tol) ? uword(1) : uword(0); }
82
83 if(count > 0)
84 {
85 out = U.head_cols(count); // out *= eT(-1);
86 }
87 else
88 {
89 out.set_size(A.n_rows, 0);
90 }
91
92 return true;
93 }
94
95
96
97 //
98
99
100
101 template<typename T1>
102 inline
103 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_null> & expr)104 op_null::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_null>& expr)
105 {
106 arma_extra_debug_sigprint();
107
108 typedef typename T1::pod_type T;
109
110 const T tol = access::tmp_real(expr.aux);
111
112 const bool status = op_null::apply_direct(out, expr.m, tol);
113
114 if(status == false)
115 {
116 out.soft_reset();
117 arma_stop_runtime_error("null(): svd failed");
118 }
119 }
120
121
122
123 template<typename T1>
124 inline
125 bool
apply_direct(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & expr,typename T1::pod_type tol)126 op_null::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr, typename T1::pod_type tol)
127 {
128 arma_extra_debug_sigprint();
129
130 typedef typename T1::elem_type eT;
131 typedef typename T1::pod_type T;
132
133 arma_debug_check((tol < T(0)), "null(): tolerance must be >= 0");
134
135 Mat<eT> A(expr.get_ref());
136
137 Mat<eT> U;
138 Col< T> s;
139 Mat<eT> V;
140
141 const bool status = auxlib::svd_dc(U, s, V, A);
142
143 U.reset();
144
145 if(status == false) { return false; }
146
147 if(s.is_empty()) { out.reset(); return true; }
148
149 const uword s_n_elem = s.n_elem;
150 const T* s_mem = s.memptr();
151
152 // set tolerance to default if it hasn't been specified
153 if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * s_mem[0] * std::numeric_limits<T>::epsilon(); }
154
155 uword count = 0;
156
157 for(uword i=0; i < s_n_elem; ++i) { count += (s_mem[i] > tol) ? uword(1) : uword(0); }
158
159 if(count < A.n_cols)
160 {
161 out = V.tail_cols(A.n_cols - count);
162
163 const uword out_n_elem = out.n_elem;
164 eT* out_mem = out.memptr();
165
166 for(uword i=0; i<out_n_elem; ++i)
167 {
168 if(std::abs(out_mem[i]) < std::numeric_limits<T>::epsilon()) { out_mem[i] = eT(0); }
169 }
170 }
171 else
172 {
173 out.set_size(A.n_cols, 0);
174 }
175
176 return true;
177 }
178
179
180
181 //! @}
182