"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/op_pinv_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.

op_pinv_meat.hpp  (armadillo-10.8.2.tar.xz):op_pinv_meat.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 60 skipping to change at line 60
typedef typename T1::pod_type T; typedef typename T1::pod_type T;
arma_debug_check((tol < T(0)), "pinv(): tolerance must be >= 0"); arma_debug_check((tol < T(0)), "pinv(): tolerance must be >= 0");
// method_id = 0 -> default setting // method_id = 0 -> default setting
// method_id = 1 -> use standard algorithm // method_id = 1 -> use standard algorithm
// method_id = 2 -> use divide and conquer algorithm // method_id = 2 -> use divide and conquer algorithm
Mat<eT> A(expr.get_ref()); Mat<eT> A(expr.get_ref());
const uword n_rows = A.n_rows; if(A.is_empty()) { out.set_size(A.n_cols,A.n_rows); return true; }
const uword n_cols = A.n_cols;
if(A.is_empty()) { out.set_size(n_cols,n_rows); return true; }
if(is_op_diagmat<T1>::value || A.is_diagmat()) if(is_op_diagmat<T1>::value || A.is_diagmat())
{ {
arma_extra_debug_print("op_pinv: detected diagonal matrix"); arma_extra_debug_print("op_pinv: detected diagonal matrix");
return op_pinv::apply_diag(out, A, tol); return op_pinv::apply_diag(out, A, tol);
} }
#if defined(ARMA_OPTIMISE_SYMPD) bool do_sym = false;
bool do_sym = false; bool do_sympd = false;
bool do_sympd = false;
const bool is_sym_size_ok = (n_rows > (is_cx<eT>::yes ? uword(20) : uword(40 const bool is_sym_size_ok = (A.n_rows > (is_cx<eT>::yes ? uword(20) : uword(40
))); )));
const bool is_arg_default = ((tol == T(0)) && (method_id == uword(0))); const bool is_arg_default = ((tol == T(0)) && (method_id == uword(0)));
if( (auxlib::crippled_lapack(A) == false) && (is_arg_default || is_sym_size_ if( (arma_config::optimise_sympd) && (auxlib::crippled_lapack(A) == false) &&
ok) ) (is_arg_default || is_sym_size_ok) )
{ {
bool is_approx_sym = false; bool is_approx_sym = false;
bool is_approx_sympd = false; bool is_approx_sympd = false;
sympd_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A); sympd_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A);
do_sym = is_sym_size_ok && ((is_cx<eT>::no) ? (is_approx_sym) : (is_appr do_sym = is_sym_size_ok && ((is_cx<eT>::no) ? (is_approx_sym) : (is_approx
ox_sym && is_approx_sympd)); _sym && is_approx_sympd));
do_sympd = is_arg_default && is_approx_sympd; do_sympd = is_arg_default && is_approx_sympd;
} }
#else
const bool do_sym = false;
const bool do_sympd = false;
#endif
if(do_sympd) if(do_sympd)
{ {
arma_extra_debug_print("op_pinv: attempting sympd optimisation"); arma_extra_debug_print("op_pinv: attempting sympd optimisation");
out = A; out = A;
const T rcond_threshold = T((std::max)(uword(100), uword(A.n_rows))) * std:: bool is_sympd_junk = false;
numeric_limits<T>::epsilon(); T rcond_calc = T(0);
const T rcond_threshold = T((std::max)(uword(100), uword(A.n_rows))) * st
d::numeric_limits<T>::epsilon();
const bool status = auxlib::inv_sympd_rcond(out, rcond_threshold); const bool status = auxlib::inv_sympd_rcond(out, is_sympd_junk, rcond_calc, rcond_threshold);
if(status) { return true; } if(status && arma_isfinite(rcond_calc)) { return true; }
arma_extra_debug_print("op_pinv: sympd optimisation failed"); arma_extra_debug_print("op_pinv: sympd optimisation failed");
// auxlib::inv_sympd_rcond() will fail if A isn't really positive definite o r its rcond is below rcond_threshold // auxlib::inv_sympd_rcond() will fail if A isn't really positive definite o r its rcond is below rcond_threshold
} }
if(do_sym) if(do_sym)
{ {
arma_extra_debug_print("op_pinv: symmetric/hermitian optimisation"); arma_extra_debug_print("op_pinv: symmetric/hermitian optimisation");
return op_pinv::apply_sym(out, A, tol, method_id); return op_pinv::apply_sym(out, A, tol, method_id);
} }
// economical SVD decomposition return op_pinv::apply_gen(out, A, tol, method_id);
Mat<eT> U;
Col< T> s;
Mat<eT> V;
if(n_cols > n_rows) { A = trans(A); }
const bool status = ((method_id == uword(0)) || (method_id == uword(2))) ? aux
lib::svd_dc_econ(U, s, V, A) : auxlib::svd_econ(U, s, V, A, 'b');
if(status == false) { return false; }
// set tolerance to default if it hasn't been specified
if( (tol == T(0)) && (s.n_elem > 0) ) { tol = (std::max)(n_rows, n_cols) * s[
0] * std::numeric_limits<T>::epsilon(); }
uword count = 0;
for(uword i=0; i < s.n_elem; ++i) { count += (s[i] >= tol) ? uword(1) : uword
(0); }
if(count == 0) { out.zeros(n_cols, n_rows); return true; }
Col<T> s2(count, arma_nozeros_indicator());
uword count2 = 0;
for(uword i=0; i < s.n_elem; ++i)
{
const T val = s[i];
if(val >= tol) { s2[count2] = (val > T(0)) ? T(T(1) / val) : T(0); ++count2
; }
}
const Mat<eT> U_use(U.memptr(), U.n_rows, count, false);
const Mat<eT> V_use(V.memptr(), V.n_rows, count, false);
Mat<eT> tmp;
if(n_rows >= n_cols)
{
// out = ( (V.n_cols > count) ? V.cols(0,count-1) : V ) * diagmat(s2) * tran
s( (U.n_cols > count) ? U.cols(0,count-1) : U );
tmp = V_use * diagmat(s2);
out = tmp * trans(U_use);
}
else
{
// out = ( (U.n_cols > count) ? U.cols(0,count-1) : U ) * diagmat(s2) * tran
s( (V.n_cols > count) ? V.cols(0,count-1) : V );
tmp = U_use * diagmat(s2);
out = tmp * trans(V_use);
}
return true;
} }
template<typename eT> template<typename eT>
inline inline
bool bool
op_pinv::apply_diag(Mat<eT>& out, const Mat<eT>& A, typename get_pod_type<eT>::r esult tol) op_pinv::apply_diag(Mat<eT>& out, const Mat<eT>& A, typename get_pod_type<eT>::r esult tol)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T; typedef typename get_pod_type<eT>::result T;
skipping to change at line 271 skipping to change at line 212
if(abs_val >= tol) { eigval2[count2] = (val != T(0)) ? T(T(1) / val) : T(0) ; ++count2; } if(abs_val >= tol) { eigval2[count2] = (val != T(0)) ? T(T(1) / val) : T(0) ; ++count2; }
} }
const Mat<eT> eigvec_use(eigvec.memptr(), eigvec.n_rows, count, false); const Mat<eT> eigvec_use(eigvec.memptr(), eigvec.n_rows, count, false);
out = (eigvec_use * diagmat(eigval2)).eval() * eigvec_use.t(); out = (eigvec_use * diagmat(eigval2)).eval() * eigvec_use.t();
return true; return true;
} }
template<typename eT>
inline
bool
op_pinv::apply_gen(Mat<eT>& out, Mat<eT>& A, typename get_pod_type<eT>::result t
ol, const uword method_id)
{
arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T;
const uword n_rows = A.n_rows;
const uword n_cols = A.n_cols;
// economical SVD decomposition
Mat<eT> U;
Col< T> s;
Mat<eT> V;
if(n_cols > n_rows) { A = trans(A); }
const bool status = ((method_id == uword(0)) || (method_id == uword(2))) ? aux
lib::svd_dc_econ(U, s, V, A) : auxlib::svd_econ(U, s, V, A, 'b');
if(status == false) { return false; }
// set tolerance to default if it hasn't been specified
if( (tol == T(0)) && (s.n_elem > 0) ) { tol = (std::max)(n_rows, n_cols) * s[
0] * std::numeric_limits<T>::epsilon(); }
uword count = 0;
for(uword i=0; i < s.n_elem; ++i) { count += (s[i] >= tol) ? uword(1) : uword
(0); }
if(count == 0) { out.zeros(n_cols, n_rows); return true; }
Col<T> s2(count, arma_nozeros_indicator());
uword count2 = 0;
for(uword i=0; i < s.n_elem; ++i)
{
const T val = s[i];
if(val >= tol) { s2[count2] = (val > T(0)) ? T(T(1) / val) : T(0); ++count2
; }
}
const Mat<eT> U_use(U.memptr(), U.n_rows, count, false);
const Mat<eT> V_use(V.memptr(), V.n_rows, count, false);
Mat<eT> tmp;
if(n_rows >= n_cols)
{
// out = ( (V.n_cols > count) ? V.cols(0,count-1) : V ) * diagmat(s2) * tran
s( (U.n_cols > count) ? U.cols(0,count-1) : U );
tmp = V_use * diagmat(s2);
out = tmp * trans(U_use);
}
else
{
// out = ( (U.n_cols > count) ? U.cols(0,count-1) : U ) * diagmat(s2) * tran
s( (V.n_cols > count) ? V.cols(0,count-1) : V );
tmp = U_use * diagmat(s2);
out = tmp * trans(V_use);
}
return true;
}
//! @} //! @}
 End of changes. 11 change blocks. 
88 lines changed or deleted 98 lines changed or added

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