From dd5e73ae200afb58bfae982a67211145f6f5e8ec Mon Sep 17 00:00:00 2001 From: conrad Date: Mon, 11 Nov 2024 14:14:20 +1000 Subject: [PATCH] add solve_sym_rcond --- include/armadillo_bits/auxlib_bones.hpp | 6 + include/armadillo_bits/auxlib_meat.hpp | 280 +++++++++++++++++++----- 2 files changed, 237 insertions(+), 49 deletions(-) diff --git a/include/armadillo_bits/auxlib_bones.hpp b/include/armadillo_bits/auxlib_bones.hpp index c05cda88..f68ba1f5 100644 --- a/include/armadillo_bits/auxlib_bones.hpp +++ b/include/armadillo_bits/auxlib_bones.hpp @@ -287,6 +287,12 @@ class auxlib template inline static bool solve_sym_fast(Mat< std::complex >& out, Mat< std::complex >& A, const Base< std::complex, T1 >& B_expr); + template + inline static bool solve_sym_rcond(Mat& out, typename T1::pod_type& out_rcond, Mat& A, const Base& B_expr); + + template + inline static bool solve_sym_rcond(Mat< std::complex >& out, typename T1::pod_type& out_rcond, Mat< std::complex >& A, const Base< std::complex,T1>& B_expr); + // template diff --git a/include/armadillo_bits/auxlib_meat.hpp b/include/armadillo_bits/auxlib_meat.hpp index 632b32c0..173efb2d 100644 --- a/include/armadillo_bits/auxlib_meat.hpp +++ b/include/armadillo_bits/auxlib_meat.hpp @@ -322,8 +322,9 @@ auxlib::inv_sym(Mat< std::complex >& A) #if defined(ARMA_CRIPPLED_LAPACK) { - arma_ignore(A); - return false; + arma_debug_print("auxlib::inv_sym(): redirecting to auxlib::inv() due to crippled LAPACK"); + + return auxlib::inv(A); } #elif defined(ARMA_USE_LAPACK) { @@ -408,17 +409,20 @@ auxlib::inv_sym_rcond(Mat& A, eT& out_rcond) podarray ipiv(A.n_rows); podarray iwork(A.n_rows); - eT work_query[2] = {}; - blas_int lwork_query = -1; - - arma_debug_print("lapack::sytrf()"); - lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); - - if(info != 0) { return false; } - - blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); - - lwork = (std::max)(lwork_proposed, lwork); + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::sytrf()"); + lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return false; } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } podarray work( static_cast(lwork) ); @@ -475,9 +479,9 @@ auxlib::inv_sym_rcond(Mat< std::complex >& A, T& out_rcond) #if defined(ARMA_CRIPPLED_LAPACK) { - arma_ignore(A); - arma_ignore(out_rcond); - return false; + arma_debug_print("auxlib::inv_sym_rcond(): redirecting to auxlib::inv_rcond() due to crippled LAPACK"); + + return auxlib::inv_rcond(A, out_rcond); } #elif defined(ARMA_USE_LAPACK) { @@ -497,17 +501,20 @@ auxlib::inv_sym_rcond(Mat< std::complex >& A, T& out_rcond) podarray ipiv(A.n_rows); podarray lanhe_work(A.n_rows); - eT work_query[2] = {}; - blas_int lwork_query = -1; - - arma_debug_print("lapack::hetrf()"); - lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); - - if(info != 0) { return false; } - - blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); - - lwork = (std::max)(lwork_proposed, lwork); + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::hetrf()"); + lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return false; } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } podarray work( static_cast(lwork) ); @@ -4697,6 +4704,175 @@ auxlib::solve_sym_fast(Mat< std::complex >& out, Mat< std +template +inline +bool +auxlib::solve_sym_rcond(Mat& out, typename T1::pod_type& out_rcond, Mat& A, const Base& B_expr) + { + arma_debug_sigprint(); + + typedef typename T1::pod_type eT; + + out = B_expr.get_ref(); + out_rcond = eT(0); + + const uword B_n_rows = out.n_rows; + const uword B_n_cols = out.n_cols; + + arma_conform_check( (A.n_rows != B_n_rows), "solve(): number of rows in given matrices must be the same", [&](){ out.soft_reset(); } ); + + if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return true; } + + #if defined(ARMA_USE_LAPACK) + { + arma_conform_assert_blas_size(A); + + char norm_id = '1'; + char uplo = 'L'; + blas_int n = blas_int(A.n_rows); + blas_int lda = blas_int(A.n_rows); + blas_int ldb = blas_int(out.n_rows); + blas_int nrhs = blas_int(out.n_cols); + blas_int lwork = (std::max)(blas_int(podarray_prealloc_n_elem::val), 2*n); // 2*n due to lapack::sycon() requirements + blas_int info = 0; + eT norm_val = eT(0); + eT tmp_rcond = eT(0); + + podarray ipiv(A.n_rows); + podarray iwork(A.n_rows); + + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::sytrf()"); + lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return false; } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } + + podarray work( static_cast(lwork) ); + + arma_debug_print("lapack::lansy()"); + norm_val = (has_blas_float_bug::value) ? auxlib::norm1_sym(A) : lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr()); + + arma_debug_print("lapack::sytrf()"); + lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), work.memptr(), &lwork, &info); + + if(info != 0) { return false; } + + arma_debug_print("lapack::sytrs()"); + lapack::sytrs(&uplo, &n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info); + + if(info != 0) { return false; } + + arma_debug_print("lapack::sycon()"); + lapack::sycon(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &norm_val, &tmp_rcond, work.memptr(), iwork.memptr(), &info); + + out_rcond = tmp_rcond; + + return (info == 0); + } + #else + { + arma_stop_logic_error("solve(): use of LAPACK must be enabled"); + return false; + } + #endif + } + + + +template +inline +bool +auxlib::solve_sym_rcond(Mat< std::complex >& out, typename T1::pod_type& out_rcond, Mat< std::complex >& A, const Base< std::complex,T1>& B_expr) + { + arma_debug_sigprint(); + + typedef typename T1::pod_type T; + typedef typename std::complex eT; + + out = B_expr.get_ref(); + out_rcond = T(0); + + const uword B_n_rows = out.n_rows; + const uword B_n_cols = out.n_cols; + + arma_conform_check( (A.n_rows != B_n_rows), "solve(): number of rows in given matrices must be the same", [&](){ out.soft_reset(); } ); + + if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return true; } + + #if defined(ARMA_USE_LAPACK) + { + arma_conform_assert_blas_size(A); + + char norm_id = '1'; + char uplo = 'L'; + blas_int n = blas_int(A.n_rows); + blas_int lda = blas_int(A.n_rows); + blas_int ldb = blas_int(out.n_rows); + blas_int nrhs = blas_int(out.n_cols); + blas_int lwork = (std::max)(blas_int(podarray_prealloc_n_elem::val), 2*n); // 2*n due to lapack::hecon() requirements + blas_int info = 0; + T norm_val = T(0); + T tmp_rcond = T(0); + + podarray ipiv(A.n_rows); + podarray lanhe_work(A.n_rows); + + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::hetrf()"); + lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return false; } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } + + podarray work( static_cast(lwork) ); + + arma_debug_print("lapack::lanhe()"); + norm_val = (has_blas_float_bug::value) ? auxlib::norm1_sym(A) : lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &lda, lanhe_work.memptr()); + + arma_debug_print("lapack::hetrf()"); + lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), work.memptr(), &lwork, &info); + + if(info != 0) { return false; } + + arma_debug_print("lapack::hetrs()"); + lapack::hetrs(&uplo, &n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info); + + if(info != 0) { return false; } + + arma_debug_print("lapack::hecon()"); + lapack::hecon(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &norm_val, &tmp_rcond, work.memptr(), &info); + + out_rcond = tmp_rcond; + + return (info == 0); + } + #else + { + arma_stop_logic_error("solve(): use of LAPACK must be enabled"); + return false; + } + #endif + } + + + template inline bool @@ -6642,17 +6818,20 @@ auxlib::rcond_sym(Mat& A) podarray ipiv(A.n_rows); podarray iwork(A.n_rows); - eT work_query[2] = {}; - blas_int lwork_query = -1; - - arma_debug_print("lapack::sytrf()"); - lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); - - if(info != 0) { return eT(0); } - - blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); - - lwork = (std::max)(lwork_proposed, lwork); + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::sytrf()"); + lapack::sytrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return eT(0); } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } podarray work( static_cast(lwork) ); @@ -6713,17 +6892,20 @@ auxlib::rcond_sym(Mat< std::complex >& A) podarray ipiv(A.n_rows); podarray lanhe_work(A.n_rows); - eT work_query[2] = {}; - blas_int lwork_query = -1; - - arma_debug_print("lapack::hetrf()"); - lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); - - if(info != 0) { return T(0); } - - blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); - - lwork = (std::max)(lwork_proposed, lwork); + if( (2*n) > blas_int(podarray_prealloc_n_elem::val) ) + { + eT work_query[2] = {}; + blas_int lwork_query = -1; + + arma_debug_print("lapack::hetrf()"); + lapack::hetrf(&uplo, &n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_query, &info); + + if(info != 0) { return T(0); } + + blas_int lwork_proposed = static_cast( access::tmp_real(work_query[0]) ); + + lwork = (std::max)(lwork_proposed, lwork); + } podarray work( static_cast(lwork) );