"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/auxlib_meat.hpp" between
armadillo-10.8.2.tar.xz and armadillo-11.0.0.tar.xz

About: Armadillo 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

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)