is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use.
| auxlib_meat.hpp (armadillo-10.8.2.tar.xz) | : | auxlib_meat.hpp (armadillo-11.0.0.tar.xz) |
|---|
| | |
| skipping to change at line 30 | | skipping to change at line 30 |
|---|
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::inv(Mat<eT>& A) | | auxlib::inv(Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| if(A.is_empty()) { return true; } | | if(A.is_empty()) { return true; } |
| | |
|
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_LAPACK) |
| { | | |
| arma_debug_assert_atlas_size(A); | | |
| | |
| podarray<int> ipiv(A.n_rows); | | |
| | |
| int info = 0; | | |
| | |
| arma_extra_debug_print("atlas::clapack_getrf()"); | | |
| info = atlas::clapack_getrf(atlas::CblasColMajor, A.n_rows, A.n_cols, A.memp | | |
| tr(), A.n_rows, ipiv.memptr()); | | |
| | |
| if(info != 0) { return false; } | | |
| | |
| arma_extra_debug_print("atlas::clapack_getri()"); | | |
| info = atlas::clapack_getri(atlas::CblasColMajor, A.n_rows, A.memptr(), A.n_ | | |
| rows, ipiv.memptr()); | | |
| | |
| return (info == 0); | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int lwork = (std::max)(blas_int(podarray_prealloc_n_elem::val), n); | | blas_int lwork = (std::max)(blas_int(podarray_prealloc_n_elem::val), n); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<blas_int> ipiv(A.n_rows); | | podarray<blas_int> ipiv(A.n_rows); |
| | |
|
| | arma_extra_debug_print("lapack::getrf()"); |
| | lapack::getrf(&n, &n, A.memptr(), &lda, ipiv.memptr(), &info); |
| | |
| | if(info != 0) { return false; } |
| | |
| if(n > 16) | | if(n > 16) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::getri()"); | | arma_extra_debug_print("lapack::getri()"); |
| lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_
query, &info); | | lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_
query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_que
ry[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_que
ry[0]) ); |
| | |
| lwork = (std::max)(lwork_proposed, lwork); | | lwork = (std::max)(lwork_proposed, lwork); |
| } | | } |
| | |
| podarray<eT> work( static_cast<uword>(lwork) ); | | podarray<eT> work( static_cast<uword>(lwork) ); |
| | |
|
| arma_extra_debug_print("lapack::getrf()"); | | |
| lapack::getrf(&n, &n, A.memptr(), &lda, ipiv.memptr(), &info); | | |
| | |
| if(info != 0) { return false; } | | |
| | | |
| arma_extra_debug_print("lapack::getri()"); | | arma_extra_debug_print("lapack::getri()"); |
| lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), work.memptr(), &lwork, &i
nfo); | | lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), work.memptr(), &lwork, &i
nfo); |
| | |
| return (info == 0); | | return (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
|
| arma_stop_logic_error("inv(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("inv(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::inv(Mat<eT>& out, const Mat<eT>& X) | | auxlib::inv(Mat<eT>& out, const Mat<eT>& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| out = X; | | out = X; |
| | |
| return auxlib::inv(out); | | return auxlib::inv(out); |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
|
| | auxlib::inv_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond) |
| | { |
| | arma_extra_debug_sigprint(); |
| | |
| | typedef typename get_pod_type<eT>::result T; |
| | |
| | out_rcond = T(0); |
| | |
| | if(A.is_empty()) { return true; } |
| | |
| | #if defined(ARMA_USE_LAPACK) |
| | { |
| | arma_debug_assert_blas_size(A); |
| | |
| | char norm_id = '1'; |
| | blas_int n = blas_int(A.n_rows); |
| | blas_int lda = blas_int(A.n_rows); |
| | blas_int lwork = (std::max)(blas_int(podarray_prealloc_n_elem::val), n); |
| | blas_int info = 0; |
| | T norm_val = T(0); |
| | |
| | podarray<T> junk(1); |
| | podarray<blas_int> ipiv(A.n_rows); |
| | |
| | arma_extra_debug_print("lapack::lange()"); |
| | norm_val = lapack::lange<eT>(&norm_id, &n, &n, A.memptr(), &lda, junk.memptr |
| | ()); |
| | |
| | arma_extra_debug_print("lapack::getrf()"); |
| | lapack::getrf(&n, &n, A.memptr(), &lda, ipiv.memptr(), &info); |
| | |
| | if(info != 0) { return false; } |
| | |
| | out_rcond = auxlib::lu_rcond<T>(A, norm_val); |
| | |
| | if(n > 16) |
| | { |
| | eT work_query[2] = {}; |
| | blas_int lwork_query = -1; |
| | |
| | arma_extra_debug_print("lapack::getri()"); |
| | lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), &work_query[0], &lwork_ |
| | query, &info); |
| | |
| | if(info != 0) { return false; } |
| | |
| | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_que |
| | ry[0]) ); |
| | |
| | lwork = (std::max)(lwork_proposed, lwork); |
| | } |
| | |
| | podarray<eT> work( static_cast<uword>(lwork) ); |
| | |
| | arma_extra_debug_print("lapack::getri()"); |
| | lapack::getri(&n, A.memptr(), &lda, ipiv.memptr(), work.memptr(), &lwork, &i |
| | nfo); |
| | |
| | return (info == 0); |
| | } |
| | #else |
| | { |
| | arma_ignore(A); |
| | arma_stop_logic_error("inv_rcond(): use of LAPACK must be enabled"); |
| | return false; |
| | } |
| | #endif |
| | } |
| | |
| | template<typename eT> |
| | inline |
| | bool |
| auxlib::inv_tr(Mat<eT>& A, const uword layout) | | auxlib::inv_tr(Mat<eT>& A, const uword layout) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| if(A.is_empty()) { return true; } | | if(A.is_empty()) { return true; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; |
| char diag = 'N'; | | char diag = 'N'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| arma_extra_debug_print("lapack::trtri()"); | | arma_extra_debug_print("lapack::trtri()"); |
| lapack::trtri(&uplo, &diag, &n, A.memptr(), &n, &info); | | lapack::trtri(&uplo, &diag, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
|
| if(layout == 0) | | |
| { | | |
| A = trimatu(A); // upper triangular | | |
| } | | |
| else | | |
| { | | |
| A = trimatl(A); // lower triangular | | |
| } | | |
| | | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(layout); | | arma_ignore(layout); |
| arma_stop_logic_error("inv(): use of LAPACK must be enabled"); | | arma_stop_logic_error("inv(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
|
| auxlib::inv_sympd(Mat<eT>& A) | | auxlib::inv_tr_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond, c
onst uword layout) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
|
| if(A.is_empty()) { return true; } | | #if defined(ARMA_USE_LAPACK) |
| | | |
| #if defined(ARMA_USE_ATLAS) | | |
| { | | { |
|
| arma_debug_assert_atlas_size(A); | | typedef typename get_pod_type<eT>::result T; |
| | |
|
| int info = 0; | | if(A.is_empty()) { return true; } |
| | |
|
| arma_extra_debug_print("atlas::clapack_potrf()"); | | out_rcond = auxlib::rcond_trimat(A, layout); |
| info = atlas::clapack_potrf(atlas::CblasColMajor, atlas::CblasLower, A.n_row | | |
| s, A.memptr(), A.n_rows); | | |
| | |
|
| if(info != 0) { return false; } | | arma_debug_assert_blas_size(A); |
| | |
|
| arma_extra_debug_print("atlas::clapack_potri()"); | | char uplo = (layout == 0) ? 'U' : 'L'; |
| info = atlas::clapack_potri(atlas::CblasColMajor, atlas::CblasLower, A.n_row | | char diag = 'N'; |
| s, A.memptr(), A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| | blas_int info = 0; |
| | |
|
| if(info != 0) { return false; } | | arma_extra_debug_print("lapack::trtri()"); |
| | lapack::trtri(&uplo, &diag, &n, A.memptr(), &n, &info); |
| | |
|
| A = symmatl(A); | | if(info != 0) { out_rcond = T(0); return false; } |
| | |
| return true; | | return true; |
| } | | } |
|
| #elif defined(ARMA_USE_LAPACK) | | #else |
| | { |
| | arma_ignore(A); |
| | arma_ignore(out_rcond); |
| | arma_ignore(layout); |
| | arma_stop_logic_error("inv(): use of LAPACK must be enabled"); |
| | return false; |
| | } |
| | #endif |
| | } |
| | |
| | template<typename eT> |
| | inline |
| | bool |
| | auxlib::inv_sympd(Mat<eT>& A, bool& out_sympd_state) |
| | { |
| | arma_extra_debug_sigprint(); |
| | |
| | out_sympd_state = false; |
| | |
| | if(A.is_empty()) { return true; } |
| | |
| | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| // NOTE: for complex matrices, zpotrf() assumes the matrix is hermitian (not
simply symmetric) | | // NOTE: for complex matrices, zpotrf() assumes the matrix is hermitian (not
simply symmetric) |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
|
| | out_sympd_state = true; |
| | |
| arma_extra_debug_print("lapack::potri()"); | | arma_extra_debug_print("lapack::potri()"); |
| lapack::potri(&uplo, &n, A.memptr(), &n, &info); | | lapack::potri(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| A = symmatl(A); | | A = symmatl(A); |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
|
| arma_stop_logic_error("inv_sympd(): use of ATLAS or LAPACK must be enabled") | | arma_ignore(out_sympd_state); |
| ; | | arma_stop_logic_error("inv_sympd(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::inv_sympd(Mat<eT>& out, const Mat<eT>& X) | | auxlib::inv_sympd(Mat<eT>& out, const Mat<eT>& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| out = X; | | out = X; |
| | |
|
| return auxlib::inv_sympd(out); | | bool sympd_state_junk = false; |
| | |
| | return auxlib::inv_sympd(out, sympd_state_junk); |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
|
| auxlib::inv_sympd_rcond(Mat<eT>& A, const eT rcond_threshold) | | auxlib::inv_sympd_rcond(Mat<eT>& A, bool& out_sympd_state, eT& out_rcond, const
eT rcond_threshold) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
|
| | out_sympd_state = false; |
| | |
| if(A.is_empty()) { return true; } | | if(A.is_empty()) { return true; } |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename get_pod_type<eT>::result T; | | typedef typename get_pod_type<eT>::result T; |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char norm_id = '1'; | | char norm_id = '1'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| | |
| skipping to change at line 254 | | skipping to change at line 324 |
|---|
| T norm_val = T(0); | | T norm_val = T(0); |
| | |
| podarray<T> work(A.n_rows); | | podarray<T> work(A.n_rows); |
| | |
| arma_extra_debug_print("lapack::lansy()"); | | arma_extra_debug_print("lapack::lansy()"); |
| norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; | | norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, A.memptr(), &n, &info); |
| | |
|
| if(info != 0) { return false; } | | if(info != 0) { out_rcond = eT(0); return false; } |
| | |
|
| const T rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); | | out_sympd_state = true; |
| | |
|
| if(rcond < rcond_threshold) { return false; } | | out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); |
| | |
| | if( (rcond_threshold > eT(0)) && (out_rcond < rcond_threshold) ) { return f |
| | alse; } |
| | |
| arma_extra_debug_print("lapack::potri()"); | | arma_extra_debug_print("lapack::potri()"); |
| lapack::potri(&uplo, &n, A.memptr(), &n, &info); | | lapack::potri(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| A = symmatl(A); | | A = symmatl(A); |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
|
| | arma_ignore(out_sympd_state); |
| | arma_ignore(out_rcond); |
| arma_ignore(rcond_threshold); | | arma_ignore(rcond_threshold); |
| arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); | | arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename T> | | template<typename T> |
| inline | | inline |
| bool | | bool |
|
| auxlib::inv_sympd_rcond(Mat< std::complex<T> >& A, const T rcond_threshold) | | auxlib::inv_sympd_rcond(Mat< std::complex<T> >& A, bool& out_sympd_state, T& out
_rcond, const T rcond_threshold) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
|
| | out_sympd_state = false; |
| | |
| if(A.is_empty()) { return true; } | | if(A.is_empty()) { return true; } |
| | |
| #if defined(ARMA_CRIPPLED_LAPACK) | | #if defined(ARMA_CRIPPLED_LAPACK) |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
|
| | arma_ignore(out_sympd_state); |
| | arma_ignore(out_rcond); |
| arma_ignore(rcond_threshold); | | arma_ignore(rcond_threshold); |
| return false; | | return false; |
| } | | } |
| #elif defined(ARMA_USE_LAPACK) | | #elif defined(ARMA_USE_LAPACK) |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char norm_id = '1'; | | char norm_id = '1'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| | |
| skipping to change at line 312 | | skipping to change at line 390 |
|---|
| T norm_val = T(0); | | T norm_val = T(0); |
| | |
| podarray<T> work(A.n_rows); | | podarray<T> work(A.n_rows); |
| | |
| arma_extra_debug_print("lapack::lanhe()"); | | arma_extra_debug_print("lapack::lanhe()"); |
| norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; | | norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, A.memptr(), &n, &info); |
| | |
|
| if(info != 0) { return false; } | | if(info != 0) { out_rcond = T(0); return false; } |
| | |
| | out_sympd_state = true; |
| | |
|
| const T rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); | | out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); |
| | |
|
| if(rcond < rcond_threshold) { return false; } | | if( (rcond_threshold > T(0)) && (out_rcond < rcond_threshold) ) { return fa
lse; } |
| | |
| arma_extra_debug_print("lapack::potri()"); | | arma_extra_debug_print("lapack::potri()"); |
| lapack::potri(&uplo, &n, A.memptr(), &n, &info); | | lapack::potri(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| A = symmatl(A); | | A = symmatl(A); |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
|
| | arma_ignore(out_sympd_state); |
| | arma_ignore(out_rcond); |
| arma_ignore(rcond_threshold); | | arma_ignore(rcond_threshold); |
| arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); | | arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! determinant of a matrix | | //! determinant of a matrix |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::det(eT& out_val, Mat<eT>& A) | | auxlib::det(eT& out_val, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| if(A.is_empty()) { out_val = eT(1); return true; } | | if(A.is_empty()) { out_val = eT(1); return true; } |
| | |
|
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_LAPACK) |
| { | | |
| arma_debug_assert_atlas_size(A); | | |
| | |
| podarray<int> ipiv(A.n_rows); | | |
| | |
| arma_extra_debug_print("atlas::clapack_getrf()"); | | |
| const int info = atlas::clapack_getrf(atlas::CblasColMajor, A.n_rows, A.n_co | | |
| ls, A.memptr(), A.n_rows, ipiv.memptr()); | | |
| | |
| if(info < 0) { return false; } | | |
| | |
| // on output A appears to be L+U_alt, where U_alt is U with the main diagona | | |
| l set to zero | | |
| eT val = A.at(0,0); | | |
| for(uword i=1; i < A.n_rows; ++i) { val *= A.at(i,i); } | | |
| | |
| int sign = +1; | | |
| for(uword i=0; i < A.n_rows; ++i) | | |
| { | | |
| // NOTE: no adjustment required, as the clapack version of getrf() assumes | | |
| counting from 0 | | |
| if( int(i) != ipiv.mem[i] ) { sign *= -1; } | | |
| } | | |
| | |
| out_val = (sign < 0) ? eT(-val) : eT(val); | | |
| | |
| return true; | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| podarray<blas_int> ipiv(A.n_rows); | | podarray<blas_int> ipiv(A.n_rows); |
| | |
| blas_int info = 0; | | blas_int info = 0; |
| blas_int n_rows = blas_int(A.n_rows); | | blas_int n_rows = blas_int(A.n_rows); |
| blas_int n_cols = blas_int(A.n_cols); | | blas_int n_cols = blas_int(A.n_cols); |
| | |
| arma_extra_debug_print("lapack::getrf()"); | | arma_extra_debug_print("lapack::getrf()"); |
| | |
| skipping to change at line 407 | | skipping to change at line 463 |
|---|
| } | | } |
| | |
| out_val = (sign < 0) ? eT(-val) : eT(val); | | out_val = (sign < 0) ? eT(-val) : eT(val); |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out_val); | | arma_ignore(out_val); |
| arma_ignore(A); | | arma_ignore(A); |
|
| arma_stop_logic_error("det(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("det(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! log determinant of a matrix | | //! log determinant of a matrix |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, Mat<eT
>& A) | | auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, Mat<eT
>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| typedef typename get_pod_type<eT>::result T; | | typedef typename get_pod_type<eT>::result T; |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { out_val = eT(0); out_sign = T(1); return true; } |
| { | | |
| out_val = eT(0); | | |
| out_sign = T(1); | | |
| return true; | | |
| } | | |
| | | |
| #if defined(ARMA_USE_ATLAS) | | |
| { | | |
| arma_debug_assert_atlas_size(A); | | |
| | |
| podarray<int> ipiv(A.n_rows); | | |
| | |
| arma_extra_debug_print("atlas::clapack_getrf()"); | | |
| const int info = atlas::clapack_getrf(atlas::CblasColMajor, A.n_rows, A.n_co | | |
| ls, A.memptr(), A.n_rows, ipiv.memptr()); | | |
| | |
| if(info < 0) { return false; } | | |
| | |
| // on output A appears to be L+U_alt, where U_alt is U with the main diagona | | |
| l set to zero | | |
| | |
| sword sign = (is_cx<eT>::no) ? ( (access::tmp_real( A.at(0,0) ) < T(0)) ? -1 | | |
| : +1 ) : +1; | | |
| eT val = (is_cx<eT>::no) ? std::log( (access::tmp_real( A.at(0,0) ) < T( | | |
| 0)) ? A.at(0,0)*T(-1) : A.at(0,0) ) : std::log( A.at(0,0) ); | | |
| | |
| for(uword i=1; i < A.n_rows; ++i) | | |
| { | | |
| const eT x = A.at(i,i); | | |
| | |
| sign *= (is_cx<eT>::no) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1; | | |
| val += (is_cx<eT>::no) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) | | |
| : x ) : std::log(x); | | |
| } | | |
| | |
| for(uword i=0; i < A.n_rows; ++i) | | |
| { | | |
| if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the clapa | | |
| ck version of getrf() assumes counting from 0 | | |
| { | | |
| sign *= -1; | | |
| } | | |
| } | | |
| | |
| out_val = val; | | |
| out_sign = T(sign); | | |
| | |
|
| return true; | | #if defined(ARMA_USE_LAPACK) |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| podarray<blas_int> ipiv(A.n_rows); | | podarray<blas_int> ipiv(A.n_rows); |
| | |
| blas_int info = 0; | | blas_int info = 0; |
| blas_int n_rows = blas_int(A.n_rows); | | blas_int n_rows = blas_int(A.n_rows); |
| blas_int n_cols = blas_int(A.n_cols); | | blas_int n_cols = blas_int(A.n_cols); |
| | |
| arma_extra_debug_print("lapack::getrf()"); | | arma_extra_debug_print("lapack::getrf()"); |
| | |
| skipping to change at line 513 | | skipping to change at line 527 |
|---|
| out_val = val; | | out_val = val; |
| out_sign = T(sign); | | out_sign = T(sign); |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(out_val); | | arma_ignore(out_val); |
| arma_ignore(out_sign); | | arma_ignore(out_sign); |
|
| arma_stop_logic_error("log_det(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("log_det(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::log_det_sympd(typename get_pod_type<eT>::result& out_val, Mat<eT>& A) | | auxlib::log_det_sympd(typename get_pod_type<eT>::result& out_val, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| typedef typename get_pod_type<eT>::result T; | | typedef typename get_pod_type<eT>::result T; |
| | |
| if(A.is_empty()) { out_val = T(0); return true; } | | if(A.is_empty()) { out_val = T(0); return true; } |
| | |
|
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_LAPACK) |
| { | | |
| arma_debug_assert_atlas_size(A); | | |
| | |
| int info = 0; | | |
| | |
| arma_extra_debug_print("atlas::clapack_potrf()"); | | |
| info = atlas::clapack_potrf(atlas::CblasColMajor, atlas::CblasLower, A.n_row | | |
| s, A.memptr(), A.n_rows); | | |
| | |
| if(info != 0) { return false; } | | |
| | |
| T val = std::log( access::tmp_real(A.at(0,0)) ); | | |
| | |
| for(uword i=1; i < A.n_rows; ++i) { val += std::log( access::tmp_real(A.at( | | |
| i,i)) ); } | | |
| | |
| out_val = T(2) * val; | | |
| | |
| return true; | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
|
| T val = std::log( access::tmp_real(A.at(0,0)) ); | | T val = T(0); |
| | |
|
| for(uword i=1; i < A.n_rows; ++i) { val += std::log( access::tmp_real(A.at(
i,i)) ); } | | for(uword i=0; i < A.n_rows; ++i) { val += std::log( access::tmp_real(A.at(
i,i)) ); } |
| | |
| out_val = T(2) * val; | | out_val = T(2) * val; |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out_val); | | arma_ignore(out_val); |
| arma_ignore(A); | | arma_ignore(A); |
|
| arma_stop_logic_error("det(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("log_det_sympd(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! LU decomposition of a matrix | | //! LU decomposition of a matrix |
| template<typename eT, typename T1> | | template<typename eT, typename T1> |
| inline | | inline |
| bool | | bool |
| auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>&
X) | | auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>&
X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| U = X.get_ref(); | | U = X.get_ref(); |
| | |
| const uword U_n_rows = U.n_rows; | | const uword U_n_rows = U.n_rows; |
| const uword U_n_cols = U.n_cols; | | const uword U_n_cols = U.n_cols; |
| | |
|
| if(U.is_empty()) | | if(U.is_empty()) { L.set_size(U_n_rows, 0); U.set_size(0, U_n_cols); ipiv.res |
| { | | et(); return true; } |
| L.set_size(U_n_rows, 0); | | |
| U.set_size(0, U_n_cols); | | |
| ipiv.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
|
| #if defined(ARMA_USE_ATLAS) | | arma_debug_assert_blas_size(U); |
| { | | |
| arma_debug_assert_atlas_size(U); | | |
| | |
|
| ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | | ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); |
| | |
|
| arma_extra_debug_print("atlas::clapack_getrf()"); | | blas_int info = 0; |
| int info = atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_cols, | | |
| U.memptr(), U_n_rows, ipiv.memptr()); | | |
| | |
| if(info < 0) { return false; } | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | |
| arma_debug_assert_blas_size(U); | | |
| | |
| ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | | |
| | | |
| blas_int info = 0; | | |
| | |
|
| blas_int n_rows = blas_int(U_n_rows); | | blas_int n_rows = blas_int(U_n_rows); |
| blas_int n_cols = blas_int(U_n_cols); | | blas_int n_cols = blas_int(U_n_cols); |
| | |
|
| arma_extra_debug_print("lapack::getrf()"); | | arma_extra_debug_print("lapack::getrf()"); |
| lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info) | | lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info); |
| ; | | |
| | |
|
| if(info < 0) { return false; } | | if(info < 0) { return false; } |
| | |
|
| // take into account that Fortran counts from 1 | | // take into account that Fortran counts from 1 |
| arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem); | | arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem); |
| } | | |
| #endif | | |
| | |
| L.copy_size(U); | | L.copy_size(U); |
| | |
| for(uword col=0; col < U_n_cols; ++col) | | for(uword col=0; col < U_n_cols; ++col) |
| { | | { |
| for(uword row=0; (row < col) && (row < U_n_rows); ++row) | | for(uword row=0; (row < col) && (row < U_n_rows); ++row) |
| { | | { |
| L.at(row,col) = eT(0); | | L.at(row,col) = eT(0); |
| } | | } |
| | |
| | |
| skipping to change at line 660 | | skipping to change at line 634 |
|---|
| { | | { |
| L.at(row,col) = U.at(row,col); | | L.at(row,col) = U.at(row,col); |
| U.at(row,col) = eT(0); | | U.at(row,col) = eT(0); |
| } | | } |
| } | | } |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
|
| arma_stop_logic_error("lu(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("lu(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> |
| inline | | inline |
| bool | | bool |
| auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X) | | auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X) |
| { | | { |
| | |
| skipping to change at line 807 | | skipping to change at line 781 |
|---|
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
| Mat<T> X = expr.get_ref(); | | Mat<T> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| Mat<T> tmp(1, 1, arma_nozeros_indicator()); | | Mat<T> tmp(1, 1, arma_nozeros_indicator()); |
| | |
| if(vecs_on) | | if(vecs_on) |
| { | | { |
| vecs.set_size(X.n_rows, X.n_rows); | | vecs.set_size(X.n_rows, X.n_rows); |
| tmp.set_size(X.n_rows, X.n_rows); | | tmp.set_size(X.n_rows, X.n_rows); |
| } | | } |
| | |
| skipping to change at line 918 | | skipping to change at line 887 |
|---|
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| Mat<eT> X = expr.get_ref(); | | Mat<eT> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| if(vecs_on) { vecs.set_size(X.n_rows, X.n_rows); } | | if(vecs_on) { vecs.set_size(X.n_rows, X.n_rows); } |
| | |
| podarray<eT> junk(1); | | podarray<eT> junk(1); |
| | |
| char jobvl = 'N'; | | char jobvl = 'N'; |
| char jobvr = (vecs_on) ? 'V' : 'N'; | | char jobvr = (vecs_on) ? 'V' : 'N'; |
| blas_int N = blas_int(X.n_rows); | | blas_int N = blas_int(X.n_rows); |
| | |
| skipping to change at line 988 | | skipping to change at line 952 |
|---|
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
| Mat<T> X = expr.get_ref(); | | Mat<T> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| Mat<T> tmp(1, 1, arma_nozeros_indicator()); | | Mat<T> tmp(1, 1, arma_nozeros_indicator()); |
| | |
| if(vecs_on) | | if(vecs_on) |
| { | | { |
| vecs.set_size(X.n_rows, X.n_rows); | | vecs.set_size(X.n_rows, X.n_rows); |
| tmp.set_size(X.n_rows, X.n_rows); | | tmp.set_size(X.n_rows, X.n_rows); |
| } | | } |
| | |
| skipping to change at line 1115 | | skipping to change at line 1074 |
|---|
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| Mat<eT> X = expr.get_ref(); | | Mat<eT> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| if(vecs_on) { vecs.set_size(X.n_rows, X.n_rows); } | | if(vecs_on) { vecs.set_size(X.n_rows, X.n_rows); } |
| | |
| podarray<eT> junk(1); | | podarray<eT> junk(1); |
| | |
| char bal = 'B'; | | char bal = 'B'; |
| char jobvl = 'N'; | | char jobvl = 'N'; |
| char jobvr = (vecs_on) ? 'V' : 'N'; | | char jobvr = (vecs_on) ? 'V' : 'N'; |
| | |
| skipping to change at line 1194 | | skipping to change at line 1148 |
|---|
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
| Mat<T> X = expr.get_ref(); | | Mat<T> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| lvecs.set_size(X.n_rows, X.n_rows); | | lvecs.set_size(X.n_rows, X.n_rows); |
| rvecs.set_size(X.n_rows, X.n_rows); | | rvecs.set_size(X.n_rows, X.n_rows); |
| | |
| Mat<T> ltmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); | | Mat<T> ltmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); |
| Mat<T> rtmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); | | Mat<T> rtmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); |
| | |
| char jobvl = 'V'; | | char jobvl = 'V'; |
| | |
| skipping to change at line 1299 | | skipping to change at line 1247 |
|---|
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| Mat<eT> X = expr.get_ref(); | | Mat<eT> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| lvecs.set_size(X.n_rows, X.n_rows); | | lvecs.set_size(X.n_rows, X.n_rows); |
| rvecs.set_size(X.n_rows, X.n_rows); | | rvecs.set_size(X.n_rows, X.n_rows); |
| | |
| char jobvl = 'V'; | | char jobvl = 'V'; |
| char jobvr = 'V'; | | char jobvr = 'V'; |
| blas_int N = blas_int(X.n_rows); | | blas_int N = blas_int(X.n_rows); |
| blas_int ldvl = blas_int(lvecs.n_rows); | | blas_int ldvl = blas_int(lvecs.n_rows); |
| | |
| skipping to change at line 1367 | | skipping to change at line 1309 |
|---|
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
| Mat<T> X = expr.get_ref(); | | Mat<T> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| lvecs.set_size(X.n_rows, X.n_rows); | | lvecs.set_size(X.n_rows, X.n_rows); |
| rvecs.set_size(X.n_rows, X.n_rows); | | rvecs.set_size(X.n_rows, X.n_rows); |
| | |
| Mat<T> ltmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); | | Mat<T> ltmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); |
| Mat<T> rtmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); | | Mat<T> rtmp(X.n_rows, X.n_rows, arma_nozeros_indicator()); |
| | |
| char bal = 'B'; | | char bal = 'B'; |
| | |
| skipping to change at line 1488 | | skipping to change at line 1424 |
|---|
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| Mat<eT> X = expr.get_ref(); | | Mat<eT> X = expr.get_ref(); |
| | |
| arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_gen(): given matrix must be
square sized" ); |
| | |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
|
| if(X.is_empty()) | | if(X.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(X.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && X.has_nonfinite()) { return false; } |
| | |
| vals.set_size(X.n_rows, 1); | | vals.set_size(X.n_rows, 1); |
| | |
| lvecs.set_size(X.n_rows, X.n_rows); | | lvecs.set_size(X.n_rows, X.n_rows); |
| rvecs.set_size(X.n_rows, X.n_rows); | | rvecs.set_size(X.n_rows, X.n_rows); |
| | |
| char bal = 'B'; | | char bal = 'B'; |
| char jobvl = 'V'; | | char jobvl = 'V'; |
| char jobvr = 'V'; | | char jobvr = 'V'; |
| char sense = 'N'; | | char sense = 'N'; |
| | |
| skipping to change at line 1570 | | skipping to change at line 1500 |
|---|
| | |
| Mat<T> A(A_expr.get_ref()); | | Mat<T> A(A_expr.get_ref()); |
| Mat<T> B(B_expr.get_ref()); | | Mat<T> B(B_expr.get_ref()); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(A.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| if(B.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| vals.set_size(A.n_rows, 1); | | vals.set_size(A.n_rows, 1); |
| | |
| Mat<T> tmp(1, 1, arma_nozeros_indicator()); | | Mat<T> tmp(1, 1, arma_nozeros_indicator()); |
| | |
| if(vecs_on) | | if(vecs_on) |
| { | | { |
| vecs.set_size(A.n_rows, A.n_rows); | | vecs.set_size(A.n_rows, A.n_rows); |
| tmp.set_size(A.n_rows, A.n_rows); | | tmp.set_size(A.n_rows, A.n_rows); |
| } | | } |
| | |
| skipping to change at line 1711 | | skipping to change at line 1636 |
|---|
| | |
| Mat<eT> A(A_expr.get_ref()); | | Mat<eT> A(A_expr.get_ref()); |
| Mat<eT> B(B_expr.get_ref()); | | Mat<eT> B(B_expr.get_ref()); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { vals.reset(); vecs.reset(); return true; } |
| { | | |
| vals.reset(); | | |
| vecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(A.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| if(B.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| vals.set_size(A.n_rows, 1); | | vals.set_size(A.n_rows, 1); |
| | |
| if(vecs_on) { vecs.set_size(A.n_rows, A.n_rows); } | | if(vecs_on) { vecs.set_size(A.n_rows, A.n_rows); } |
| | |
| podarray<eT> junk(1); | | podarray<eT> junk(1); |
| | |
| char jobvl = 'N'; | | char jobvl = 'N'; |
| char jobvr = (vecs_on) ? 'V' : 'N'; | | char jobvr = (vecs_on) ? 'V' : 'N'; |
| blas_int N = blas_int(A.n_rows); | | blas_int N = blas_int(A.n_rows); |
| | |
| skipping to change at line 1811 | | skipping to change at line 1731 |
|---|
| | |
| Mat<T> A(A_expr.get_ref()); | | Mat<T> A(A_expr.get_ref()); |
| Mat<T> B(B_expr.get_ref()); | | Mat<T> B(B_expr.get_ref()); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(A.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| if(B.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| vals.set_size(A.n_rows, 1); | | vals.set_size(A.n_rows, 1); |
| | |
| lvecs.set_size(A.n_rows, A.n_rows); | | lvecs.set_size(A.n_rows, A.n_rows); |
| rvecs.set_size(A.n_rows, A.n_rows); | | rvecs.set_size(A.n_rows, A.n_rows); |
| | |
| Mat<T> ltmp(A.n_rows, A.n_rows, arma_nozeros_indicator()); | | Mat<T> ltmp(A.n_rows, A.n_rows, arma_nozeros_indicator()); |
| Mat<T> rtmp(A.n_rows, A.n_rows, arma_nozeros_indicator()); | | Mat<T> rtmp(A.n_rows, A.n_rows, arma_nozeros_indicator()); |
| | |
| char jobvl = 'V'; | | char jobvl = 'V'; |
| | |
| skipping to change at line 1946 | | skipping to change at line 1860 |
|---|
| | |
| Mat<eT> A(A_expr.get_ref()); | | Mat<eT> A(A_expr.get_ref()); |
| Mat<eT> B(B_expr.get_ref()); | | Mat<eT> B(B_expr.get_ref()); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "e
ig_pair(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "eig_pair(): given matrices must h
ave the same size" ); |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { vals.reset(); lvecs.reset(); rvecs.reset(); return true; |
| { | | } |
| vals.reset(); | | |
| lvecs.reset(); | | |
| rvecs.reset(); | | |
| return true; | | |
| } | | |
| | |
|
| if(A.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| if(B.is_finite() == false) { return false; } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| vals.set_size(A.n_rows, 1); | | vals.set_size(A.n_rows, 1); |
| | |
| lvecs.set_size(A.n_rows, A.n_rows); | | lvecs.set_size(A.n_rows, A.n_rows); |
| rvecs.set_size(A.n_rows, A.n_rows); | | rvecs.set_size(A.n_rows, A.n_rows); |
| | |
| char jobvl = 'V'; | | char jobvl = 'V'; |
| char jobvr = 'V'; | | char jobvr = 'V'; |
| blas_int N = blas_int(A.n_rows); | | blas_int N = blas_int(A.n_rows); |
| blas_int ldvl = blas_int(lvecs.n_rows); | | blas_int ldvl = blas_int(lvecs.n_rows); |
| | |
| skipping to change at line 2029 | | skipping to change at line 1937 |
|---|
| auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& A) | | auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
| if(A.is_empty()) { eigval.reset(); return true; } | | if(A.is_empty()) { eigval.reset(); return true; } |
| | |
|
| // if(auxlib::rudimentary_sym_check(A) == false) | | |
| // { | | |
| // arma_debug_warn_level(1, "eig_sym(): given matrix is not symmetric"); | | |
| // return false; | | |
| // } | | |
| | | |
| if((arma_config::debug) && (auxlib::rudimentary_sym_check(A) == false)) | | if((arma_config::debug) && (auxlib::rudimentary_sym_check(A) == false)) |
| { | | { |
| arma_debug_warn_level(1, "eig_sym(): given matrix is not symmetric"); | | arma_debug_warn_level(1, "eig_sym(): given matrix is not symmetric"); |
| } | | } |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(A)) { |
| | return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| eigval.set_size(A.n_rows); | | eigval.set_size(A.n_rows); |
| | |
| char jobz = 'N'; | | char jobz = 'N'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(A.n_rows); | | blas_int N = blas_int(A.n_rows); |
| blas_int lwork = (64+2)*N; // lwork_min = (std::max)(blas_int(1), 3*N-1) | | blas_int lwork = (64+2)*N; // lwork_min = (std::max)(blas_int(1), 3*N-1) |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| skipping to change at line 2084 | | skipping to change at line 1988 |
|---|
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
| if(A.is_empty()) { eigval.reset(); return true; } | | if(A.is_empty()) { eigval.reset(); return true; } |
| | |
|
| // if(auxlib::rudimentary_sym_check(A) == false) | | |
| // { | | |
| // arma_debug_warn_level(1, "eig_sym(): given matrix is not hermitian"); | | |
| // return false; | | |
| // } | | |
| | | |
| if((arma_config::debug) && (auxlib::rudimentary_sym_check(A) == false)) | | if((arma_config::debug) && (auxlib::rudimentary_sym_check(A) == false)) |
| { | | { |
| arma_debug_warn_level(1, "eig_sym(): given matrix is not hermitian"); | | arma_debug_warn_level(1, "eig_sym(): given matrix is not hermitian"); |
| } | | } |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(A)) { |
| | return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| eigval.set_size(A.n_rows); | | eigval.set_size(A.n_rows); |
| | |
| char jobz = 'N'; | | char jobz = 'N'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(A.n_rows); | | blas_int N = blas_int(A.n_rows); |
| blas_int lwork = (64+1)*N; // lwork_min = (std::max)(blas_int(1), 2*N-1) | | blas_int lwork = (64+1)*N; // lwork_min = (std::max)(blas_int(1), 2*N-1) |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| skipping to change at line 2136 | | skipping to change at line 2036 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X) | | auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(X)) { |
| | return false; } |
| | |
| eigvec = X; | | eigvec = X; |
| | |
|
| if(eigvec.is_empty()) | | if(eigvec.is_empty()) { eigval.reset(); eigvec.reset(); return true; } |
| { | | |
| eigval.reset(); | | |
| eigvec.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(eigvec); | | arma_debug_assert_blas_size(eigvec); |
| | |
| eigval.set_size(eigvec.n_rows); | | eigval.set_size(eigvec.n_rows); |
| | |
| char jobz = 'V'; | | char jobz = 'V'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(eigvec.n_rows); | | blas_int N = blas_int(eigvec.n_rows); |
| blas_int lwork = (64+2)*N; // lwork_min = (std::max)(blas_int(1), 3*N-1) | | blas_int lwork = (64+2)*N; // lwork_min = (std::max)(blas_int(1), 3*N-1) |
| | |
| skipping to change at line 2188 | | skipping to change at line 2085 |
|---|
| auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< std::
complex<T> >& X) | | auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< std::
complex<T> >& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(X)) { |
| | return false; } |
| | |
| eigvec = X; | | eigvec = X; |
| | |
|
| if(eigvec.is_empty()) | | if(eigvec.is_empty()) { eigval.reset(); eigvec.reset(); return true; } |
| { | | |
| eigval.reset(); | | |
| eigvec.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(eigvec); | | arma_debug_assert_blas_size(eigvec); |
| | |
| eigval.set_size(eigvec.n_rows); | | eigval.set_size(eigvec.n_rows); |
| | |
| char jobz = 'V'; | | char jobz = 'V'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(eigvec.n_rows); | | blas_int N = blas_int(eigvec.n_rows); |
| blas_int lwork = (64+1)*N; // lwork_min = (std::max)(blas_int(1), 2*N-1) | | blas_int lwork = (64+1)*N; // lwork_min = (std::max)(blas_int(1), 2*N-1) |
| | |
| skipping to change at line 2239 | | skipping to change at line 2133 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::eig_sym_dc(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X) | | auxlib::eig_sym_dc(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(X)) { |
| | return false; } |
| | |
| eigvec = X; | | eigvec = X; |
| | |
|
| if(eigvec.is_empty()) | | if(eigvec.is_empty()) { eigval.reset(); eigvec.reset(); return true; } |
| { | | |
| eigval.reset(); | | |
| eigvec.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(eigvec); | | arma_debug_assert_blas_size(eigvec); |
| | |
| eigval.set_size(eigvec.n_rows); | | eigval.set_size(eigvec.n_rows); |
| | |
| char jobz = 'V'; | | char jobz = 'V'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(eigvec.n_rows); | | blas_int N = blas_int(eigvec.n_rows); |
| blas_int lwork_min = 1 + 6*N + 2*(N*N); | | blas_int lwork_min = 1 + 6*N + 2*(N*N); |
| blas_int liwork_min = 3 + 5*N; | | blas_int liwork_min = 3 + 5*N; |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| blas_int liwork_proposed = 0; | | blas_int liwork_proposed = 0; |
| | |
| if(N >= 32) | | if(N >= 32) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int iwork_query[2]; | | blas_int iwork_query[2] = {}; |
| | |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| blas_int liwork_query = -1; | | blas_int liwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::syevd()"); | | arma_extra_debug_print("lapack::syevd()"); |
| lapack::syevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), &wor
k_query[0], &lwork_query, &iwork_query[0], &liwork_query, &info); | | lapack::syevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), &wor
k_query[0], &lwork_query, &iwork_query[0], &liwork_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( work_query[0] ); | | lwork_proposed = static_cast<blas_int>( work_query[0] ); |
| | |
| skipping to change at line 2316 | | skipping to change at line 2207 |
|---|
| auxlib::eig_sym_dc(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< st
d::complex<T> >& X) | | auxlib::eig_sym_dc(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< st
d::complex<T> >& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); | | arma_debug_check( (X.is_square() == false), "eig_sym(): given matrix must be
square sized" ); |
| | |
|
| | if(arma_config::check_nonfinite && trimat_helper::has_nonfinite_triu(X)) { |
| | return false; } |
| | |
| eigvec = X; | | eigvec = X; |
| | |
|
| if(eigvec.is_empty()) | | if(eigvec.is_empty()) { eigval.reset(); eigvec.reset(); return true; } |
| { | | |
| eigval.reset(); | | |
| eigvec.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(eigvec); | | arma_debug_assert_blas_size(eigvec); |
| | |
| eigval.set_size(eigvec.n_rows); | | eigval.set_size(eigvec.n_rows); |
| | |
| char jobz = 'V'; | | char jobz = 'V'; |
| char uplo = 'U'; | | char uplo = 'U'; |
| | |
| blas_int N = blas_int(eigvec.n_rows); | | blas_int N = blas_int(eigvec.n_rows); |
| blas_int lwork_min = 2*N + N*N; | | blas_int lwork_min = 2*N + N*N; |
| blas_int lrwork_min = 1 + 5*N + 2*(N*N); | | blas_int lrwork_min = 1 + 5*N + 2*(N*N); |
| blas_int liwork_min = 3 + 5*N; | | blas_int liwork_min = 3 + 5*N; |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| blas_int lrwork_proposed = 0; | | blas_int lrwork_proposed = 0; |
| blas_int liwork_proposed = 0; | | blas_int liwork_proposed = 0; |
| | |
| if(N >= 32) | | if(N >= 32) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| T rwork_query[2]; | | T rwork_query[2] = {}; |
| blas_int iwork_query[2]; | | blas_int iwork_query[2] = {}; |
| | |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| blas_int lrwork_query = -1; | | blas_int lrwork_query = -1; |
| blas_int liwork_query = -1; | | blas_int liwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::heevd()"); | | arma_extra_debug_print("lapack::heevd()"); |
| lapack::heevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), &wor
k_query[0], &lwork_query, &rwork_query[0], &lrwork_query, &iwork_query[0], &liwo
rk_query, &info); | | lapack::heevd(&jobz, &uplo, &N, eigvec.memptr(), &N, eigval.memptr(), &wor
k_query[0], &lwork_query, &rwork_query[0], &lrwork_query, &iwork_query[0], &liwo
rk_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| | |
| skipping to change at line 2393 | | skipping to change at line 2281 |
|---|
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::chol_simple(Mat<eT>& X) | | auxlib::chol_simple(Mat<eT>& X) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
|
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_LAPACK) |
| { | | |
| arma_debug_assert_atlas_size(X); | | |
| | |
| int info = 0; | | |
| | |
| arma_extra_debug_print("atlas::clapack_potrf()"); | | |
| info = atlas::clapack_potrf(atlas::CblasColMajor, atlas::CblasUpper, X.n_row | | |
| s, X.memptr(), X.n_rows); | | |
| | |
| return (info == 0); | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
| char uplo = 'U'; | | char uplo = 'U'; |
| blas_int n = blas_int(X.n_rows); | | blas_int n = blas_int(X.n_rows); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, X.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, X.memptr(), &n, &info); |
| | |
| return (info == 0); | | return (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(X); | | arma_ignore(X); |
| | |
|
| arma_stop_logic_error("chol(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("chol(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::chol(Mat<eT>& X, const uword layout) | | auxlib::chol(Mat<eT>& X, const uword layout) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
|
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_LAPACK) |
| { | | |
| arma_debug_assert_atlas_size(X); | | |
| | |
| int info = 0; | | |
| | |
| arma_extra_debug_print("atlas::clapack_potrf()"); | | |
| info = atlas::clapack_potrf(atlas::CblasColMajor, ((layout == 0) ? atlas::Cb | | |
| lasUpper : atlas::CblasLower), X.n_rows, X.memptr(), X.n_rows); | | |
| | |
| if(info != 0) { return false; } | | |
| | |
| X = (layout == 0) ? trimatu(X) : trimatl(X); // trimatu() and trimatl() ret | | |
| urn the same type | | |
| | |
| return true; | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(X); | | arma_debug_assert_blas_size(X); |
| | |
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; |
| blas_int n = blas_int(X.n_rows); | | blas_int n = blas_int(X.n_rows); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf(&uplo, &n, X.memptr(), &n, &info); | | lapack::potrf(&uplo, &n, X.memptr(), &n, &info); |
| | |
| | |
| skipping to change at line 2471 | | skipping to change at line 2333 |
|---|
| | |
| X = (layout == 0) ? trimatu(X) : trimatl(X); // trimatu() and trimatl() ret
urn the same type | | X = (layout == 0) ? trimatu(X) : trimatl(X); // trimatu() and trimatl() ret
urn the same type |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(X); | | arma_ignore(X); |
| arma_ignore(layout); | | arma_ignore(layout); |
| | |
|
| arma_stop_logic_error("chol(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("chol(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::chol_band(Mat<eT>& X, const uword KD, const uword layout) | | auxlib::chol_band(Mat<eT>& X, const uword KD, const uword layout) |
| { | | { |
| | |
| skipping to change at line 2622 | | skipping to change at line 2484 |
|---|
| auxlib::hess(Mat<eT>& H, const Base<eT,T1>& X, Col<eT>& tao) | | auxlib::hess(Mat<eT>& H, const Base<eT,T1>& X, Col<eT>& tao) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| H = X.get_ref(); | | H = X.get_ref(); |
| | |
| arma_debug_check( (H.is_square() == false), "hess(): given matrix must be sq
uare sized" ); | | arma_debug_check( (H.is_square() == false), "hess(): given matrix must be sq
uare sized" ); |
| | |
|
| if(H.is_empty()) | | if(H.is_empty()) { return true; } |
| { | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(H); | | arma_debug_assert_blas_size(H); |
| | |
| if(H.n_rows > 2) | | if(H.n_rows > 2) |
| { | | { |
| tao.set_size(H.n_rows-1); | | tao.set_size(H.n_rows-1); |
| | |
| blas_int n = blas_int(H.n_rows); | | blas_int n = blas_int(H.n_rows); |
| blas_int ilo = 1; | | blas_int ilo = 1; |
| blas_int ihi = blas_int(H.n_rows); | | blas_int ihi = blas_int(H.n_rows); |
| | |
| skipping to change at line 2675 | | skipping to change at line 2534 |
|---|
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| R = X.get_ref(); | | R = X.get_ref(); |
| | |
| const uword R_n_rows = R.n_rows; | | const uword R_n_rows = R.n_rows; |
| const uword R_n_cols = R.n_cols; | | const uword R_n_cols = R.n_cols; |
| | |
|
| if(R.is_empty()) | | if(R.is_empty()) { Q.eye(R_n_rows, R_n_rows); return true; } |
| { | | |
| Q.eye(R_n_rows, R_n_rows); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(R); | | arma_debug_assert_blas_size(R); |
| | |
| blas_int m = static_cast<blas_int>(R_n_rows); | | blas_int m = static_cast<blas_int>(R_n_rows); |
| blas_int n = static_cast<blas_int>(R_n_cols); | | blas_int n = static_cast<blas_int>(R_n_cols); |
| blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n)); // take into
account requirements of geqrf() _and_ orgqr()/ungqr() | | blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n)); // take into
account requirements of geqrf() _and_ orgqr()/ungqr() |
| blas_int k = (std::min)(m,n); | | blas_int k = (std::min)(m,n); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<eT> tau( static_cast<uword>(k) ); | | podarray<eT> tau( static_cast<uword>(k) ); |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::geqrf()"); | | arma_extra_debug_print("lapack::geqrf()"); |
| lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), &work_query[0], &lwork_q
uery, &info); | | lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), &work_query[0], &lwork_q
uery, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| podarray<eT> work( static_cast<uword>(lwork_final) ); | | podarray<eT> work( static_cast<uword>(lwork_final) ); |
| | |
| skipping to change at line 2763 | | skipping to change at line 2618 |
|---|
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| if(is_Mat<T1>::value) | | if(is_Mat<T1>::value) |
| { | | { |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); |
| const Mat<eT>& M = tmp.M; | | const Mat<eT>& M = tmp.M; |
| | |
|
| if(M.n_rows < M.n_cols) | | if(M.n_rows < M.n_cols) { return auxlib::qr(Q, R, X); } |
| { | | |
| return auxlib::qr(Q, R, X); | | |
| } | | |
| } | | } |
| | |
| Q = X.get_ref(); | | Q = X.get_ref(); |
| | |
| const uword Q_n_rows = Q.n_rows; | | const uword Q_n_rows = Q.n_rows; |
| const uword Q_n_cols = Q.n_cols; | | const uword Q_n_cols = Q.n_cols; |
| | |
|
| if( Q_n_rows <= Q_n_cols ) | | if( Q_n_rows <= Q_n_cols ) { return auxlib::qr(Q, R, Q); } |
| { | | |
| return auxlib::qr(Q, R, Q); | | |
| } | | |
| | |
|
| if(Q.is_empty()) | | if(Q.is_empty()) { Q.set_size(Q_n_rows, 0); R.set_size(0, Q_n_cols); return |
| { | | true; } |
| Q.set_size(Q_n_rows, 0 ); | | |
| R.set_size(0, Q_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(Q); | | arma_debug_assert_blas_size(Q); |
| | |
| blas_int m = static_cast<blas_int>(Q_n_rows); | | blas_int m = static_cast<blas_int>(Q_n_rows); |
| blas_int n = static_cast<blas_int>(Q_n_cols); | | blas_int n = static_cast<blas_int>(Q_n_cols); |
| blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n)); // take into
account requirements of geqrf() _and_ orgqr()/ungqr() | | blas_int lwork_min = (std::max)(blas_int(1), (std::max)(m,n)); // take into
account requirements of geqrf() _and_ orgqr()/ungqr() |
| blas_int k = (std::min)(m,n); | | blas_int k = (std::min)(m,n); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<eT> tau( static_cast<uword>(k) ); | | podarray<eT> tau( static_cast<uword>(k) ); |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::geqrf()"); | | arma_extra_debug_print("lapack::geqrf()"); |
| lapack::geqrf(&m, &n, Q.memptr(), &m, tau.memptr(), &work_query[0], &lwork_q
uery, &info); | | lapack::geqrf(&m, &n, Q.memptr(), &m, tau.memptr(), &work_query[0], &lwork_q
uery, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| podarray<eT> work( static_cast<uword>(lwork_final) ); | | podarray<eT> work( static_cast<uword>(lwork_final) ); |
| | |
| skipping to change at line 2890 | | skipping to change at line 2734 |
|---|
| blas_int n = static_cast<blas_int>(R_n_cols); | | blas_int n = static_cast<blas_int>(R_n_cols); |
| blas_int lwork_min = (std::max)(blas_int(3*n + 1), (std::max)(m,n)); // tak
e into account requirements of geqp3() and orgqr() | | blas_int lwork_min = (std::max)(blas_int(3*n + 1), (std::max)(m,n)); // tak
e into account requirements of geqp3() and orgqr() |
| blas_int k = (std::min)(m,n); | | blas_int k = (std::min)(m,n); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<eT> tau( static_cast<uword>(k) ); | | podarray<eT> tau( static_cast<uword>(k) ); |
| podarray<blas_int> jpvt( R_n_cols ); | | podarray<blas_int> jpvt( R_n_cols ); |
| | |
| jpvt.zeros(); | | jpvt.zeros(); |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::geqp3()"); | | arma_extra_debug_print("lapack::geqp3()"); |
| lapack::geqp3(&m, &n, R.memptr(), &m, jpvt.memptr(), tau.memptr(), &work_que
ry[0], &lwork_query, &info); | | lapack::geqp3(&m, &n, R.memptr(), &m, jpvt.memptr(), tau.memptr(), &work_que
ry[0], &lwork_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| podarray<eT> work( static_cast<uword>(lwork_final) ); | | podarray<eT> work( static_cast<uword>(lwork_final) ); |
| | |
| skipping to change at line 2982 | | skipping to change at line 2826 |
|---|
| blas_int lwork_min = (std::max)(blas_int(3*n + 1), (std::max)(m,n)); // tak
e into account requirements of geqp3() and ungqr() | | blas_int lwork_min = (std::max)(blas_int(3*n + 1), (std::max)(m,n)); // tak
e into account requirements of geqp3() and ungqr() |
| blas_int k = (std::min)(m,n); | | blas_int k = (std::min)(m,n); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<eT> tau( static_cast<uword>(k) ); | | podarray<eT> tau( static_cast<uword>(k) ); |
| podarray< T> rwork( 2*R_n_cols ); | | podarray< T> rwork( 2*R_n_cols ); |
| podarray<blas_int> jpvt( R_n_cols ); | | podarray<blas_int> jpvt( R_n_cols ); |
| | |
| jpvt.zeros(); | | jpvt.zeros(); |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::geqp3()"); | | arma_extra_debug_print("lapack::geqp3()"); |
| lapack::cx_geqp3(&m, &n, R.memptr(), &m, jpvt.memptr(), tau.memptr(), &work_
query[0], &lwork_query, rwork.memptr(), &info); | | lapack::cx_geqp3(&m, &n, R.memptr(), &m, jpvt.memptr(), tau.memptr(), &work_
query[0], &lwork_query, rwork.memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| podarray<eT> work( static_cast<uword>(lwork_final) ); | | podarray<eT> work( static_cast<uword>(lwork_final) ); |
| | |
| skipping to change at line 3044 | | skipping to change at line 2888 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd(Col<eT>& S, Mat<eT>& A) | | auxlib::svd(Col<eT>& S, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| if(A.is_empty()) { S.reset(); return true; } | | if(A.is_empty()) { S.reset(); return true; } |
| | |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| Mat<eT> U(1, 1, arma_nozeros_indicator()); | | Mat<eT> U(1, 1, arma_nozeros_indicator()); |
| Mat<eT> V(1, A.n_cols, arma_nozeros_indicator()); | | Mat<eT> V(1, A.n_cols, arma_nozeros_indicator()); |
| | |
| char jobu = 'N'; | | char jobu = 'N'; |
| char jobvt = 'N'; | | char jobvt = 'N'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| | |
| skipping to change at line 3065 | | skipping to change at line 2911 |
|---|
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldu = blas_int(U.n_rows); | | blas_int ldu = blas_int(U.n_rows); |
| blas_int ldvt = blas_int(V.n_rows); | | blas_int ldvt = blas_int(V.n_rows); |
| blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); | | blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::gesvd()"); | | arma_extra_debug_print("lapack::gesvd()"); |
| lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); | | lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( work_query[0] ); | | lwork_proposed = static_cast<blas_int>( work_query[0] ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3110 | | skipping to change at line 2956 |
|---|
| auxlib::svd(Col<T>& S, Mat< std::complex<T> >& A) | | auxlib::svd(Col<T>& S, Mat< std::complex<T> >& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
| if(A.is_empty()) { S.reset(); return true; } | | if(A.is_empty()) { S.reset(); return true; } |
| | |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| Mat<eT> U(1, 1, arma_nozeros_indicator()); | | Mat<eT> U(1, 1, arma_nozeros_indicator()); |
| Mat<eT> V(1, A.n_cols, arma_nozeros_indicator()); | | Mat<eT> V(1, A.n_cols, arma_nozeros_indicator()); |
| | |
| char jobu = 'N'; | | char jobu = 'N'; |
| char jobvt = 'N'; | | char jobvt = 'N'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| | |
| skipping to change at line 3133 | | skipping to change at line 2981 |
|---|
| blas_int ldvt = blas_int(V.n_rows); | | blas_int ldvt = blas_int(V.n_rows); |
| blas_int lwork_min = (std::max)( blas_int(1), 2*min_mn+(std::max)(m,n) ); | | blas_int lwork_min = (std::max)( blas_int(1), 2*min_mn+(std::max)(m,n) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<T> rwork( static_cast<uword>(5*min_mn) ); | | podarray<T> rwork( static_cast<uword>(5*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; // query to find optimum size of workspace | | blas_int lwork_query = -1; // query to find optimum size of workspace |
| | |
| arma_extra_debug_print("lapack::cx_gesvd()"); | | arma_extra_debug_print("lapack::cx_gesvd()"); |
| lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); | | lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3174 | | skipping to change at line 3022 |
|---|
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) | | auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(A.n_rows, A.n_rows); S.reset(); V.eye(A.n_cols, A. |
| { | | n_cols); return true; } |
| U.eye(A.n_rows, A.n_rows); | | |
| S.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| V.eye(A.n_cols, A.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); |
| | |
| char jobu = 'A'; | | char jobu = 'A'; |
| char jobvt = 'A'; | | char jobvt = 'A'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| | |
| skipping to change at line 3203 | | skipping to change at line 3047 |
|---|
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldu = blas_int(U.n_rows); | | blas_int ldu = blas_int(U.n_rows); |
| blas_int ldvt = blas_int(V.n_rows); | | blas_int ldvt = blas_int(V.n_rows); |
| blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std:
:max)(m,n)), 5*min_mn ) ); | | blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std:
:max)(m,n)), 5*min_mn ) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
| // query to find optimum size of workspace | | // query to find optimum size of workspace |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::gesvd()"); | | arma_extra_debug_print("lapack::gesvd()"); |
| lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); | | lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( work_query[0] ); | | lwork_proposed = static_cast<blas_int>( work_query[0] ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3253 | | skipping to change at line 3097 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat
< std::complex<T> >& A) | | auxlib::svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, Mat
< std::complex<T> >& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(A.n_rows, A.n_rows); S.reset(); V.eye(A.n_cols, A. |
| { | | n_cols); return true; } |
| U.eye(A.n_rows, A.n_rows); | | |
| S.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| V.eye(A.n_cols, A.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); |
| | |
| char jobu = 'A'; | | char jobu = 'A'; |
| char jobvt = 'A'; | | char jobvt = 'A'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| | |
| skipping to change at line 3284 | | skipping to change at line 3124 |
|---|
| blas_int ldvt = blas_int(V.n_rows); | | blas_int ldvt = blas_int(V.n_rows); |
| blas_int lwork_min = (std::max)( blas_int(1), 2*min_mn + (std::max)(m,n) ); | | blas_int lwork_min = (std::max)( blas_int(1), 2*min_mn + (std::max)(m,n) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<T> rwork( static_cast<uword>(5*min_mn) ); | | podarray<T> rwork( static_cast<uword>(5*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; // query to find optimum size of workspace | | blas_int lwork_query = -1; // query to find optimum size of workspace |
| | |
| arma_extra_debug_print("lapack::cx_gesvd()"); | | arma_extra_debug_print("lapack::cx_gesvd()"); |
| lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); | | lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3331 | | skipping to change at line 3171 |
|---|
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A, const char mode
) | | auxlib::svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A, const char mode
) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(); S.reset(); V.eye(); return true; } |
| { | | |
| U.eye(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| S.reset(); | | |
| V.eye(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| | |
| skipping to change at line 3395 | | skipping to change at line 3231 |
|---|
| | |
| U.set_size( static_cast<uword>(ldu), static_cast<uword>(min_mn) ); | | U.set_size( static_cast<uword>(ldu), static_cast<uword>(min_mn) ); |
| V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n ) ); | | V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n ) ); |
| } | | } |
| | |
| blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); | | blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; // query to find optimum size of workspace | | blas_int lwork_query = -1; // query to find optimum size of workspace |
| | |
| arma_extra_debug_print("lapack::gesvd()"); | | arma_extra_debug_print("lapack::gesvd()"); |
| lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); | | lapack::gesvd<eT>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.m
emptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>(work_query[0]); | | lwork_proposed = static_cast<blas_int>(work_query[0]); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3445 | | skipping to change at line 3281 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V
, Mat< std::complex<T> >& A, const char mode) | | auxlib::svd_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V
, Mat< std::complex<T> >& A, const char mode) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(); S.reset(); V.eye(); return true; } |
| { | | |
| U.eye(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| S.reset(); | | |
| V.eye(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| | |
| skipping to change at line 3511 | | skipping to change at line 3343 |
|---|
| V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) ); | | V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) ); |
| } | | } |
| | |
| blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); | | blas_int lwork_min = (std::max)( blas_int(1), (std::max)( (3*min_mn + (std::
max)(m,n)), 5*min_mn ) ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| podarray<T> rwork( static_cast<uword>(5*min_mn) ); | | podarray<T> rwork( static_cast<uword>(5*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; // query to find optimum size of workspace | | blas_int lwork_query = -1; // query to find optimum size of workspace |
| | |
| arma_extra_debug_print("lapack::cx_gesvd()"); | | arma_extra_debug_print("lapack::cx_gesvd()"); |
| lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); | | lapack::cx_gesvd<T>(&jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U
.memptr(), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr()
, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3561 | | skipping to change at line 3393 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd_dc(Col<eT>& S, Mat<eT>& A) | | auxlib::svd_dc(Col<eT>& S, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| if(A.is_empty()) { S.reset(); return true; } | | if(A.is_empty()) { S.reset(); return true; } |
| | |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| Mat<eT> U(1, 1, arma_nozeros_indicator()); | | Mat<eT> U(1, 1, arma_nozeros_indicator()); |
| Mat<eT> V(1, 1, arma_nozeros_indicator()); | | Mat<eT> V(1, 1, arma_nozeros_indicator()); |
| | |
| char jobz = 'N'; | | char jobz = 'N'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| | |
| skipping to change at line 3584 | | skipping to change at line 3418 |
|---|
| blas_int ldvt = blas_int(V.n_rows); | | blas_int ldvt = blas_int(V.n_rows); |
| blas_int lwork_min = 3*min_mn + (std::max)( max_mn, 7*min_mn ); | | blas_int lwork_min = 3*min_mn + (std::max)( max_mn, 7*min_mn ); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::gesdd()"); | | arma_extra_debug_print("lapack::gesdd()"); |
| lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); | | lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( work_query[0] ); | | lwork_proposed = static_cast<blas_int>( work_query[0] ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3629 | | skipping to change at line 3463 |
|---|
| auxlib::svd_dc(Col<T>& S, Mat< std::complex<T> >& A) | | auxlib::svd_dc(Col<T>& S, Mat< std::complex<T> >& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
| if(A.is_empty()) { S.reset(); return true; } | | if(A.is_empty()) { S.reset(); return true; } |
| | |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| Mat<eT> U(1, 1, arma_nozeros_indicator()); | | Mat<eT> U(1, 1, arma_nozeros_indicator()); |
| Mat<eT> V(1, 1, arma_nozeros_indicator()); | | Mat<eT> V(1, 1, arma_nozeros_indicator()); |
| | |
| char jobz = 'N'; | | char jobz = 'N'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| | |
| skipping to change at line 3653 | | skipping to change at line 3489 |
|---|
| blas_int lwork_min = 2*min_mn + max_mn; | | blas_int lwork_min = 2*min_mn + max_mn; |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<T> rwork( static_cast<uword>(7*min_mn) ); // from LAPACK 3.
8 docs: LAPACK <= v3.6 needs 7*mn | | podarray<T> rwork( static_cast<uword>(7*min_mn) ); // from LAPACK 3.
8 docs: LAPACK <= v3.6 needs 7*mn |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::cx_gesdd()"); | | arma_extra_debug_print("lapack::cx_gesdd()"); |
| lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); | | lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3694 | | skipping to change at line 3530 |
|---|
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::svd_dc(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) | | auxlib::svd_dc(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(A.n_rows, A.n_rows); S.reset(); V.eye(A.n_cols, A. |
| { | | n_cols); return true; } |
| U.eye(A.n_rows, A.n_rows); | | |
| S.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| V.eye(A.n_cols, A.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); |
| | |
| char jobz = 'A'; | | char jobz = 'A'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| | |
| skipping to change at line 3727 | | skipping to change at line 3559 |
|---|
| blas_int lwork2 = 4*min_mn*min_mn + 6*min_mn + max_mn; // as per LAPACK
3.8 docs; consistent with LAPACK 3.4 docs | | blas_int lwork2 = 4*min_mn*min_mn + 6*min_mn + max_mn; // as per LAPACK
3.8 docs; consistent with LAPACK 3.4 docs |
| blas_int lwork_min = (std::max)(lwork1, lwork2); // due to differences bet
ween LAPACK 3.2 and 3.8 | | blas_int lwork_min = (std::max)(lwork1, lwork2); // due to differences bet
ween LAPACK 3.2 and 3.8 |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::gesdd()"); | | arma_extra_debug_print("lapack::gesdd()"); |
| lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); | | lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>(work_query[0]); | | lwork_proposed = static_cast<blas_int>(work_query[0]); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3776 | | skipping to change at line 3608 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd_dc(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V,
Mat< std::complex<T> >& A) | | auxlib::svd_dc(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V,
Mat< std::complex<T> >& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { U.eye(A.n_rows, A.n_rows); S.reset(); V.eye(A.n_cols, A. |
| { | | n_cols); return true; } |
| U.eye(A.n_rows, A.n_rows); | | |
| S.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| V.eye(A.n_cols, A.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); |
| | |
| char jobz = 'A'; | | char jobz = 'A'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| | |
| skipping to change at line 3809 | | skipping to change at line 3637 |
|---|
| blas_int lrwork = min_mn * ((std::max)(5*min_mn+7, 2*max_mn + 2*min_mn+1)
); // as per LAPACK 3.4 docs; LAPACK 3.8 uses 5*min_mn+5 instead of 5*min_mn+7 | | blas_int lrwork = min_mn * ((std::max)(5*min_mn+7, 2*max_mn + 2*min_mn+1)
); // as per LAPACK 3.4 docs; LAPACK 3.8 uses 5*min_mn+5 instead of 5*min_mn+7 |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| podarray<T> rwork( static_cast<uword>(lrwork ) ); | | podarray<T> rwork( static_cast<uword>(lrwork ) ); |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::cx_gesdd()"); | | arma_extra_debug_print("lapack::cx_gesdd()"); |
| lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); | | lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3856 | | skipping to change at line 3684 |
|---|
| | |
| template<typename eT> | | template<typename eT> |
| inline | | inline |
| bool | | bool |
| auxlib::svd_dc_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) | | auxlib::svd_dc_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, Mat<eT>& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char jobz = 'S'; | | char jobz = 'S'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int max_mn = (std::max)(m,n); | | blas_int max_mn = (std::max)(m,n); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldu = m; | | blas_int ldu = m; |
| | |
| skipping to change at line 3890 | | skipping to change at line 3720 |
|---|
| S.set_size( static_cast<uword>(min_mn) ); | | S.set_size( static_cast<uword>(min_mn) ); |
| | |
| U.set_size( static_cast<uword>(m), static_cast<uword>(min_mn) ); | | U.set_size( static_cast<uword>(m), static_cast<uword>(min_mn) ); |
| | |
| V.set_size( static_cast<uword>(min_mn), static_cast<uword>(n) ); | | V.set_size( static_cast<uword>(min_mn), static_cast<uword>(n) ); |
| | |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 1024) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::gesdd()"); | | arma_extra_debug_print("lapack::gesdd()"); |
| lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); | | lapack::gesdd<eT>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(),
&ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, iwork.memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>(work_query[0]); | | lwork_proposed = static_cast<blas_int>(work_query[0]); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 3939 | | skipping to change at line 3769 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::svd_dc_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >
& V, Mat< std::complex<T> >& A) | | auxlib::svd_dc_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >
& V, Mat< std::complex<T> >& A) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
|
| | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char jobz = 'S'; | | char jobz = 'S'; |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int max_mn = (std::max)(m,n); | | blas_int max_mn = (std::max)(m,n); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldu = m; | | blas_int ldu = m; |
| | |
| skipping to change at line 3973 | | skipping to change at line 3805 |
|---|
| | |
| U.set_size( static_cast<uword>(m), static_cast<uword>(min_mn) ); | | U.set_size( static_cast<uword>(m), static_cast<uword>(min_mn) ); |
| | |
| V.set_size( static_cast<uword>(min_mn), static_cast<uword>(n) ); | | V.set_size( static_cast<uword>(min_mn), static_cast<uword>(n) ); |
| | |
| podarray<T> rwork( static_cast<uword>(lrwork ) ); | | podarray<T> rwork( static_cast<uword>(lrwork ) ); |
| podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); | | podarray<blas_int> iwork( static_cast<uword>(8*min_mn) ); |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= 256) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::cx_gesdd()"); | | arma_extra_debug_print("lapack::cx_gesdd()"); |
| lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); | | lapack::cx_gesdd<T>(&jobz, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(
), &ldu, V.memptr(), &ldvt, &work_query[0], &lwork_query, rwork.memptr(), iwork.
memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 4011 | | skipping to change at line 3843 |
|---|
| arma_ignore(U); | | arma_ignore(U); |
| arma_ignore(S); | | arma_ignore(S); |
| arma_ignore(V); | | arma_ignore(V); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_stop_logic_error("svd(): use of LAPACK must be enabled"); | | arma_stop_logic_error("svd(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
|
| //! solve a system of linear equations via explicit inverse (tiny matrices) | | |
| template<typename T1> | | |
| arma_cold | | |
| inline | | |
| bool | | |
| auxlib::solve_square_tiny(Mat<typename T1::elem_type>& out, const Mat<typename T | | |
| 1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr) | | |
| { | | |
| arma_extra_debug_sigprint(); | | |
| | |
| // NOTE: assuming A has size <= 4x4 | | |
| | |
| typedef typename T1::elem_type eT; | | |
| | |
| const uword A_n_rows = A.n_rows; | | |
| | |
| Mat<eT> A_inv(A_n_rows, A_n_rows, arma_nozeros_indicator()); | | |
| | |
| const bool status = op_inv::apply_tiny_noalias(A_inv, A); | | |
| | |
| if(status == false) { return false; } | | |
| | |
| const quasi_unwrap<T1> UB(B_expr.get_ref()); | | |
| const Mat<eT>& B = UB.M; | | |
| | |
| const uword B_n_rows = B.n_rows; | | |
| const uword B_n_cols = B.n_cols; | | |
| | |
| arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the give | | |
| n matrices must be the same" ); | | |
| | |
| if(A.is_empty() || B.is_empty()) | | |
| { | | |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| if(UB.is_alias(out)) | | |
| { | | |
| Mat<eT> tmp(A_n_rows, B_n_cols, arma_nozeros_indicator()); | | |
| | |
| gemm_emul<false,false,false,false>::apply(tmp, A_inv, B); | | |
| | |
| out.steal_mem(tmp); | | |
| } | | |
| else | | |
| { | | |
| out.set_size(A_n_rows, B_n_cols); | | |
| | |
| gemm_emul<false,false,false,false>::apply(out, A_inv, B); | | |
| } | | |
| | |
| return true; | | |
| } | | |
| | | |
| //! solve a system of linear equations via LU decomposition | | //! solve a system of linear equations via LU decomposition |
| template<typename T1> | | template<typename T1> |
| inline | | inline |
| bool | | bool |
| auxlib::solve_square_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::ele
m_type>& A, const Base<typename T1::elem_type,T1>& B_expr) | | auxlib::solve_square_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::ele
m_type>& A, const Base<typename T1::elem_type,T1>& B_expr) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| | |
| const uword A_n_rows = A.n_rows; | | const uword A_n_rows = A.n_rows; |
| | |
|
| if((A_n_rows <= 4) && is_cx<eT>::no) | | |
| { | | |
| const bool status = auxlib::solve_square_tiny(out, A, B_expr.get_ref()); | | |
| | |
| if(status) { return true; } | | |
| } | | |
| | | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the give
n matrices must be the same" ); | | arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the give
n matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | | |
| #if defined(ARMA_USE_ATLAS) | | |
| { | | |
| arma_debug_assert_atlas_size(A); | | |
| | |
| podarray<int> ipiv(A_n_rows + 2); // +2 for paranoia: old versions of Atlas | | |
| might be trashing memory | | |
| | |
| arma_extra_debug_print("atlas::clapack_gesv()"); | | |
| int info = atlas::clapack_gesv<eT>(atlas::CblasColMajor, A_n_rows, B_n_cols, | | |
| A.memptr(), A_n_rows, ipiv.memptr(), out.memptr(), A_n_rows); | | |
| | |
|
| return (info == 0); | | #if defined(ARMA_USE_LAPACK) |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| blas_int n = blas_int(A_n_rows); // assuming A is square | | blas_int n = blas_int(A_n_rows); // assuming A is square |
| blas_int lda = blas_int(A_n_rows); | | blas_int lda = blas_int(A_n_rows); |
| blas_int ldb = blas_int(B_n_rows); | | blas_int ldb = blas_int(B_n_rows); |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| podarray<blas_int> ipiv(A_n_rows + 2); // +2 for paranoia: some versions of
Lapack might be trashing memory | | podarray<blas_int> ipiv(A_n_rows + 2); // +2 for paranoia: some versions of
Lapack might be trashing memory |
| | |
| arma_extra_debug_print("lapack::gesv()"); | | arma_extra_debug_print("lapack::gesv()"); |
| lapack::gesv<eT>(&n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &
ldb, &info); | | lapack::gesv<eT>(&n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &
ldb, &info); |
| | |
| return (info == 0); | | return (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
|
| arma_stop_logic_error("solve(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("solve(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! solve a system of linear equations via LU decomposition with rcond estimate | | //! solve a system of linear equations via LU decomposition with rcond estimate |
| template<typename T1> | | template<typename T1> |
| inline | | inline |
| bool | | bool |
| auxlib::solve_square_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_ty
pe& out_rcond, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type
,T1>& B_expr, const bool allow_ugly) | | auxlib::solve_square_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_ty
pe& out_rcond, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type
,T1>& B_expr, const bool allow_ugly) |
| | |
| skipping to change at line 4154 | | skipping to change at line 3911 |
|---|
| | |
| out_rcond = T(0); | | out_rcond = T(0); |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| char norm_id = '1'; | | char norm_id = '1'; |
| char trans = 'N'; | | char trans = 'N'; |
| blas_int n = blas_int(A.n_rows); // assuming A is square | | blas_int n = blas_int(A.n_rows); // assuming A is square |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldb = blas_int(B_n_rows); | | blas_int ldb = blas_int(B_n_rows); |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| skipping to change at line 4232 | | skipping to change at line 3985 |
|---|
| const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow | | const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow |
| | |
| const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); | | const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); |
| | |
| Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } | | Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } |
| | |
| const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; | | const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| out.set_size(A.n_rows, B.n_cols); | | out.set_size(A.n_rows, B.n_cols); |
| | |
| char fact = (equilibrate) ? 'E' : 'N'; | | char fact = (equilibrate) ? 'E' : 'N'; |
| char trans = 'N'; | | char trans = 'N'; |
| char equed = char(0); | | char equed = char(0); |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| | |
| skipping to change at line 4332 | | skipping to change at line 4081 |
|---|
| const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow | | const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow |
| | |
| const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); | | const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); |
| | |
| Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } | | Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } |
| | |
| const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; | | const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| out.set_size(A.n_rows, B.n_cols); | | out.set_size(A.n_rows, B.n_cols); |
| | |
| char fact = (equilibrate) ? 'E' : 'N'; | | char fact = (equilibrate) ? 'E' : 'N'; |
| char trans = 'N'; | | char trans = 'N'; |
| char equed = char(0); | | char equed = char(0); |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| | |
| skipping to change at line 4436 | | skipping to change at line 4181 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::solve_sympd_fast_common(Mat<typename T1::elem_type>& out, Mat<typename T
1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr) | | auxlib::solve_sympd_fast_common(Mat<typename T1::elem_type>& out, Mat<typename T
1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| | |
| const uword A_n_rows = A.n_rows; | | const uword A_n_rows = A.n_rows; |
| | |
|
| if((A_n_rows <= 4) && is_cx<eT>::no) | | |
| { | | |
| const bool status = auxlib::solve_square_tiny(out, A, B_expr.get_ref()); | | |
| | |
| if(status) { return true; } | | |
| } | | |
| | | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the give
n matrices must be the same" ); | | arma_debug_check( (A_n_rows != B_n_rows), "solve(): number of rows in the give
n matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | | |
| #if defined(ARMA_USE_ATLAS) | | |
| { | | |
| arma_debug_assert_atlas_size(A, out); | | |
| | |
| int info = 0; | | |
| | |
|
| arma_extra_debug_print("atlas::clapack_posv()"); | | #if defined(ARMA_USE_LAPACK) |
| info = atlas::clapack_posv<eT>(atlas::CblasColMajor, atlas::CblasLower, A_n_ | | |
| rows, B_n_cols, A.memptr(), A_n_rows, out.memptr(), B_n_rows); | | |
| | |
| return (info == 0); | | |
| } | | |
| #elif defined(ARMA_USE_LAPACK) | | |
| { | | { |
| arma_debug_assert_blas_size(A, out); | | arma_debug_assert_blas_size(A, out); |
| | |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A_n_rows); // assuming A is square | | blas_int n = blas_int(A_n_rows); // assuming A is square |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int lda = blas_int(A_n_rows); | | blas_int lda = blas_int(A_n_rows); |
| blas_int ldb = blas_int(B_n_rows); | | blas_int ldb = blas_int(B_n_rows); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| arma_extra_debug_print("lapack::posv()"); | | arma_extra_debug_print("lapack::posv()"); |
| lapack::posv<eT>(&uplo, &n, &nrhs, A.memptr(), &lda, out.memptr(), &ldb, &in
fo); | | lapack::posv<eT>(&uplo, &n, &nrhs, A.memptr(), &lda, out.memptr(), &ldb, &in
fo); |
| | |
| return (info == 0); | | return (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out); | | arma_ignore(out); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(B_expr); | | arma_ignore(B_expr); |
|
| arma_stop_logic_error("solve(): use of ATLAS or LAPACK must be enabled"); | | arma_stop_logic_error("solve(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! solve a system of linear equations via Cholesky decomposition with rcond est
imate (real matrices) | | //! solve a system of linear equations via Cholesky decomposition with rcond est
imate (real matrices) |
| template<typename T1> | | template<typename T1> |
| inline | | inline |
| bool | | bool |
|
| auxlib::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) | | auxlib::solve_sympd_rcond(Mat<typename T1::pod_type>& out, bool& out_sympd_state
, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const Base<ty
pename T1::pod_type,T1>& B_expr, const bool allow_ugly) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
|
| out_rcond = T(0); | | out_sympd_state = false; |
| | out_rcond = T(0); |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A, out); | | arma_debug_assert_blas_size(A, out); |
| | |
| char norm_id = '1'; | | char norm_id = '1'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A.n_rows); // assuming A is square | | blas_int n = blas_int(A.n_rows); // assuming A is square |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| T norm_val = T(0); | | T norm_val = T(0); |
| | |
| podarray<T> work(A.n_rows); | | podarray<T> work(A.n_rows); |
| | |
| arma_extra_debug_print("lapack::lansy()"); | | arma_extra_debug_print("lapack::lansy()"); |
| norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; | | norm_val = lapack::lansy(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
|
| | out_sympd_state = true; |
| | |
| arma_extra_debug_print("lapack::potrs()"); | | arma_extra_debug_print("lapack::potrs()"); |
| lapack::potrs<eT>(&uplo, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info)
; | | lapack::potrs<eT>(&uplo, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info)
; |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); | | out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); |
| | |
| if( (allow_ugly == false) && (out_rcond < auxlib::epsilon_lapack(A)) ) { re
turn false; } | | if( (allow_ugly == false) && (out_rcond < auxlib::epsilon_lapack(A)) ) { re
turn false; } |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out); | | arma_ignore(out); |
|
| | arma_ignore(out_sympd_state); |
| arma_ignore(out_rcond); | | arma_ignore(out_rcond); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(B_expr); | | arma_ignore(B_expr); |
| arma_ignore(allow_ugly); | | arma_ignore(allow_ugly); |
| arma_stop_logic_error("solve(): use of LAPACK must be enabled"); | | arma_stop_logic_error("solve(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| //! solve a system of linear equations via Cholesky decomposition with rcond est
imate (complex matrices) | | //! solve a system of linear equations via Cholesky decomposition with rcond est
imate (complex matrices) |
| template<typename T1> | | template<typename T1> |
| inline | | inline |
| bool | | bool |
|
| auxlib::solve_sympd_rcond(Mat< std::complex<typename T1::pod_type> >& out, typen
ame T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, cons
t Base< std::complex<typename T1::pod_type>,T1>& B_expr, const bool allow_ugly) | | auxlib::solve_sympd_rcond(Mat< std::complex<typename T1::pod_type> >& out, bool&
out_sympd_state, typename T1::pod_type& out_rcond, Mat< std::complex<typename T
1::pod_type> >& A, const Base< std::complex<typename T1::pod_type>,T1>& B_expr,
const bool allow_ugly) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_CRIPPLED_LAPACK) | | #if defined(ARMA_CRIPPLED_LAPACK) |
| { | | { |
| arma_extra_debug_print("auxlib::solve_sympd_rcond(): redirecting to auxlib::
solve_square_rcond() due to crippled LAPACK"); | | arma_extra_debug_print("auxlib::solve_sympd_rcond(): redirecting to auxlib::
solve_square_rcond() due to crippled LAPACK"); |
| | |
|
| | out_sympd_state = false; |
| | |
| return auxlib::solve_square_rcond(out, out_rcond, A, B_expr, allow_ugly); | | return auxlib::solve_square_rcond(out, out_rcond, A, B_expr, allow_ugly); |
| } | | } |
| #elif defined(ARMA_USE_LAPACK) | | #elif defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
|
| out_rcond = T(0); | | out_sympd_state = false; |
| | out_rcond = T(0); |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A, out); | | arma_debug_assert_blas_size(A, out); |
| | |
| char norm_id = '1'; | | char norm_id = '1'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| blas_int n = blas_int(A.n_rows); // assuming A is square | | blas_int n = blas_int(A.n_rows); // assuming A is square |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| T norm_val = T(0); | | T norm_val = T(0); |
| | |
| podarray<T> work(A.n_rows); | | podarray<T> work(A.n_rows); |
| | |
| arma_extra_debug_print("lapack::lanhe()"); | | arma_extra_debug_print("lapack::lanhe()"); |
| norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; | | norm_val = lapack::lanhe(&norm_id, &uplo, &n, A.memptr(), &n, work.memptr())
; |
| | |
| arma_extra_debug_print("lapack::potrf()"); | | arma_extra_debug_print("lapack::potrf()"); |
| lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info); | | lapack::potrf<eT>(&uplo, &n, A.memptr(), &n, &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
|
| | out_sympd_state = true; |
| | |
| arma_extra_debug_print("lapack::potrs()"); | | arma_extra_debug_print("lapack::potrs()"); |
| lapack::potrs<eT>(&uplo, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info)
; | | lapack::potrs<eT>(&uplo, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info)
; |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); | | out_rcond = auxlib::lu_rcond_sympd<T>(A, norm_val); |
| | |
| if( (allow_ugly == false) && (out_rcond < auxlib::epsilon_lapack(A)) ) { re
turn false; } | | if( (allow_ugly == false) && (out_rcond < auxlib::epsilon_lapack(A)) ) { re
turn false; } |
| | |
| return true; | | return true; |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out); | | arma_ignore(out); |
|
| | arma_ignore(out_sympd_state); |
| arma_ignore(out_rcond); | | arma_ignore(out_rcond); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(B_expr); | | arma_ignore(B_expr); |
| arma_ignore(allow_ugly); | | arma_ignore(allow_ugly); |
| arma_stop_logic_error("solve(): use of LAPACK must be enabled"); | | arma_stop_logic_error("solve(): use of LAPACK must be enabled"); |
| return false; | | return false; |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| | |
| skipping to change at line 4668 | | skipping to change at line 4393 |
|---|
| const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow | | const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow |
| | |
| const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); | | const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); |
| | |
| Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } | | Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } |
| | |
| const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; | | const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| out.set_size(A.n_rows, B.n_cols); | | out.set_size(A.n_rows, B.n_cols); |
| | |
| char fact = (equilibrate) ? 'E' : 'N'; | | char fact = (equilibrate) ? 'E' : 'N'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| char equed = char(0); | | char equed = char(0); |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| | |
| skipping to change at line 4704 | | skipping to change at line 4425 |
|---|
| podarray<eT> BERR( B.n_cols); | | podarray<eT> BERR( B.n_cols); |
| podarray<eT> WORK(3*A.n_rows); | | podarray<eT> WORK(3*A.n_rows); |
| podarray<blas_int> IWORK( A.n_rows); | | podarray<blas_int> IWORK( A.n_rows); |
| | |
| arma_extra_debug_print("lapack::posvx()"); | | arma_extra_debug_print("lapack::posvx()"); |
| lapack::posvx(&fact, &uplo, &n, &nrhs, A.memptr(), &lda, AF.memptr(), &ldaf,
&equed, S.memptr(), const_cast<eT*>(B.memptr()), &ldb, out.memptr(), &ldx, &rco
nd, FERR.memptr(), BERR.memptr(), WORK.memptr(), IWORK.memptr(), &info); | | lapack::posvx(&fact, &uplo, &n, &nrhs, A.memptr(), &lda, AF.memptr(), &ldaf,
&equed, S.memptr(), const_cast<eT*>(B.memptr()), &ldb, out.memptr(), &ldx, &rco
nd, FERR.memptr(), BERR.memptr(), WORK.memptr(), IWORK.memptr(), &info); |
| | |
| // NOTE: using const_cast<eT*>(B.memptr()) to allow B to be overwritten for
equilibration; | | // NOTE: using const_cast<eT*>(B.memptr()) to allow B to be overwritten for
equilibration; |
| // NOTE: B is created as a copy of B_expr if equilibration is enabled; other
wise B is a reference to B_expr | | // NOTE: B is created as a copy of B_expr if equilibration is enabled; other
wise B is a reference to B_expr |
| | |
|
| | // NOTE: lapack::posvx() sets rcond to zero if A is not sympd |
| out_rcond = rcond; | | out_rcond = rcond; |
| | |
| return (allow_ugly) ? ((info == 0) || (info == (n+1))) : (info == 0); | | return (allow_ugly) ? ((info == 0) || (info == (n+1))) : (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out); | | arma_ignore(out); |
| arma_ignore(out_rcond); | | arma_ignore(out_rcond); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(B_expr); | | arma_ignore(B_expr); |
| | |
| skipping to change at line 4755 | | skipping to change at line 4477 |
|---|
| const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow | | const Mat<eT>& UB_M_as_Mat = UB.M; // so we don't confuse the ?: operator b
elow |
| | |
| const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); | | const bool use_copy = ((equilibrate && UB.is_const) || UB.is_alias(out)); |
| | |
| Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } | | Mat<eT> B_tmp; if(use_copy) { B_tmp = UB_M_as_Mat; } |
| | |
| const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; | | const Mat<eT>& B = (use_copy) ? B_tmp : UB_M_as_Mat; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| out.set_size(A.n_rows, B.n_cols); | | out.set_size(A.n_rows, B.n_cols); |
| | |
| char fact = (equilibrate) ? 'E' : 'N'; | | char fact = (equilibrate) ? 'E' : 'N'; |
| char uplo = 'L'; | | char uplo = 'L'; |
| char equed = char(0); | | char equed = char(0); |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| | |
| skipping to change at line 4791 | | skipping to change at line 4509 |
|---|
| podarray< T> BERR( B.n_cols); | | podarray< T> BERR( B.n_cols); |
| podarray<eT> WORK(2*A.n_rows); | | podarray<eT> WORK(2*A.n_rows); |
| podarray< T> RWORK( A.n_rows); | | podarray< T> RWORK( A.n_rows); |
| | |
| arma_extra_debug_print("lapack::cx_posvx()"); | | arma_extra_debug_print("lapack::cx_posvx()"); |
| lapack::cx_posvx(&fact, &uplo, &n, &nrhs, A.memptr(), &lda, AF.memptr(), &ld
af, &equed, S.memptr(), const_cast<eT*>(B.memptr()), &ldb, out.memptr(), &ldx, &
rcond, FERR.memptr(), BERR.memptr(), WORK.memptr(), RWORK.memptr(), &info); | | lapack::cx_posvx(&fact, &uplo, &n, &nrhs, A.memptr(), &lda, AF.memptr(), &ld
af, &equed, S.memptr(), const_cast<eT*>(B.memptr()), &ldb, out.memptr(), &ldx, &
rcond, FERR.memptr(), BERR.memptr(), WORK.memptr(), RWORK.memptr(), &info); |
| | |
| // NOTE: using const_cast<eT*>(B.memptr()) to allow B to be overwritten for
equilibration; | | // NOTE: using const_cast<eT*>(B.memptr()) to allow B to be overwritten for
equilibration; |
| // NOTE: B is created as a copy of B_expr if equilibration is enabled; other
wise B is a reference to B_expr | | // NOTE: B is created as a copy of B_expr if equilibration is enabled; other
wise B is a reference to B_expr |
| | |
|
| | // NOTE: lapack::cx_posvx() sets rcond to zero if A is not sympd |
| out_rcond = rcond; | | out_rcond = rcond; |
| | |
| return (allow_ugly) ? ((info == 0) || (info == (n+1))) : (info == 0); | | return (allow_ugly) ? ((info == 0) || (info == (n+1))) : (info == 0); |
| } | | } |
| #else | | #else |
| { | | { |
| arma_ignore(out); | | arma_ignore(out); |
| arma_ignore(out_rcond); | | arma_ignore(out_rcond); |
| arma_ignore(A); | | arma_ignore(A); |
| arma_ignore(B_expr); | | arma_ignore(B_expr); |
| | |
| skipping to change at line 4826 | | skipping to change at line 4545 |
|---|
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| | |
| const unwrap<T1> U(B_expr.get_ref()); | | const unwrap<T1> U(B_expr.get_ref()); |
| const Mat<eT>& B = U.M; | | const Mat<eT>& B = U.M; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_cols, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); | | Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); |
| | |
| if(arma::size(tmp) == arma::size(B)) | | if(arma::size(tmp) == arma::size(B)) |
| { | | { |
| tmp = B; | | tmp = B; |
| } | | } |
| else | | else |
| | |
| skipping to change at line 4858 | | skipping to change at line 4573 |
|---|
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldb = blas_int(tmp.n_rows); | | blas_int ldb = blas_int(tmp.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int lwork_min = (std::max)(blas_int(1), min_mn + (std::max)(min_mn, nr
hs)); | | blas_int lwork_min = (std::max)(blas_int(1), min_mn + (std::max)(min_mn, nr
hs)); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= ((is_cx<eT>::yes) ? uword(256) : uword(1024))) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::gels()"); | | arma_extra_debug_print("lapack::gels()"); |
| lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &
ldb, &work_query[0], &lwork_query, &info ); | | lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &
ldb, &work_query[0], &lwork_query, &info ); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 4922 | | skipping to change at line 4637 |
|---|
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| | |
| out_rcond = T(0); | | out_rcond = T(0); |
| | |
| const unwrap<T1> U(B_expr.get_ref()); | | const unwrap<T1> U(B_expr.get_ref()); |
| const Mat<eT>& B = U.M; | | const Mat<eT>& B = U.M; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_cols, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); | | Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); |
| | |
| if(arma::size(tmp) == arma::size(B)) | | if(arma::size(tmp) == arma::size(B)) |
| { | | { |
| tmp = B; | | tmp = B; |
| } | | } |
| else | | else |
| | |
| skipping to change at line 4954 | | skipping to change at line 4665 |
|---|
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldb = blas_int(tmp.n_rows); | | blas_int ldb = blas_int(tmp.n_rows); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| blas_int min_mn = (std::min)(m,n); | | blas_int min_mn = (std::min)(m,n); |
| blas_int lwork_min = (std::max)(blas_int(1), min_mn + (std::max)(min_mn, nr
hs)); | | blas_int lwork_min = (std::max)(blas_int(1), min_mn + (std::max)(min_mn, nr
hs)); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| blas_int lwork_proposed = 0; | | blas_int lwork_proposed = 0; |
| | |
|
| if((m*n) >= 1024) | | if(A.n_elem >= ((is_cx<eT>::yes) ? uword(256) : uword(1024))) |
| { | | { |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = -1; | | blas_int lwork_query = -1; |
| | |
| arma_extra_debug_print("lapack::gels()"); | | arma_extra_debug_print("lapack::gels()"); |
| lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &
ldb, &work_query[0], &lwork_query, &info ); | | lapack::gels<eT>( &trans, &m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &
ldb, &work_query[0], &lwork_query, &info ); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); | | lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query[0]) ); |
| } | | } |
| | |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 5062 | | skipping to change at line 4773 |
|---|
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type eT; | | typedef typename T1::pod_type eT; |
| | |
| const unwrap<T1> U(B_expr.get_ref()); | | const unwrap<T1> U(B_expr.get_ref()); |
| const Mat<eT>& B = U.M; | | const Mat<eT>& B = U.M; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_cols, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B.n_cols); | | |
| return true; | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); | | Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); |
| | |
| if(arma::size(tmp) == arma::size(B)) | | if(arma::size(tmp) == arma::size(B)) |
| { | | { |
| tmp = B; | | tmp = B; |
| } | | } |
| else | | else |
| | |
| skipping to change at line 5088 | | skipping to change at line 4798 |
|---|
| tmp.zeros(); | | tmp.zeros(); |
| tmp(0,0, arma::size(B)) = B; | | tmp(0,0, arma::size(B)) = B; |
| } | | } |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m, n); | | blas_int min_mn = (std::min)(m, n); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldb = blas_int(tmp.n_rows); | | blas_int ldb = blas_int(tmp.n_rows); |
|
| eT rcond = eT(-1); // -1 means "use machine precision" | | //eT rcond = eT(-1); // -1 means "use machine precision" |
| | eT rcond = (std::max)(A.n_rows, A.n_cols) * std::numeric_limits<eT>:: |
| | epsilon(); |
| blas_int rank = blas_int(0); | | blas_int rank = blas_int(0); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| podarray<eT> S( static_cast<uword>(min_mn) ); | | podarray<eT> S( static_cast<uword>(min_mn) ); |
| | |
| // NOTE: with LAPACK 3.8, can use the workspace query to also obtain liwork, | | // NOTE: with LAPACK 3.8, can use the workspace query to also obtain liwork, |
| // NOTE: which makes the call to lapack::laenv() redundant | | // NOTE: which makes the call to lapack::laenv() redundant |
| | |
| blas_int ispec = blas_int(9); | | blas_int ispec = blas_int(9); |
| | |
| | |
| skipping to change at line 5122 | | skipping to change at line 4833 |
|---|
| blas_int smlsiz = (std::max)( blas_int(25), laenv_result ); | | blas_int smlsiz = (std::max)( blas_int(25), laenv_result ); |
| blas_int smlsiz_p1 = blas_int(1) + smlsiz; | | blas_int smlsiz_p1 = blas_int(1) + smlsiz; |
| | |
| blas_int nlvl = (std::max)( blas_int(0), blas_int(1) + blas_int( std::log(
double(min_mn) / double(smlsiz_p1))/double(0.69314718055994530942) ) ); | | blas_int nlvl = (std::max)( blas_int(0), blas_int(1) + blas_int( std::log(
double(min_mn) / double(smlsiz_p1))/double(0.69314718055994530942) ) ); |
| blas_int liwork = (std::max)( blas_int(1), (blas_int(3)*min_mn*nlvl + blas_i
nt(11)*min_mn) ); | | blas_int liwork = (std::max)( blas_int(1), (blas_int(3)*min_mn*nlvl + blas_i
nt(11)*min_mn) ); |
| | |
| podarray<blas_int> iwork( static_cast<uword>(liwork) ); | | podarray<blas_int> iwork( static_cast<uword>(liwork) ); |
| | |
| blas_int lwork_min = blas_int(12)*min_mn + blas_int(2)*min_mn*smlsiz + blas_
int(8)*min_mn*nlvl + min_mn*nrhs + smlsiz_p1*smlsiz_p1; | | blas_int lwork_min = blas_int(12)*min_mn + blas_int(2)*min_mn*smlsiz + blas_
int(8)*min_mn*nlvl + min_mn*nrhs + smlsiz_p1*smlsiz_p1; |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::gelsd()"); | | arma_extra_debug_print("lapack::gelsd()"); |
| lapack::gelsd(&m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, S.memptr(
), &rcond, &rank, &work_query[0], &lwork_query, iwork.memptr(), &info); | | lapack::gelsd(&m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, S.memptr(
), &rcond, &rank, &work_query[0], &lwork_query, iwork.memptr(), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| // NOTE: in LAPACK 3.8, iwork[0] returns the minimum liwork | | // NOTE: in LAPACK 3.8, iwork[0] returns the minimum liwork |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real(work_query
[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| skipping to change at line 5181 | | skipping to change at line 4892 |
|---|
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| const unwrap<T1> U(B_expr.get_ref()); | | const unwrap<T1> U(B_expr.get_ref()); |
| const Mat<eT>& B = U.M; | | const Mat<eT>& B = U.M; |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_cols, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_cols, B.n_cols); | | |
| return true; | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| } | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| | |
| arma_debug_assert_blas_size(A,B); | | arma_debug_assert_blas_size(A,B); |
| | |
| Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); | | Mat<eT> tmp( (std::max)(A.n_rows, A.n_cols), B.n_cols, arma_nozeros_indicato
r() ); |
| | |
| if(arma::size(tmp) == arma::size(B)) | | if(arma::size(tmp) == arma::size(B)) |
| { | | { |
| tmp = B; | | tmp = B; |
| } | | } |
| else | | else |
| | |
| skipping to change at line 5207 | | skipping to change at line 4917 |
|---|
| tmp.zeros(); | | tmp.zeros(); |
| tmp(0,0, arma::size(B)) = B; | | tmp(0,0, arma::size(B)) = B; |
| } | | } |
| | |
| blas_int m = blas_int(A.n_rows); | | blas_int m = blas_int(A.n_rows); |
| blas_int n = blas_int(A.n_cols); | | blas_int n = blas_int(A.n_cols); |
| blas_int min_mn = (std::min)(m, n); | | blas_int min_mn = (std::min)(m, n); |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); |
| blas_int lda = blas_int(A.n_rows); | | blas_int lda = blas_int(A.n_rows); |
| blas_int ldb = blas_int(tmp.n_rows); | | blas_int ldb = blas_int(tmp.n_rows); |
|
| T rcond = T(-1); // -1 means "use machine precision" | | //T rcond = T(-1); // -1 means "use machine precision" |
| | T rcond = (std::max)(A.n_rows, A.n_cols) * std::numeric_limits<T>::e |
| | psilon(); |
| blas_int rank = blas_int(0); | | blas_int rank = blas_int(0); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| podarray<T> S( static_cast<uword>(min_mn) ); | | podarray<T> S( static_cast<uword>(min_mn) ); |
| | |
| blas_int ispec = blas_int(9); | | blas_int ispec = blas_int(9); |
| | |
| const char* const_name = (is_float<T>::value) ? "CGELSD" : "ZGELSD"; | | const char* const_name = (is_float<T>::value) ? "CGELSD" : "ZGELSD"; |
| const char* const_opts = " "; | | const char* const_opts = " "; |
| | |
| | |
| skipping to change at line 5244 | | skipping to change at line 4955 |
|---|
| ? blas_int(10)*n + blas_int(2)*n*smlsiz + blas_int(8)*n*nlvl + blas_int(3)
*smlsiz*nrhs + (std::max)( (smlsiz_p1)*(smlsiz_p1), n*(blas_int(1)+nrhs) + blas_
int(2)*nrhs ) | | ? blas_int(10)*n + blas_int(2)*n*smlsiz + blas_int(8)*n*nlvl + blas_int(3)
*smlsiz*nrhs + (std::max)( (smlsiz_p1)*(smlsiz_p1), n*(blas_int(1)+nrhs) + blas_
int(2)*nrhs ) |
| : blas_int(10)*m + blas_int(2)*m*smlsiz + blas_int(8)*m*nlvl + blas_int(3)
*smlsiz*nrhs + (std::max)( (smlsiz_p1)*(smlsiz_p1), n*(blas_int(1)+nrhs) + blas_
int(2)*nrhs ); | | : blas_int(10)*m + blas_int(2)*m*smlsiz + blas_int(8)*m*nlvl + blas_int(3)
*smlsiz*nrhs + (std::max)( (smlsiz_p1)*(smlsiz_p1), n*(blas_int(1)+nrhs) + blas_
int(2)*nrhs ); |
| | |
| blas_int liwork = (std::max)( blas_int(1), (blas_int(3)*blas_int(min_mn)*nlv
l + blas_int(11)*blas_int(min_mn)) ); | | blas_int liwork = (std::max)( blas_int(1), (blas_int(3)*blas_int(min_mn)*nlv
l + blas_int(11)*blas_int(min_mn)) ); |
| | |
| podarray<T> rwork( static_cast<uword>(lrwork) ); | | podarray<T> rwork( static_cast<uword>(lrwork) ); |
| podarray<blas_int> iwork( static_cast<uword>(liwork) ); | | podarray<blas_int> iwork( static_cast<uword>(liwork) ); |
| | |
| blas_int lwork_min = 2*min_mn + min_mn*nrhs; | | blas_int lwork_min = 2*min_mn + min_mn*nrhs; |
| | |
|
| eT work_query[2]; | | eT work_query[2] = {}; |
| blas_int lwork_query = blas_int(-1); | | blas_int lwork_query = blas_int(-1); |
| | |
| arma_extra_debug_print("lapack::cx_gelsd()"); | | arma_extra_debug_print("lapack::cx_gelsd()"); |
| lapack::cx_gelsd(&m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, S.memp
tr(), &rcond, &rank, &work_query[0], &lwork_query, rwork.memptr(), iwork.memptr(
), &info); | | lapack::cx_gelsd(&m, &n, &nrhs, A.memptr(), &lda, tmp.memptr(), &ldb, S.memp
tr(), &rcond, &rank, &work_query[0], &lwork_query, rwork.memptr(), iwork.memptr(
), &info); |
| | |
| if(info != 0) { return false; } | | if(info != 0) { return false; } |
| | |
| blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real( work_quer
y[0]) ); | | blas_int lwork_proposed = static_cast<blas_int>( access::tmp_real( work_quer
y[0]) ); |
| blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); | | blas_int lwork_final = (std::max)(lwork_proposed, lwork_min); |
| | |
| podarray<eT> work( static_cast<uword>(lwork_final) ); | | podarray<eT> work( static_cast<uword>(lwork_final) ); |
| | |
| skipping to change at line 5300 | | skipping to change at line 5011 |
|---|
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,out); | | arma_debug_assert_blas_size(A,out); |
| | |
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; |
| char trans = 'N'; | | char trans = 'N'; |
| char diag = 'N'; | | char diag = 'N'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| | |
| skipping to change at line 5352 | | skipping to change at line 5059 |
|---|
| | |
| out_rcond = T(0); | | out_rcond = T(0); |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_cols, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_cols, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A,out); | | arma_debug_assert_blas_size(A,out); |
| | |
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; |
| char trans = 'N'; | | char trans = 'N'; |
| char diag = 'N'; | | char diag = 'N'; |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int info = 0; | | blas_int info = 0; |
| | |
| | |
| skipping to change at line 5447 | | skipping to change at line 5150 |
|---|
| { | | { |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_rows, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_rows, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| // for gbsv, matrix AB size: 2*KL+KU+1 x N; band representation of A stored
in rows KL+1 to 2*KL+KU+1 (note: fortran counts from 1) | | // for gbsv, matrix AB size: 2*KL+KU+1 x N; band representation of A stored
in rows KL+1 to 2*KL+KU+1 (note: fortran counts from 1) |
| | |
| Mat<eT> AB; | | Mat<eT> AB; |
| band_helper::compress(AB, A, KL, KU, true); | | band_helper::compress(AB, A, KL, KU, true); |
| | |
| const uword N = AB.n_cols; // order of the original square matrix A | | const uword N = AB.n_cols; // order of the original square matrix A |
| | |
| arma_debug_assert_blas_size(AB,out); | | arma_debug_assert_blas_size(AB,out); |
| | |
| | |
| skipping to change at line 5549 | | skipping to change at line 5248 |
|---|
| | |
| out_rcond = T(0); | | out_rcond = T(0); |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_rows, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_rows, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| // for gbtrf, matrix AB size: 2*KL+KU+1 x N; band representation of A stored
in rows KL+1 to 2*KL+KU+1 (note: fortran counts from 1) | | // for gbtrf, matrix AB size: 2*KL+KU+1 x N; band representation of A stored
in rows KL+1 to 2*KL+KU+1 (note: fortran counts from 1) |
| | |
| Mat<eT> AB; | | Mat<eT> AB; |
| band_helper::compress(AB, A, KL, KU, true); | | band_helper::compress(AB, A, KL, KU, true); |
| | |
| const uword N = AB.n_cols; // order of the original square matrix A | | const uword N = AB.n_cols; // order of the original square matrix A |
| | |
| arma_debug_assert_blas_size(AB,out); | | arma_debug_assert_blas_size(AB,out); |
| | |
| | |
| skipping to change at line 5628 | | skipping to change at line 5323 |
|---|
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type eT; | | typedef typename T1::pod_type eT; |
| | |
| Mat<eT> B = B_expr.get_ref(); // B is overwritten | | Mat<eT> B = B_expr.get_ref(); // B is overwritten |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| // for gbsvx, matrix AB size: KL+KU+1 x N; band representation of A stored i
n rows 1 to KL+KU+1 (note: fortran counts from 1) | | // for gbsvx, matrix AB size: KL+KU+1 x N; band representation of A stored i
n rows 1 to KL+KU+1 (note: fortran counts from 1) |
| | |
| Mat<eT> AB; | | Mat<eT> AB; |
| band_helper::compress(AB, A, KL, KU, false); | | band_helper::compress(AB, A, KL, KU, false); |
| | |
| const uword N = AB.n_cols; | | const uword N = AB.n_cols; |
| | |
| arma_debug_assert_blas_size(AB,B); | | arma_debug_assert_blas_size(AB,B); |
| | |
| | |
| skipping to change at line 5735 | | skipping to change at line 5426 |
|---|
| } | | } |
| #elif defined(ARMA_USE_LAPACK) | | #elif defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| Mat<eT> B = B_expr.get_ref(); // B is overwritten | | Mat<eT> B = B_expr.get_ref(); // B is overwritten |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) { out.zeros(A.n_rows, B.n_cols); return tr |
| { | | ue; } |
| out.zeros(A.n_rows, B.n_cols); | | |
| return true; | | |
| } | | |
| | |
| // for gbsvx, matrix AB size: KL+KU+1 x N; band representation of A stored i
n rows 1 to KL+KU+1 (note: fortran counts from 1) | | // for gbsvx, matrix AB size: KL+KU+1 x N; band representation of A stored i
n rows 1 to KL+KU+1 (note: fortran counts from 1) |
| | |
| Mat<eT> AB; | | Mat<eT> AB; |
| band_helper::compress(AB, A, KL, KU, false); | | band_helper::compress(AB, A, KL, KU, false); |
| | |
| const uword N = AB.n_cols; | | const uword N = AB.n_cols; |
| | |
| arma_debug_assert_blas_size(AB,B); | | arma_debug_assert_blas_size(AB,B); |
| | |
| | |
| skipping to change at line 5867 | | skipping to change at line 5554 |
|---|
| { | | { |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; |
| | |
| out = B_expr.get_ref(); | | out = B_expr.get_ref(); |
| | |
| const uword B_n_rows = out.n_rows; | | const uword B_n_rows = out.n_rows; |
| const uword B_n_cols = out.n_cols; | | const uword B_n_cols = out.n_cols; |
| | |
| arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); | | arma_debug_check( (A.n_rows != B_n_rows), "solve(): number of rows in the gi
ven matrices must be the same" ); |
| | |
|
| if(A.is_empty() || out.is_empty()) | | if(A.is_empty() || out.is_empty()) { out.zeros(A.n_rows, B_n_cols); return |
| { | | true; } |
| out.zeros(A.n_rows, B_n_cols); | | |
| return true; | | |
| } | | |
| | |
| Mat<eT> tridiag; | | Mat<eT> tridiag; |
| band_helper::extract_tridiag(tridiag, A); | | band_helper::extract_tridiag(tridiag, A); |
| | |
| arma_debug_assert_blas_size(tridiag, out); | | arma_debug_assert_blas_size(tridiag, out); |
| | |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); |
| blas_int nrhs = blas_int(B_n_cols); | | blas_int nrhs = blas_int(B_n_cols); |
| blas_int ldb = blas_int(B_n_rows); | | blas_int ldb = blas_int(B_n_rows); |
| blas_int info = blas_int(0); | | blas_int info = blas_int(0); |
| | |
| skipping to change at line 5915 | | skipping to change at line 5598 |
|---|
| auxlib::schur(Mat<eT>& U, Mat<eT>& S, const Base<eT,T1>& X, const bool calc_U) | | auxlib::schur(Mat<eT>& U, Mat<eT>& S, const Base<eT,T1>& X, const bool calc_U) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| S = X.get_ref(); | | S = X.get_ref(); |
| | |
| arma_debug_check( (S.is_square() == false), "schur(): given matrix must be s
quare sized" ); | | arma_debug_check( (S.is_square() == false), "schur(): given matrix must be s
quare sized" ); |
| | |
|
| if(S.is_empty()) | | if(S.is_empty()) { U.reset(); S.reset(); return true; } |
| { | | |
| U.reset(); | | |
| S.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(S); | | arma_debug_assert_blas_size(S); |
| | |
| const uword S_n_rows = S.n_rows; | | const uword S_n_rows = S.n_rows; |
| | |
| if(calc_U) { U.set_size(S_n_rows, S_n_rows); } else { U.set_size(1,1); } | | if(calc_U) { U.set_size(S_n_rows, S_n_rows); } else { U.set_size(1,1); } |
| | |
| char jobvs = calc_U ? 'V' : 'N'; | | char jobvs = calc_U ? 'V' : 'N'; |
| char sort = 'N'; | | char sort = 'N'; |
| void* select = 0; | | void* select = 0; |
| | |
| skipping to change at line 5985 | | skipping to change at line 5663 |
|---|
| inline | | inline |
| bool | | bool |
| auxlib::schur(Mat< std::complex<T> >& U, Mat< std::complex<T> >& S, const bool c
alc_U) | | auxlib::schur(Mat< std::complex<T> >& U, Mat< std::complex<T> >& S, const bool c
alc_U) |
| { | | { |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); |
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; |
| | |
|
| if(S.is_empty()) | | if(S.is_empty()) { U.reset(); S.reset(); return true; } |
| { | | |
| U.reset(); | | |
| S.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(S); | | arma_debug_assert_blas_size(S); |
| | |
| const uword S_n_rows = S.n_rows; | | const uword S_n_rows = S.n_rows; |
| | |
| if(calc_U) { U.set_size(S_n_rows, S_n_rows); } else { U.set_size(1,1); } | | if(calc_U) { U.set_size(S_n_rows, S_n_rows); } else { U.set_size(1,1); } |
| | |
| char jobvs = calc_U ? 'V' : 'N'; | | char jobvs = calc_U ? 'V' : 'N'; |
| char sort = 'N'; | | char sort = 'N'; |
| void* select = 0; | | void* select = 0; |
| | |
| skipping to change at line 6107 | | skipping to change at line 5780 |
|---|
| | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) |
| { | | { |
| A = X_expr.get_ref(); | | A = X_expr.get_ref(); |
| B = Y_expr.get_ref(); | | B = Y_expr.get_ref(); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "q
z(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "q
z(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "qz(): given matrices must have th
e same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "qz(): given matrices must have th
e same size" ); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { A.reset(); B.reset(); vsl.reset(); vsr.reset(); return |
| { | | true; } |
| A.reset(); | | |
| B.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| vsl.reset(); | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| vsr.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| vsl.set_size(A.n_rows, A.n_rows); | | vsl.set_size(A.n_rows, A.n_rows); |
| vsr.set_size(A.n_rows, A.n_rows); | | vsr.set_size(A.n_rows, A.n_rows); |
| | |
| char jobvsl = 'V'; | | char jobvsl = 'V'; |
| char jobvsr = 'V'; | | char jobvsr = 'V'; |
| char eigsort = 'N'; | | char eigsort = 'N'; |
| void* selctg = 0; | | void* selctg = 0; |
| | |
| skipping to change at line 6196 | | skipping to change at line 5865 |
|---|
| { | | { |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; |
| | |
| A = X_expr.get_ref(); | | A = X_expr.get_ref(); |
| B = Y_expr.get_ref(); | | B = Y_expr.get_ref(); |
| | |
| arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "q
z(): given matrices must be square sized" ); | | arma_debug_check( ((A.is_square() == false) || (B.is_square() == false)), "q
z(): given matrices must be square sized" ); |
| | |
| arma_debug_check( (A.n_rows != B.n_rows), "qz(): given matrices must have th
e same size" ); | | arma_debug_check( (A.n_rows != B.n_rows), "qz(): given matrices must have th
e same size" ); |
| | |
|
| if(A.is_empty()) | | if(A.is_empty()) { A.reset(); B.reset(); vsl.reset(); vsr.reset(); return t |
| { | | rue; } |
| A.reset(); | | |
| B.reset(); | | if(arma_config::check_nonfinite && A.has_nonfinite()) { return false; } |
| vsl.reset(); | | if(arma_config::check_nonfinite && B.has_nonfinite()) { return false; } |
| vsr.reset(); | | |
| return true; | | |
| } | | |
| | |
| arma_debug_assert_blas_size(A); | | arma_debug_assert_blas_size(A); |
| | |
| vsl.set_size(A.n_rows, A.n_rows); | | vsl.set_size(A.n_rows, A.n_rows); |
| vsr.set_size(A.n_rows, A.n_rows); | | vsr.set_size(A.n_rows, A.n_rows); |
| | |
| char jobvsl = 'V'; | | char jobvsl = 'V'; |
| char jobvsr = 'V'; | | char jobvsr = 'V'; |
| char eigsort = 'N'; | | char eigsort = 'N'; |
| void* selctg = 0; | | void* selctg = 0; |
| | | |
| End of changes. 189 change blocks. |
|---|
| 727 lines changed or deleted | | 413 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 