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