"Fossies" - the Fresh Open Source Software Archive  

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

[ To the main armadillo source changes report ]

sp_auxlib_meat.hpp  (armadillo-10.8.2.tar.xz):sp_auxlib_meat.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 72 skipping to change at line 72
const unwrap_spmat<T1> U(X.get_ref()); const unwrap_spmat<T1> U(X.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must b e square sized" );
if((arma_config::debug) && (sp_auxlib::rudimentary_sym_check(U.M) == false)) if((arma_config::debug) && (sp_auxlib::rudimentary_sym_check(U.M) == false))
{ {
if(is_cx<eT>::no ) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not symmetric"); } if(is_cx<eT>::no ) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not symmetric"); }
if(is_cx<eT>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); } if(is_cx<eT>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); }
} }
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_sym(): detected non-finite elements");
return false;
}
// TODO: investigate optional redirection of "sm" to ARPACK as it's capable of shift-invert; // TODO: investigate optional redirection of "sm" to ARPACK as it's capable of shift-invert;
// TODO: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0) // TODO: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0)
#if defined(ARMA_USE_NEWARP) #if defined(ARMA_USE_NEWARP)
{ {
return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts); return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts);
} }
#elif defined(ARMA_USE_ARPACK) #elif defined(ARMA_USE_ARPACK)
{ {
constexpr eT sigma = eT(0); constexpr eT sigma = eT(0);
skipping to change at line 117 skipping to change at line 123
const unwrap_spmat<T1> U(X.get_ref()); const unwrap_spmat<T1> U(X.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must b e square sized" );
if((arma_config::debug) && (sp_auxlib::rudimentary_sym_check(U.M) == false)) if((arma_config::debug) && (sp_auxlib::rudimentary_sym_check(U.M) == false))
{ {
if(is_cx<eT>::no ) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not symmetric"); } if(is_cx<eT>::no ) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not symmetric"); }
if(is_cx<eT>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); } if(is_cx<eT>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); }
} }
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_sym(): detected non-finite elements");
return false;
}
#if (defined(ARMA_USE_NEWARP) && defined(ARMA_USE_SUPERLU)) #if (defined(ARMA_USE_NEWARP) && defined(ARMA_USE_SUPERLU))
{ {
return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, sigma, opt s); return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, sigma, opt s);
} }
#elif (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) #elif (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU))
{ {
constexpr form_type form_val = form_sigma; constexpr form_type form_val = form_sigma;
return sp_auxlib::eigs_sym_arpack<eT,true>(eigval, eigvec, U.M, n_eigvals, f orm_val, sigma, opts); return sp_auxlib::eigs_sym_arpack<eT,true>(eigval, eigvec, U.M, n_eigvals, f orm_val, sigma, opts);
} }
skipping to change at line 473 skipping to change at line 485
// The process has converged, and now we need to recover the actual eigenvec tors using seupd() // The process has converged, and now we need to recover the actual eigenvec tors using seupd()
blas_int rvec = 1; // .TRUE blas_int rvec = 1; // .TRUE
blas_int nev = blas_int(n_eigvals); blas_int nev = blas_int(n_eigvals);
char howmny = 'A'; char howmny = 'A';
char bmat = 'I'; // We are considering the standard eigenvalue problem. char bmat = 'I'; // We are considering the standard eigenvalue problem.
podarray<blas_int> select(ncv); // Logical array of dimension NCV. podarray<blas_int> select(ncv); // Logical array of dimension NCV.
blas_int ldz = n; blas_int ldz = n;
select.zeros();
// seupd() will output directly into the eigval and eigvec objects. // seupd() will output directly into the eigval and eigvec objects.
eigval.zeros(n_eigvals); eigval.zeros(n_eigvals);
eigvec.zeros(n, n_eigvals); eigvec.zeros(n, n_eigvals);
arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.mempt r(), &ldz, (eT*) &sigma, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.m emptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, &info); arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.mempt r(), &ldz, (eT*) &sigma, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.m emptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, &info);
// Check for errors. // Check for errors.
if(info != 0) { arma_debug_warn_level(1, "eigs_sym(): ARPACK error ", info, " in seupd()"); return false; } if(info != 0) { arma_debug_warn_level(1, "eigs_sym(): ARPACK error ", info, " in seupd()"); return false; }
return (info == 0); return (info == 0);
skipping to change at line 511 skipping to change at line 525
inline inline
bool bool
sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase<T, T1>& X, const uword n_eigvals, const form_type form_val, con st eigs_opts& opts) sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase<T, T1>& X, const uword n_eigvals, const form_type form_val, con st eigs_opts& opts)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
const unwrap_spmat<T1> U(X.get_ref()); const unwrap_spmat<T1> U(X.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" );
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements");
return false;
}
// TODO: investigate optional redirection of "sm" to ARPACK as it's capable of shift-invert; // TODO: investigate optional redirection of "sm" to ARPACK as it's capable of shift-invert;
// TODO: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0) // TODO: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0)
#if defined(ARMA_USE_NEWARP) #if defined(ARMA_USE_NEWARP)
{ {
return sp_auxlib::eigs_gen_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts); return sp_auxlib::eigs_gen_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts);
} }
#elif defined(ARMA_USE_ARPACK) #elif defined(ARMA_USE_ARPACK)
{ {
constexpr std::complex<T> sigma = T(0); constexpr std::complex<T> sigma = T(0);
skipping to change at line 550 skipping to change at line 570
inline inline
bool bool
sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase<T, T1>& X, const uword n_eigvals, const std::complex<T> sigma, const eigs_opts& opts) sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase<T, T1>& X, const uword n_eigvals, const std::complex<T> sigma, const eigs_opts& opts)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
const unwrap_spmat<T1> U(X.get_ref()); const unwrap_spmat<T1> U(X.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" );
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements");
return false;
}
#if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) #if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU))
{ {
constexpr form_type form_val = form_sigma; constexpr form_type form_val = form_sigma;
return sp_auxlib::eigs_gen_arpack<T,true>(eigval, eigvec, U.M, n_eigvals, fo rm_val, sigma, opts); return sp_auxlib::eigs_gen_arpack<T,true>(eigval, eigvec, U.M, n_eigvals, fo rm_val, sigma, opts);
} }
#else #else
{ {
arma_ignore(eigval); arma_ignore(eigval);
arma_ignore(eigvec); arma_ignore(eigvec);
skipping to change at line 840 skipping to change at line 866
char howmny = 'A'; char howmny = 'A';
char bmat = 'I'; // We are considering the standard eigenvalue problem. char bmat = 'I'; // We are considering the standard eigenvalue problem.
podarray<blas_int> select(ncv); // logical array of dimension NCV podarray<blas_int> select(ncv); // logical array of dimension NCV
podarray<T> dr(nev + 1); // real array of dimension NEV + 1 podarray<T> dr(nev + 1); // real array of dimension NEV + 1
podarray<T> di(nev + 1); // real array of dimension NEV + 1 podarray<T> di(nev + 1); // real array of dimension NEV + 1
podarray<T> z(n * (nev + 1)); // real N by NEV array if HOWMNY = 'A' podarray<T> z(n * (nev + 1)); // real N by NEV array if HOWMNY = 'A'
blas_int ldz = n; blas_int ldz = n;
podarray<T> workev(3 * ncv); podarray<T> workev(3 * ncv);
select.zeros();
dr.zeros(); dr.zeros();
di.zeros(); di.zeros();
z.zeros(); z.zeros();
workev.zeros();
arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr(), z.m emptr(), &ldz, (T*) &sigmar, (T*) &sigmai, workev.memptr(), &bmat, &n, which, &n ev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr( ), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr(), z.m emptr(), &ldz, (T*) &sigmar, (T*) &sigmai, workev.memptr(), &bmat, &n, which, &n ev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr( ), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info);
// Check for errors. // Check for errors.
if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; } if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; }
// Put it into the outputs. // Put it into the outputs.
eigval.set_size(n_eigvals); eigval.set_size(n_eigvals);
eigvec.zeros(n, n_eigvals); eigvec.zeros(n, n_eigvals);
skipping to change at line 920 skipping to change at line 948
inline inline
bool bool
sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase< std::complex<T>, T1>& X_expr, const uword n_eigvals, const for m_type form_val, const eigs_opts& opts) sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase< std::complex<T>, T1>& X_expr, const uword n_eigvals, const for m_type form_val, const eigs_opts& opts)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
const unwrap_spmat<T1> U(X_expr.get_ref()); const unwrap_spmat<T1> U(X_expr.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" );
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements");
return false;
}
constexpr std::complex<T> sigma = T(0); constexpr std::complex<T> sigma = T(0);
return sp_auxlib::eigs_gen<T, false>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); return sp_auxlib::eigs_gen<T, false>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts);
} }
//! immediate eigendecomposition of non-symmetric complex sparse object //! immediate eigendecomposition of non-symmetric complex sparse object
template<typename T, typename T1> template<typename T, typename T1>
inline inline
bool bool
sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase< std::complex<T>, T1>& X, const uword n_eigvals, const std::com plex<T> sigma, const eigs_opts& opts) sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigv ec, const SpBase< std::complex<T>, T1>& X, const uword n_eigvals, const std::com plex<T> sigma, const eigs_opts& opts)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
const unwrap_spmat<T1> U(X.get_ref()); const unwrap_spmat<T1> U(X.get_ref());
arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" ); arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must b e square sized" );
if(arma_config::check_nonfinite && U.M.has_nonfinite())
{
arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements");
return false;
}
#if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) #if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU))
{ {
constexpr form_type form_val = form_sigma; constexpr form_type form_val = form_sigma;
return sp_auxlib::eigs_gen<T, true>(eigval, eigvec, U.M, n_eigvals, form_val , sigma, opts); return sp_auxlib::eigs_gen<T, true>(eigval, eigvec, U.M, n_eigvals, form_val , sigma, opts);
} }
#else #else
{ {
arma_ignore(eigval); arma_ignore(eigval);
arma_ignore(eigvec); arma_ignore(eigvec);
skipping to change at line 1066 skipping to change at line 1106
char howmny = 'A'; char howmny = 'A';
char bmat = 'I'; // We are considering the standard eigenvalue problem. char bmat = 'I'; // We are considering the standard eigenvalue problem.
podarray<blas_int> select(ncv); // logical array of dimension NCV podarray<blas_int> select(ncv); // logical array of dimension NCV
podarray<std::complex<T>> d(nev + 1); // complex array of dimension NEV + 1 podarray<std::complex<T>> d(nev + 1); // complex array of dimension NEV + 1
podarray<std::complex<T>> z(n * nev); // complex N by NEV array if HOWMNY = 'A' podarray<std::complex<T>> z(n * nev); // complex N by NEV array if HOWMNY = 'A'
blas_int ldz = n; blas_int ldz = n;
podarray<std::complex<T>> workev(2 * ncv); podarray<std::complex<T>> workev(2 * ncv);
select.zeros();
d.zeros();
z.zeros();
workev.zeros();
// Prepare the outputs; neupd() will write directly to them. // Prepare the outputs; neupd() will write directly to them.
eigval.zeros(n_eigvals); eigval.zeros(n_eigvals);
eigvec.zeros(n, n_eigvals); eigvec.zeros(n, n_eigvals);
arpack::neupd(&rvec, &howmny, select.memptr(), eigval.memptr(), arpack::neupd(&rvec, &howmny, select.memptr(), eigval.memptr(),
(std::complex<T>*) NULL, eigvec.memptr(), &ldz, (std::complex<T>*) &sigma, (std: :complex<T>*) NULL, workev.memptr(), &bmat, &n, which, &nev, &tol, resid.memptr( ), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), work l.memptr(), &lworkl, rwork.memptr(), &info); (std::complex<T>*) NULL, eigvec.memptr(), &ldz, (std::complex<T>*) &sigma, (std: :complex<T>*) NULL, workev.memptr(), &bmat, &n, which, &nev, &tol, resid.memptr( ), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), work l.memptr(), &lworkl, rwork.memptr(), &info);
// Check for errors. // Check for errors.
if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; } if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; }
skipping to change at line 1137 skipping to change at line 1182
arma_debug_check( (A.n_rows != X.n_rows), "spsolve(): number of rows in the given objects must be the same" ); arma_debug_check( (A.n_rows != X.n_rows), "spsolve(): number of rows in the given objects must be the same" );
if(A.is_empty() || X.is_empty()) if(A.is_empty() || X.is_empty())
{ {
X.zeros(A.n_cols, X.n_cols); X.zeros(A.n_cols, X.n_cols);
return true; return true;
} }
if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; } if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; }
if(arma_config::check_nonfinite && (A.has_nonfinite() || X.has_nonfinite()))
{
arma_debug_warn_level(3, "spsolve(): detected non-finite elements");
return false;
}
if(arma_config::debug) if(arma_config::debug)
{ {
bool overflow = false; bool overflow = false;
overflow = (A.n_nonzero > INT_MAX); overflow = (A.n_nonzero > INT_MAX);
overflow = (A.n_rows > INT_MAX) || overflow; overflow = (A.n_rows > INT_MAX) || overflow;
overflow = (A.n_cols > INT_MAX) || overflow; overflow = (A.n_cols > INT_MAX) || overflow;
overflow = (X.n_rows > INT_MAX) || overflow; overflow = (X.n_rows > INT_MAX) || overflow;
overflow = (X.n_cols > INT_MAX) || overflow; overflow = (X.n_cols > INT_MAX) || overflow;
skipping to change at line 1256 skipping to change at line 1307
{ {
arma_stop_logic_error("spsolve(): solving under-determined systems current ly not supported"); arma_stop_logic_error("spsolve(): solving under-determined systems current ly not supported");
X.soft_reset(); X.soft_reset();
return false; return false;
} }
arma_debug_check( (A.n_rows != B.n_rows), "spsolve(): number of rows in the given objects must be the same" ); arma_debug_check( (A.n_rows != B.n_rows), "spsolve(): number of rows in the given objects must be the same" );
X.zeros(A.n_cols, B.n_cols); // set the elements to zero, as we don't trust the SuperLU spaghetti code X.zeros(A.n_cols, B.n_cols); // set the elements to zero, as we don't trust the SuperLU spaghetti code
if(A.is_empty() || B.is_empty()) if(A.is_empty() || B.is_empty()) { return true; }
{
return true;
}
if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; } if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; }
if(arma_config::check_nonfinite && (A.has_nonfinite() || X.has_nonfinite()))
{
arma_debug_warn_level(3, "spsolve(): detected non-finite elements");
return false;
}
if(arma_config::debug) if(arma_config::debug)
{ {
bool overflow; bool overflow;
overflow = (A.n_nonzero > INT_MAX); overflow = (A.n_nonzero > INT_MAX);
overflow = (A.n_rows > INT_MAX) || overflow; overflow = (A.n_rows > INT_MAX) || overflow;
overflow = (A.n_cols > INT_MAX) || overflow; overflow = (A.n_cols > INT_MAX) || overflow;
overflow = (B.n_rows > INT_MAX) || overflow; overflow = (B.n_rows > INT_MAX) || overflow;
overflow = (B.n_cols > INT_MAX) || overflow; overflow = (B.n_cols > INT_MAX) || overflow;
overflow = (X.n_rows > INT_MAX) || overflow; overflow = (X.n_rows > INT_MAX) || overflow;
skipping to change at line 1314 skipping to change at line 1368
superlu_array_wrangler<T> berr(B.n_cols+1); superlu_array_wrangler<T> berr(B.n_cols+1);
superlu::GlobalLU_t glu; superlu::GlobalLU_t glu;
arrayops::inplace_set(reinterpret_cast<char*>(&glu), char(0), sizeof(superlu ::GlobalLU_t)); arrayops::inplace_set(reinterpret_cast<char*>(&glu), char(0), sizeof(superlu ::GlobalLU_t));
superlu::mem_usage_t mu; superlu::mem_usage_t mu;
arrayops::inplace_set(reinterpret_cast<char*>(&mu), char(0), sizeof(superlu: :mem_usage_t)); arrayops::inplace_set(reinterpret_cast<char*>(&mu), char(0), sizeof(superlu: :mem_usage_t));
superlu_stat_wrangler stat; superlu_stat_wrangler stat;
char equed[8]; // extra characters for paranoia char equed[8] = {}; // extra characters for paranoia
T rpg = T(0); T rpg = T(0);
T rcond = T(0); T rcond = T(0);
int info = int(0); // Return code. int info = int(0); // Return code.
char work[8]; char work[8] = {};
int lwork = int(0); // 0 means superlu will allocate memory int lwork = int(0); // 0 means superlu will allocate memory
arma_extra_debug_print("superlu::gssvx()"); arma_extra_debug_print("superlu::gssvx()");
superlu::gssvx<eT>(&options, a.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr() , etree.get_ptr(), equed, R.get_ptr(), C.get_ptr(), l.get_ptr(), u.get_ptr(), &w ork[0], lwork, b.get_ptr(), x.get_ptr(), &rpg, &rcond, ferr.get_ptr(), berr.get_ ptr(), &glu, &mu, stat.get_ptr(), &info); superlu::gssvx<eT>(&options, a.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr() , etree.get_ptr(), equed, R.get_ptr(), C.get_ptr(), l.get_ptr(), u.get_ptr(), &w ork[0], lwork, b.get_ptr(), x.get_ptr(), &rpg, &rcond, ferr.get_ptr(), berr.get_ ptr(), &glu, &mu, stat.get_ptr(), &info);
bool status = false; bool status = false;
// Process the return code. // Process the return code.
if(info == 0) if(info == 0)
{ {
status = true; status = true;
skipping to change at line 1903 skipping to change at line 1957
// basically means that we call saupd()/naupd() and it tells us with some // basically means that we call saupd()/naupd() and it tells us with some
// return code what we need to do next (usually a matrix-vector product) and // return code what we need to do next (usually a matrix-vector product) and
// then call it again. So this results in some type of iterative process // then call it again. So this results in some type of iterative process
// where we call saupd()/naupd() many times. // where we call saupd()/naupd() many times.
blas_int ido = 0; // This must be 0 for the first call. blas_int ido = 0; // This must be 0 for the first call.
char bmat = 'I'; // We are considering the standard eigenvalue problem. char bmat = 'I'; // We are considering the standard eigenvalue problem.
n = X.n_rows; // The size of the matrix (should already be set outside). n = X.n_rows; // The size of the matrix (should already be set outside).
blas_int nev = n_eigvals; blas_int nev = n_eigvals;
resid.set_size(n); resid.zeros(n);
// Two contraints on NCV: (NCV > NEV) for sym problems or // Two contraints on NCV: (NCV > NEV) for sym problems or
// (NCV > NEV + 2) for gen problems and (NCV <= N) // (NCV > NEV + 2) for gen problems and (NCV <= N)
// //
// We're calling either arpack::saupd() or arpack::naupd(), // We're calling either arpack::saupd() or arpack::naupd(),
// which have slighly different minimum constraint and recommended value for NCV: // which have slighly different minimum constraint and recommended value for NCV:
// http://www.caam.rice.edu/software/ARPACK/UG/node136.html // http://www.caam.rice.edu/software/ARPACK/UG/node136.html
// http://www.caam.rice.edu/software/ARPACK/UG/node138.html // http://www.caam.rice.edu/software/ARPACK/UG/node138.html
if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); } if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); }
if(ncv > n ) { ncv = n; } if(ncv > n ) { ncv = n; }
v.set_size(n * ncv); // Array N by NCV (output). v.zeros(n * ncv); // Array N by NCV (output).
rwork.set_size(ncv); // Work array of size NCV for complex calls. rwork.zeros(ncv); // Work array of size NCV for complex calls.
ldv = n; // "Leading dimension of V exactly as declared in the calling progr am." ldv = n; // "Leading dimension of V exactly as declared in the calling progr am."
// IPARAM: integer array of length 11. // IPARAM: integer array of length 11.
iparam.zeros(11); iparam.zeros(11);
iparam(0) = 1; // Exact shifts (not provided by us). iparam(0) = 1; // Exact shifts (not provided by us).
iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but th ey were written in the ancient times. iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but th ey were written in the ancient times.
iparam(6) = 1; // Mode 1: A * x = lambda * x. iparam(6) = 1; // Mode 1: A * x = lambda * x.
// IPNTR: integer array of length 14 (output). // IPNTR: integer array of length 14 (output).
ipntr.set_size(14); ipntr.zeros(14);
// Real work array used in the basic Arnoldi iteration for reverse communica tion. // Real work array used in the basic Arnoldi iteration for reverse communica tion.
workd.set_size(3 * n); workd.zeros(3 * n);
// lworkl must be at least 3 * NCV^2 + 6 * NCV. // lworkl must be at least 3 * NCV^2 + 6 * NCV.
lworkl = 3 * (ncv * ncv) + 6 * ncv; lworkl = 3 * (ncv * ncv) + 6 * ncv;
// Real work array of length lworkl. // Real work array of length lworkl.
workl.set_size(lworkl); workl.zeros(lworkl);
info = 0; // Set to 0 initially to use random initial vector. info = 0; // Set to 0 initially to use random initial vector.
// All the parameters have been set or created. Time to loop a lot. // All the parameters have been set or created. Time to loop a lot.
while(ido != 99) while(ido != 99)
{ {
// Call saupd() or naupd() with the current parameters. // Call saupd() or naupd() with the current parameters.
if(sym) if(sym)
{ {
arma_extra_debug_print("arpack::saupd()"); arma_extra_debug_print("arpack::saupd()");
skipping to change at line 2073 skipping to change at line 2127
{ {
char which_lm[3] = "LM"; char which_lm[3] = "LM";
char* which = which_lm; // NOTE: which_lm is the assumed operation when usi ng shift-invert char* which = which_lm; // NOTE: which_lm is the assumed operation when usi ng shift-invert
blas_int ido = 0; // This must be 0 for the first call. blas_int ido = 0; // This must be 0 for the first call.
char bmat = 'I'; // We are considering the standard eigenvalue problem. char bmat = 'I'; // We are considering the standard eigenvalue problem.
n = X.n_rows; // The size of the matrix (should already be set outside). n = X.n_rows; // The size of the matrix (should already be set outside).
blas_int nev = n_eigvals; blas_int nev = n_eigvals;
resid.set_size(n); resid.zeros(n);
// Two contraints on NCV: (NCV > NEV) for sym problems or // Two contraints on NCV: (NCV > NEV) for sym problems or
// (NCV > NEV + 2) for gen problems and (NCV <= N) // (NCV > NEV + 2) for gen problems and (NCV <= N)
// //
// We're calling either arpack::saupd() or arpack::naupd(), // We're calling either arpack::saupd() or arpack::naupd(),
// which have slighly different minimum constraint and recommended value for NCV: // which have slighly different minimum constraint and recommended value for NCV:
// http://www.caam.rice.edu/software/ARPACK/UG/node136.html // http://www.caam.rice.edu/software/ARPACK/UG/node136.html
// http://www.caam.rice.edu/software/ARPACK/UG/node138.html // http://www.caam.rice.edu/software/ARPACK/UG/node138.html
if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); } if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); }
if(ncv > n ) { ncv = n; } if(ncv > n ) { ncv = n; }
v.set_size(n * ncv); // Array N by NCV (output). v.zeros(n * ncv); // Array N by NCV (output).
rwork.set_size(ncv); // Work array of size NCV for complex calls. rwork.zeros(ncv); // Work array of size NCV for complex calls.
ldv = n; // "Leading dimension of V exactly as declared in the calling progr am." ldv = n; // "Leading dimension of V exactly as declared in the calling progr am."
// IPARAM: integer array of length 11. // IPARAM: integer array of length 11.
iparam.zeros(11); iparam.zeros(11);
iparam(0) = 1; // Exact shifts (not provided by us). iparam(0) = 1; // Exact shifts (not provided by us).
iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but th ey were written in the ancient times. iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but th ey were written in the ancient times.
// iparam(6) = 1; // Mode 1: A * x = lambda * x. // iparam(6) = 1; // Mode 1: A * x = lambda * x.
// Change IPARAM for shift-invert // Change IPARAM for shift-invert
iparam(6) = 3; // Mode 3: A * x = lambda * M * x, M symmetric semi-definite . OP = inv[A - sigma*M]*M (A complex) or Real_Part{ inv[A - sigma*M]*M } (A real) and B = M. iparam(6) = 3; // Mode 3: A * x = lambda * M * x, M symmetric semi-definite . OP = inv[A - sigma*M]*M (A complex) or Real_Part{ inv[A - sigma*M]*M } (A real) and B = M.
// IPNTR: integer array of length 14 (output). // IPNTR: integer array of length 14 (output).
ipntr.set_size(14); ipntr.zeros(14);
// Real work array used in the basic Arnoldi iteration for reverse communica tion. // Real work array used in the basic Arnoldi iteration for reverse communica tion.
workd.set_size(3 * n); workd.zeros(3 * n);
// lworkl must be at least 3 * NCV^2 + 6 * NCV. // lworkl must be at least 3 * NCV^2 + 6 * NCV.
lworkl = 3 * (ncv * ncv) + 6 * ncv; lworkl = 3 * (ncv * ncv) + 6 * ncv;
// Real work array of length lworkl. // Real work array of length lworkl.
workl.set_size(lworkl); workl.zeros(lworkl);
info = 0; // Set to 0 initially to use random initial vector. info = 0; // Set to 0 initially to use random initial vector.
superlu_opts superlu_opts_default; superlu_opts superlu_opts_default;
superlu::superlu_options_t options; superlu::superlu_options_t options;
sp_auxlib::set_superlu_opts(options, superlu_opts_default); sp_auxlib::set_superlu_opts(options, superlu_opts_default);
int lwork = 0; int lwork = 0;
superlu::trans_t trans = superlu::NOTRANS; superlu::trans_t trans = superlu::NOTRANS;
superlu::GlobalLU_t Glu; /* Not needed on return. */ superlu::GlobalLU_t Glu; /* Not needed on return. */
 End of changes. 25 change blocks. 
22 lines changed or deleted 76 lines changed or added
To the Header Files section of the armadillo source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)