"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.

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

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