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 auxlib
20 //! @{
21 
22 
23 //! interface functions for accessing decompositions in LAPACK and ATLAS
24 class auxlib
25   {
26   public:
27 
28   //
29   // inv
30 
31   template<typename eT>
32   inline static bool inv(Mat<eT>& A);
33 
34   template<typename eT>
35   inline static bool inv(Mat<eT>& out, const Mat<eT>& X);
36 
37   template<typename eT>
38   inline static bool inv_tr(Mat<eT>& A, const uword layout);
39 
40   template<typename eT>
41   inline static bool inv_sympd(Mat<eT>& A);
42 
43   template<typename eT>
44   inline static bool inv_sympd(Mat<eT>& out, const Mat<eT>& X);
45 
46   template<typename eT>
47   inline static bool inv_sympd_rcond(Mat<eT>& A, const eT rcond_threshold);
48 
49   template<typename T>
50   inline static bool inv_sympd_rcond(Mat< std::complex<T> >& A, const T rcond_threshold);
51 
52 
53   //
54   // det and log_det
55 
56   template<typename eT>
57   inline static bool det(eT& out_val, Mat<eT>& A);
58 
59   template<typename eT>
60   inline static bool log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, Mat<eT>& A);
61 
62   template<typename eT>
63   inline static bool log_det_sympd(typename get_pod_type<eT>::result& out_val, Mat<eT>& A);
64 
65 
66   //
67   // lu
68 
69   template<typename eT, typename T1>
70   inline static bool lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>& X);
71 
72   template<typename eT, typename T1>
73   inline static bool lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X);
74 
75   template<typename eT, typename T1>
76   inline static bool lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X);
77 
78 
79   //
80   // eig_gen
81 
82   template<typename T1>
83   inline static bool eig_gen(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base<typename T1::pod_type,T1>& expr);
84 
85   template<typename T1>
86   inline static bool eig_gen(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base< std::complex<typename T1::pod_type>, T1 >& expr);
87 
88 
89   //
90   // eig_gen_balance
91 
92   template<typename T1>
93   inline static bool eig_gen_balance(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base<typename T1::pod_type,T1>& expr);
94 
95   template<typename T1>
96   inline static bool eig_gen_balance(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base< std::complex<typename T1::pod_type>, T1 >& expr);
97 
98 
99   //
100   // eig_gen_twosided
101 
102   template<typename T1>
103   inline static bool eig_gen_twosided(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base<typename T1::pod_type,T1>& expr);
104 
105   template<typename T1>
106   inline static bool eig_gen_twosided(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base< std::complex<typename T1::pod_type>, T1 >& expr);
107 
108 
109   //
110   // eig_gen_twosided_balance
111 
112   template<typename T1>
113   inline static bool eig_gen_twosided_balance(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base<typename T1::pod_type,T1>& expr);
114 
115   template<typename T1>
116   inline static bool eig_gen_twosided_balance(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base< std::complex<typename T1::pod_type>, T1 >& expr);
117 
118 
119   //
120   // eig_pair
121 
122   template<typename T1, typename T2>
123   inline static bool eig_pair(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base<typename T1::pod_type,T1>& A_expr, const Base<typename T1::pod_type,T2>& B_expr);
124 
125   template<typename T1, typename T2>
126   inline static bool eig_pair(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& vecs, const bool vecs_on, const Base< std::complex<typename T1::pod_type>, T1 >& A_expr, const Base< std::complex<typename T1::pod_type>, T2 >& B_expr);
127 
128 
129   //
130   // eig_pair_twosided
131 
132   template<typename T1, typename T2>
133   inline static bool eig_pair_twosided(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base<typename T1::pod_type,T1>& A_expr, const Base<typename T1::pod_type,T2>& B_expr);
134 
135   template<typename T1, typename T2>
136   inline static bool eig_pair_twosided(Mat< std::complex<typename T1::pod_type> >& vals, Mat< std::complex<typename T1::pod_type> >& lvecs, Mat< std::complex<typename T1::pod_type> >& rvecs, const Base< std::complex<typename T1::pod_type>, T1 >& A_expr, const Base< std::complex<typename T1::pod_type>, T2 >& B_expr);
137 
138 
139   //
140   // eig_sym
141 
142   template<typename eT, typename T1>
143   inline static bool eig_sym(Col<eT>& eigval, const Base<eT,T1>& X);
144 
145   template<typename T, typename T1>
146   inline static bool eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1>& X);
147 
148   template<typename eT>
149   inline static bool eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X);
150 
151   template<typename T>
152   inline static bool eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< std::complex<T> >& X);
153 
154   template<typename eT>
155   inline static bool eig_sym_dc(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X);
156 
157   template<typename T>
158   inline static bool eig_sym_dc(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< std::complex<T> >& X);
159 
160 
161   //
162   // chol
163 
164   template<typename eT>
165   inline static bool chol_simple(Mat<eT>& X);
166 
167   template<typename eT>
168   inline static bool chol(Mat<eT>& X, const uword layout);
169 
170   template<typename eT>
171   inline static bool chol_band(Mat<eT>& X, const uword KD, const uword layout);
172 
173   template<typename  T>
174   inline static bool chol_band(Mat< std::complex<T> >& X, const uword KD, const uword layout);
175 
176   template<typename eT>
177   inline static bool chol_band_common(Mat<eT>& X, const uword KD, const uword layout);
178 
179   template<typename eT>
180   inline static bool chol_pivot(Mat<eT>& X, Mat<uword>& P, const uword layout);
181 
182 
183   //
184   // hessenberg decomposition
185 
186   template<typename eT, typename T1>
187   inline static bool hess(Mat<eT>& H, const Base<eT,T1>& X, Col<eT>& tao);
188 
189 
190   //
191   // qr
192 
193   template<typename eT, typename T1>
194   inline static bool qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X);
195 
196   template<typename eT, typename T1>
197   inline static bool qr_econ(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X);
198 
199   template<typename eT, typename T1>
200   inline static bool qr_pivot(Mat<eT>& Q, Mat<eT>& R, Mat<uword>& P, const Base<eT,T1>& X);
201 
202   template<typename  T, typename T1>
203   inline static bool qr_pivot(Mat< std::complex<T> >& Q, Mat< std::complex<T> >& R, Mat<uword>& P, const Base<std::complex<T>,T1>& X);
204 
205 
206   //
207   // svd
208 
209   template<typename eT>
210   inline static bool svd(Col<eT>& S, Mat<eT>& A);
211 
212   template<typename T>
213   inline static bool svd(Col<T>& S, Mat< std::complex<T> >& A);
214 
215 
216   template<typename eT>
217   inline static bool svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A);
218 
219   template<typename T>
220   inline static bool svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat< std::complex<T> >& A);
221 
222   template<typename eT>
223   inline static bool svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A, const char mode);
224 
225   template<typename T>
226   inline static bool svd_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat< std::complex<T> >& A, const char mode);
227 
228 
229   template<typename eT>
230   inline static bool svd_dc(Col<eT>& S, Mat<eT>& A);
231 
232   template<typename T>
233   inline static bool svd_dc(Col<T>& S, Mat< std::complex<T> >& A);
234 
235 
236   template<typename eT>
237   inline static bool svd_dc(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A);
238 
239   template<typename T>
240   inline static bool svd_dc(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat< std::complex<T> >& A);
241 
242   template<typename eT>
243   inline static bool svd_dc_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A);
244 
245   template<typename T>
246   inline static bool svd_dc_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat< std::complex<T> >& A);
247 
248 
249   //
250   // solve
251 
252   template<typename T1>
253   arma_cold inline static bool solve_square_tiny(Mat<typename T1::elem_type>& out, const Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
254 
255   template<typename T1>
256   inline static bool solve_square_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
257 
258   template<typename T1>
259   inline static bool solve_square_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr, const bool allow_ugly);
260 
261   template<typename T1>
262   inline static bool solve_square_refine(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
263 
264   template<typename T1>
265   inline static bool solve_square_refine(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const Base<std::complex<typename T1::pod_type>,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
266 
267   //
268 
269   template<typename T1>
270   inline static bool solve_sympd_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
271 
272   template<typename T1>
273   inline static bool solve_sympd_fast_common(Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
274 
275   template<typename T1>
276   inline static bool solve_sympd_rcond(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr, const bool allow_ugly);
277 
278   template<typename T1>
279   inline static bool solve_sympd_rcond(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const Base< std::complex<typename T1::pod_type>,T1>& B_expr, const bool allow_ugly);
280 
281   template<typename T1>
282   inline static bool solve_sympd_refine(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
283 
284   template<typename T1>
285   inline static bool solve_sympd_refine(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const Base<std::complex<typename T1::pod_type>,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
286 
287   //
288 
289   template<typename T1>
290   inline static bool solve_rect_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
291 
292   template<typename T1>
293   inline static bool solve_rect_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr, const bool allow_ugly);
294 
295   //
296 
297   template<typename T1>
298   inline static bool solve_approx_svd(Mat<typename T1::pod_type>& out, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr);
299 
300   template<typename T1>
301   inline static bool solve_approx_svd(Mat< std::complex<typename T1::pod_type> >& out, Mat< std::complex<typename T1::pod_type> >& A, const Base<std::complex<typename T1::pod_type>,T1>& B_expr);
302 
303   //
304 
305   template<typename T1>
306   inline static bool solve_trimat_fast(Mat<typename T1::elem_type>& out, const Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr, const uword layout);
307 
308   template<typename T1>
309   inline static bool solve_trimat_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_type& out_rcond, const Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr, const uword layout, const bool allow_ugly);
310 
311   //
312 
313   template<typename T1>
314   inline static bool solve_band_fast(Mat<typename T1::pod_type>& out, Mat<typename T1::pod_type>& A, const uword KL, const uword KU, const Base<typename T1::pod_type,T1>& B_expr);
315 
316   template<typename T1>
317   inline static bool solve_band_fast(Mat< std::complex<typename T1::pod_type> >& out, Mat< std::complex<typename T1::pod_type> >& A, const uword KL, const uword KU, const Base< std::complex<typename T1::pod_type>,T1>& B_expr);
318 
319   template<typename T1>
320   inline static bool solve_band_fast_common(Mat<typename T1::elem_type>& out, const Mat<typename T1::elem_type>& A, const uword KL, const uword KU, const Base<typename T1::elem_type,T1>& B_expr);
321 
322   template<typename T1>
323   inline static bool solve_band_rcond(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const uword KL, const uword KU, const Base<typename T1::pod_type,T1>& B_expr, const bool allow_ugly);
324 
325   template<typename T1>
326   inline static bool solve_band_rcond(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const uword KL, const uword KU, const Base< std::complex<typename T1::pod_type>,T1>& B_expr, const bool allow_ugly);
327 
328   template<typename T1>
329   inline static bool solve_band_rcond_common(Mat<typename T1::elem_type>& out, typename T1::pod_type& out_rcond, const Mat<typename T1::elem_type>& A, const uword KL, const uword KU, const Base<typename T1::elem_type,T1>& B_expr, const bool allow_ugly);
330 
331   template<typename T1>
332   inline static bool solve_band_refine(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const uword KL, const uword KU, const Base<typename T1::pod_type,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
333 
334   template<typename T1>
335   inline static bool solve_band_refine(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const uword KL, const uword KU, const Base<std::complex<typename T1::pod_type>,T1>& B_expr, const bool equilibrate, const bool allow_ugly);
336 
337   //
338 
339   template<typename T1>
340   inline static bool solve_tridiag_fast(Mat<typename T1::pod_type>& out, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr);
341 
342   template<typename T1>
343   inline static bool solve_tridiag_fast(Mat< std::complex<typename T1::pod_type> >& out, Mat< std::complex<typename T1::pod_type> >& A, const Base< std::complex<typename T1::pod_type>,T1>& B_expr);
344 
345   template<typename T1>
346   inline static bool solve_tridiag_fast_common(Mat<typename T1::elem_type>& out, const Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
347 
348 
349   //
350   // Schur decomposition
351 
352   template<typename eT, typename T1>
353   inline static bool schur(Mat<eT>& U, Mat<eT>& S, const Base<eT,T1>& X, const bool calc_U = true);
354 
355   template<typename  T, typename T1>
356   inline static bool schur(Mat< std::complex<T> >& U, Mat< std::complex<T> >& S, const Base<std::complex<T>,T1>& X, const bool calc_U = true);
357 
358   template<typename  T>
359   inline static bool schur(Mat< std::complex<T> >& U, Mat< std::complex<T> >& S, const bool calc_U = true);
360 
361   //
362   // solve the Sylvester equation AX + XB = C
363 
364   template<typename eT>
365   inline static bool syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C);
366 
367 
368   //
369   // QZ decomposition
370 
371   template<typename T, typename T1, typename T2>
372   inline static bool qz(Mat<T>& A, Mat<T>& B, Mat<T>& vsl, Mat<T>& vsr, const Base<T,T1>& X_expr, const Base<T,T2>& Y_expr, const char mode);
373 
374   template<typename T, typename T1, typename T2>
375   inline static bool qz(Mat< std::complex<T> >& A, Mat< std::complex<T> >& B, Mat< std::complex<T> >& vsl, Mat< std::complex<T> >& vsr, const Base< std::complex<T>, T1 >& X_expr, const Base< std::complex<T>, T2 >& Y_expr, const char mode);
376 
377 
378   //
379   // rcond
380 
381   template<typename eT>
382   inline static eT rcond(Mat<eT>& A);
383 
384   template<typename  T>
385   inline static  T rcond(Mat< std::complex<T> >& A);
386 
387   template<typename eT>
388   inline static eT rcond_sympd(Mat<eT>& A, bool& calc_ok);
389 
390   template<typename T>
391   inline static  T rcond_sympd(Mat< std::complex<T> >& A, bool& calc_ok);
392 
393   template<typename eT>
394   inline static eT rcond_trimat(const Mat<eT>& A, const uword layout);
395 
396   template<typename  T>
397   inline static  T rcond_trimat(const Mat< std::complex<T> >& A, const uword layout);
398 
399 
400   //
401   // lu_rcond (rcond from pre-computed LU decomposition)
402 
403   template<typename eT>
404   inline static eT lu_rcond(const Mat<eT>& A, const eT norm_val);
405 
406   template<typename  T>
407   inline static  T lu_rcond(const Mat< std::complex<T> >& A, const T norm_val);
408 
409   template<typename eT>
410   inline static eT lu_rcond_sympd(const Mat<eT>& A, const eT norm_val);
411 
412   template<typename  T>
413   inline static  T lu_rcond_sympd(const Mat< std::complex<T> >& A, const T norm_val);
414 
415   template<typename eT>
416   inline static eT lu_rcond_band(const Mat<eT>& AB, const uword KL, const uword KU, const podarray<blas_int>& ipiv, const eT norm_val);
417 
418   template<typename  T>
419   inline static  T lu_rcond_band(const Mat< std::complex<T> >& AB, const uword KL, const uword KU, const podarray<blas_int>& ipiv, const T norm_val);
420 
421 
422   //
423   // misc
424 
425   template<typename T1>
426   inline static bool crippled_lapack(const Base<typename T1::elem_type, T1>&);
427 
428   template<typename T1>
429   inline static typename T1::pod_type epsilon_lapack(const Base<typename T1::elem_type, T1>&);
430 
431   template<typename eT>
432   inline static bool rudimentary_sym_check(const Mat<eT>& X);
433 
434   template<typename T>
435   inline static bool rudimentary_sym_check(const Mat< std::complex<T> >& X);
436   };
437 
438 
439 
440 namespace qz_helper
441   {
442   template<typename T> inline blas_int select_lhp(const T* x_ptr, const T* y_ptr, const T* z_ptr);
443   template<typename T> inline blas_int select_rhp(const T* x_ptr, const T* y_ptr, const T* z_ptr);
444   template<typename T> inline blas_int select_iuc(const T* x_ptr, const T* y_ptr, const T* z_ptr);
445   template<typename T> inline blas_int select_ouc(const T* x_ptr, const T* y_ptr, const T* z_ptr);
446 
447   template<typename T> inline blas_int cx_select_lhp(const std::complex<T>* x_ptr, const std::complex<T>* y_ptr);
448   template<typename T> inline blas_int cx_select_rhp(const std::complex<T>* x_ptr, const std::complex<T>* y_ptr);
449   template<typename T> inline blas_int cx_select_iuc(const std::complex<T>* x_ptr, const std::complex<T>* y_ptr);
450   template<typename T> inline blas_int cx_select_ouc(const std::complex<T>* x_ptr, const std::complex<T>* y_ptr);
451 
452   template<typename T> inline void_ptr ptr_cast(blas_int (*function)(const T*, const T*, const T*));
453   template<typename T> inline void_ptr ptr_cast(blas_int (*function)(const std::complex<T>*, const std::complex<T>*));
454   }
455 
456 
457 
458 //! @}
459