diff options
Diffstat (limited to 'src/armadillo/include/armadillo_bits/sp_auxlib_meat.hpp')
| -rw-r--r-- | src/armadillo/include/armadillo_bits/sp_auxlib_meat.hpp | 2814 |
1 files changed, 2814 insertions, 0 deletions
diff --git a/src/armadillo/include/armadillo_bits/sp_auxlib_meat.hpp b/src/armadillo/include/armadillo_bits/sp_auxlib_meat.hpp new file mode 100644 index 0000000..dbfdf2d --- /dev/null +++ b/src/armadillo/include/armadillo_bits/sp_auxlib_meat.hpp @@ -0,0 +1,2814 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + +//! \addtogroup sp_auxlib +//! @{ + + +inline +sp_auxlib::form_type +sp_auxlib::interpret_form_str(const char* form_str) + { + arma_extra_debug_sigprint(); + + // the order of the 3 if statements below is important + if( form_str == nullptr ) { return form_none; } + if( form_str[0] == char(0) ) { return form_none; } + if( form_str[1] == char(0) ) { return form_none; } + + const char c1 = form_str[0]; + const char c2 = form_str[1]; + + if(c1 == 'l') + { + if(c2 == 'm') { return form_lm; } + if(c2 == 'r') { return form_lr; } + if(c2 == 'i') { return form_li; } + if(c2 == 'a') { return form_la; } + } + else + if(c1 == 's') + { + if(c2 == 'm') { return form_sm; } + if(c2 == 'r') { return form_sr; } + if(c2 == 'i') { return form_si; } + if(c2 == 'a') { return form_sa; } + } + + return form_none; + } + + + +//! immediate eigendecomposition of symmetric real sparse object +template<typename eT, typename T1> +inline +bool +sp_auxlib::eigs_sym(Col<eT>& eigval, Mat<eT>& eigvec, const SpBase<eT, T1>& X, const uword n_eigvals, const form_type form_val, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must be square sized" ); + + 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>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); } + } + + if(arma_config::check_nonfinite && U.M.internal_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: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0) + + #if defined(ARMA_USE_NEWARP) + { + return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts); + } + #elif defined(ARMA_USE_ARPACK) + { + constexpr eT sigma = eT(0); + + return sp_auxlib::eigs_sym_arpack<eT,false>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(opts); + + arma_stop_logic_error("eigs_sym(): use of NEWARP or ARPACK must be enabled"); + return false; + } + #endif + } + + + +//! immediate eigendecomposition of symmetric real sparse object +template<typename eT, typename T1> +inline +bool +sp_auxlib::eigs_sym(Col<eT>& eigval, Mat<eT>& eigvec, const SpBase<eT, T1>& X, const uword n_eigvals, const eT sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_sym(): given matrix must be square sized" ); + + 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>::yes) { arma_debug_warn_level(1, "eigs_sym(): given matrix is not hermitian"); } + } + + if(arma_config::check_nonfinite && U.M.internal_has_nonfinite()) + { + arma_debug_warn_level(3, "eigs_sym(): detected non-finite elements"); + return false; + } + + #if (defined(ARMA_USE_NEWARP) && defined(ARMA_USE_SUPERLU)) + { + return sp_auxlib::eigs_sym_newarp(eigval, eigvec, U.M, n_eigvals, sigma, opts); + } + #elif (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) + { + constexpr form_type form_val = form_sigma; + + return sp_auxlib::eigs_sym_arpack<eT,true>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(n_eigvals); + arma_ignore(sigma); + arma_ignore(opts); + + arma_stop_logic_error("eigs_sym(): use of NEWARP or ARPACK as well as SuperLU must be enabled to use 'sigma'"); + return false; + } + #endif + } + + + +template<typename eT> +inline +bool +sp_auxlib::eigs_sym_newarp(Col<eT>& eigval, Mat<eT>& eigvec, const SpMat<eT>& X, const uword n_eigvals, const form_type form_val, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_NEWARP) + { + arma_debug_check( (form_val != form_lm) && (form_val != form_sm) && (form_val != form_la) && (form_val != form_sa), "eigs_sym(): unknown form specified" ); + + if(X.is_square() == false) { return false; } + + const newarp::SparseGenMatProd<eT> op(X); + + arma_debug_check( (n_eigvals >= op.n_rows), "eigs_sym(): n_eigvals must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (op.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + uword n = op.n_rows; + + // Use max(2*k+1, 20) as default subspace dimension for the sym case; MATLAB uses max(2*k, 20), but we need to be backward-compatible. + uword ncv_default = uword( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits, otherwise cap it. + uword ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 1)) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim must be greater than k; using k+1 instead of ", opts.subdim); + ncv = uword(n_eigvals + 1); + } + else + if(opts.subdim > n) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = uword(opts.subdim); + } + } + + // Re-check that we are within the limits + if(ncv < (n_eigvals + 1)) { ncv = (n_eigvals + 1); } + if(ncv > n ) { ncv = n; } + + eT tol = (std::max)(eT(opts.tol), std::numeric_limits<eT>::epsilon()); + + uword maxiter = uword(opts.maxiter); + + // eigval.set_size(n_eigvals); + // eigvec.set_size(n, n_eigvals); + + bool status = true; + + uword nconv = 0; + + try + { + if(form_val == form_lm) + { + newarp::SymEigsSolver< eT, newarp::EigsSelect::LARGEST_MAGN, newarp::SparseGenMatProd<eT> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_sm) + { + newarp::SymEigsSolver< eT, newarp::EigsSelect::SMALLEST_MAGN, newarp::SparseGenMatProd<eT> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_la) + { + newarp::SymEigsSolver< eT, newarp::EigsSelect::LARGEST_ALGE, newarp::SparseGenMatProd<eT> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_sa) + { + newarp::SymEigsSolver< eT, newarp::EigsSelect::SMALLEST_ALGE, newarp::SparseGenMatProd<eT> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + } + catch(const std::runtime_error&) + { + status = false; + } + + if(status == true) + { + if(nconv == 0) { status = false; } + } + + return status; + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(opts); + + return false; + } + #endif + } + + + +template<typename eT> +inline +bool +sp_auxlib::eigs_sym_newarp(Col<eT>& eigval, Mat<eT>& eigvec, const SpMat<eT>& X, const uword n_eigvals, const eT sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_NEWARP) + { + if(X.is_square() == false) { return false; } + + const newarp::SparseGenRealShiftSolve<eT> op(X, sigma); + + if(op.valid == false) { return false; } + + arma_debug_check( (n_eigvals >= op.n_rows), "eigs_sym(): n_eigvals must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (op.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + uword n = op.n_rows; + + // Use max(2*k+1, 20) as default subspace dimension for the sym case; MATLAB uses max(2*k, 20), but we need to be backward-compatible. + uword ncv_default = uword( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits, otherwise cap it. + uword ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 1)) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim must be greater than k; using k+1 instead of ", opts.subdim); + ncv = uword(n_eigvals + 1); + } + else + if(opts.subdim > n) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = uword(opts.subdim); + } + } + + // Re-check that we are within the limits + if(ncv < (n_eigvals + 1)) { ncv = (n_eigvals + 1); } + if(ncv > n ) { ncv = n; } + + eT tol = (std::max)(eT(opts.tol), std::numeric_limits<eT>::epsilon()); + + uword maxiter = uword(opts.maxiter); + + // eigval.set_size(n_eigvals); + // eigvec.set_size(n, n_eigvals); + + bool status = true; + + uword nconv = 0; + + try + { + newarp::SymEigsShiftSolver< eT, newarp::EigsSelect::LARGEST_MAGN, newarp::SparseGenRealShiftSolve<eT> > eigs(op, n_eigvals, ncv, sigma); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + catch(const std::runtime_error&) + { + status = false; + } + + if(status == true) + { + if(nconv == 0) { status = false; } + } + + return status; + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(sigma); + arma_ignore(opts); + + return false; + } + #endif + } + + + +template<typename eT, bool use_sigma> +inline +bool +sp_auxlib::eigs_sym_arpack(Col<eT>& eigval, Mat<eT>& eigvec, const SpMat<eT>& X, const uword n_eigvals, const form_type form_val, const eT sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_ARPACK) + { + arma_debug_check( (form_val != form_lm) && (form_val != form_sm) && (form_val != form_la) && (form_val != form_sa) && (form_val != form_sigma), "eigs_sym(): unknown form specified" ); + + if(X.is_square() == false) { return false; } + + char which_sm[3] = "SM"; + char which_lm[3] = "LM"; + char which_sa[3] = "SA"; + char which_la[3] = "LA"; + char* which; + + switch(form_val) + { + case form_sm: which = which_sm; break; + case form_lm: which = which_lm; break; + case form_sa: which = which_sa; break; + case form_la: which = which_la; break; + + default: which = which_lm; break; + } + + // Make sure we aren't asking for every eigenvalue. + // The _saupd() functions allow asking for one more eigenvalue than the _naupd() functions. + arma_debug_check( (n_eigvals >= X.n_rows), "eigs_sym(): n_eigvals must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (X.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + // Set up variables that get used for neupd(). + blas_int n, ncv, ncv_default, ldv, lworkl, info, maxiter; + + eT tol = eT(opts.tol); + maxiter = blas_int(opts.maxiter); + + podarray<eT> resid, v, workd, workl; + podarray<blas_int> iparam, ipntr; + podarray<eT> rwork; // Not used in this case. + + n = blas_int(X.n_rows); // The size of the matrix. + + // Use max(2*k+1, 20) as default subspace dimension for the sym case; MATLAB uses max(2*k, 20), but we need to be backward-compatible. + ncv_default = blas_int( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits + ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 1)) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim must be greater than k; using k+1 instead of ", opts.subdim); + ncv = blas_int(n_eigvals + 1); + } + else + if(blas_int(opts.subdim) > n) + { + arma_debug_warn_level(1, "eigs_sym(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = blas_int(opts.subdim); + } + } + + if(use_sigma) + //if(form_val == form_sigma) + { + run_aupd_shiftinvert(n_eigvals, sigma, X, true /* sym, not gen */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + else + { + const SpMat<eT> Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, true /* sym, not gen */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + + if(info != 0) { return false; } + + // The process has converged, and now we need to recover the actual eigenvectors using seupd() + blas_int rvec = 1; // .TRUE + blas_int nev = blas_int(n_eigvals); + + char howmny = 'A'; + char bmat = 'I'; // We are considering the standard eigenvalue problem. + + podarray<blas_int> select(ncv, arma_zeros_indicator()); // Logical array of dimension NCV. + blas_int ldz = n; + + // seupd() will output directly into the eigval and eigvec objects. + eigval.zeros( n_eigvals); + eigvec.zeros(n, n_eigvals); + + arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.memptr(), &ldz, (eT*) &sigma, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, &info); + + // Check for errors. + if(info != 0) { arma_debug_warn_level(1, "eigs_sym(): ARPACK error ", info, " in seupd()"); return false; } + + return (info == 0); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(sigma); + arma_ignore(opts); + + return false; + } + #endif + } + + + +//! immediate eigendecomposition of non-symmetric real sparse object +template<typename T, typename T1> +inline +bool +sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpBase<T, T1>& X, const uword n_eigvals, const form_type form_val, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must be square sized" ); + + if(arma_config::check_nonfinite && U.M.internal_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: in shift-invert mode, "sm" maps to "lm" of the shift-inverted matrix (with sigma = 0) + + #if defined(ARMA_USE_NEWARP) + { + return sp_auxlib::eigs_gen_newarp(eigval, eigvec, U.M, n_eigvals, form_val, opts); + } + #elif defined(ARMA_USE_ARPACK) + { + constexpr std::complex<T> sigma = T(0); + + return sp_auxlib::eigs_gen_arpack<T,false>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(opts); + + arma_stop_logic_error("eigs_gen(): use of NEWARP or ARPACK must be enabled"); + return false; + } + #endif + } + + + +//! immediate eigendecomposition of non-symmetric real sparse object +template<typename T, typename T1> +inline +bool +sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpBase<T, T1>& X, const uword n_eigvals, const std::complex<T> sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must be square sized" ); + + if(arma_config::check_nonfinite && U.M.internal_has_nonfinite()) + { + arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements"); + return false; + } + + #if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) + { + constexpr form_type form_val = form_sigma; + + return sp_auxlib::eigs_gen_arpack<T,true>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(n_eigvals); + arma_ignore(sigma); + arma_ignore(opts); + + arma_stop_logic_error("eigs_gen(): use of ARPACK and SuperLU must be enabled to use 'sigma'"); + return false; + } + #endif + } + + + +template<typename T> +inline +bool +sp_auxlib::eigs_gen_newarp(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpMat<T>& X, const uword n_eigvals, const form_type form_val, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_NEWARP) + { + arma_debug_check( (form_val != form_lm) && (form_val != form_sm) && (form_val != form_lr) && (form_val != form_sr) && (form_val != form_li) && (form_val != form_si), "eigs_gen(): unknown form specified" ); + + if(X.is_square() == false) { return false; } + + const newarp::SparseGenMatProd<T> op(X); + + arma_debug_check( (n_eigvals + 1 >= op.n_rows), "eigs_gen(): n_eigvals + 1 must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (op.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + uword n = op.n_rows; + + // Use max(2*k+1, 20) as default subspace dimension for the gen case; same as MATLAB. + uword ncv_default = uword( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits + uword ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 3)) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim must be greater than k+2; using k+3 instead of ", opts.subdim); + ncv = uword(n_eigvals + 3); + } + else + if(opts.subdim > n) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = uword(opts.subdim); + } + } + + // Re-check that we are within the limits + if(ncv < (n_eigvals + 3)) { ncv = (n_eigvals + 3); } + if(ncv > n ) { ncv = n; } + + T tol = (std::max)(T(opts.tol), std::numeric_limits<T>::epsilon()); + + uword maxiter = uword(opts.maxiter); + + // eigval.set_size(n_eigvals); + // eigvec.set_size(n, n_eigvals); + + bool status = true; + + uword nconv = 0; + + try + { + if(form_val == form_lm) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::LARGEST_MAGN, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_sm) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::SMALLEST_MAGN, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_lr) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::LARGEST_REAL, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_sr) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::SMALLEST_REAL, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_li) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::LARGEST_IMAG, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + else + if(form_val == form_si) + { + newarp::GenEigsSolver< T, newarp::EigsSelect::SMALLEST_IMAG, newarp::SparseGenMatProd<T> > eigs(op, n_eigvals, ncv); + eigs.init(); + nconv = eigs.compute(maxiter, tol); + eigval = eigs.eigenvalues(); + eigvec = eigs.eigenvectors(); + } + } + catch(const std::runtime_error&) + { + status = false; + } + + if(status == true) + { + if(nconv == 0) { status = false; } + } + + return status; + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(opts); + + return false; + } + #endif + } + + + + +template<typename T, bool use_sigma> +inline +bool +sp_auxlib::eigs_gen_arpack(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpMat<T>& X, const uword n_eigvals, const form_type form_val, const std::complex<T> sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_ARPACK) + { + arma_debug_check( (form_val != form_lm) && (form_val != form_sm) && (form_val != form_lr) && (form_val != form_sr) && (form_val != form_li) && (form_val != form_si) && (form_val != form_sigma), "eigs_gen(): unknown form specified" ); + + if(X.is_square() == false) { return false; } + + char which_lm[3] = "LM"; + char which_sm[3] = "SM"; + char which_lr[3] = "LR"; + char which_sr[3] = "SR"; + char which_li[3] = "LI"; + char which_si[3] = "SI"; + + char* which; + + switch(form_val) + { + case form_lm: which = which_lm; break; + case form_sm: which = which_sm; break; + case form_lr: which = which_lr; break; + case form_sr: which = which_sr; break; + case form_li: which = which_li; break; + case form_si: which = which_si; break; + + default: which = which_lm; + } + + // Make sure we aren't asking for every eigenvalue. + arma_debug_check( (n_eigvals + 1 >= X.n_rows), "eigs_gen(): n_eigvals + 1 must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (X.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + // Set up variables that get used for neupd(). + blas_int n, ncv, ncv_default, ldv, lworkl, info, maxiter; + + T tol = T(opts.tol); + maxiter = blas_int(opts.maxiter); + + podarray<T> resid, v, workd, workl; + podarray<blas_int> iparam, ipntr; + podarray<T> rwork; // Not used in the real case. + + n = blas_int(X.n_rows); // The size of the matrix. + + // Use max(2*k+1, 20) as default subspace dimension for the gen case; same as MATLAB. + ncv_default = blas_int( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits + ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 3)) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim must be greater than k+2; using k+3 instead of ", opts.subdim); + ncv = blas_int(n_eigvals + 3); + } + else + if(blas_int(opts.subdim) > n) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = blas_int(opts.subdim); + } + } + + // WARNING!!! + // We are still not able to apply truly complex shifts to real matrices, + // in which case the OP that ARPACK wants is different (see [s/d]naupd). + // Also, if sigma contains a non-zero imaginary part, retrieving the eigenvalues + // becomes utterly messy (see [s/d]eupd, remark #3). + // We should never get to the point in which the imaginary part of sigma is non-zero; + // the user-facing functions currently convert X from real to complex if a complex sigma is detected. + // The check here is just for extra safety, and as a reminder of what's missing. + T sigmar = real(sigma); + T sigmai = imag(sigma); + + if(use_sigma) + //if(form_val == form_sigma) + { + if(sigmai != T(0)) { arma_stop_logic_error("eigs_gen(): complex 'sigma' not applicable to real matrix"); return false; } + + run_aupd_shiftinvert(n_eigvals, sigmar, X, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + else + { + const SpMat<T> Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + + if(info != 0) { return false; } + + // The process has converged, and now we need to recover the actual eigenvectors using neupd(). + blas_int rvec = 1; // .TRUE + blas_int nev = blas_int(n_eigvals); + + char howmny = 'A'; + char bmat = 'I'; // We are considering the standard eigenvalue problem. + + podarray<blas_int> select(ncv, arma_zeros_indicator()); // logical array of dimension NCV + podarray<T> dr(nev + 1, arma_zeros_indicator()); // real array of dimension NEV + 1 + podarray<T> di(nev + 1, arma_zeros_indicator()); // real array of dimension NEV + 1 + podarray<T> z(n * (nev + 1), arma_zeros_indicator()); // real N by NEV array if HOWMNY = 'A' + podarray<T> workev(3 * ncv, arma_zeros_indicator()); + + blas_int ldz = n; + + arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr(), z.memptr(), &ldz, (T*) &sigmar, (T*) &sigmai, workev.memptr(), &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); + + // Check for errors. + if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; } + + // Put it into the outputs. + eigval.set_size(n_eigvals); + eigvec.zeros(n, n_eigvals); + + for(uword i = 0; i < n_eigvals; ++i) + { + eigval[i] = std::complex<T>(dr[i], di[i]); + } + + // Now recover the eigenvectors. + for(uword i = 0; i < n_eigvals; ++i) + { + // ARPACK ?neupd lays things out kinda odd in memory; + // so does LAPACK ?geev -- see auxlib::eig_gen() + if((i < n_eigvals - 1) && (eigval[i] == std::conj(eigval[i + 1]))) + { + for(uword j = 0; j < uword(n); ++j) + { + eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1) + j]); + eigvec.at(j, i + 1) = std::complex<T>(z[n * i + j], -z[n * (i + 1) + j]); + } + ++i; // Skip the next one. + } + else + if((i == n_eigvals - 1) && (std::complex<T>(eigval[i]).imag() != 0.0)) + { + // We don't have the matched conjugate eigenvalue. + for(uword j = 0; j < uword(n); ++j) + { + eigvec.at(j, i) = std::complex<T>(z[n * i + j], z[n * (i + 1) + j]); + } + } + else + { + // The eigenvector is entirely real. + for(uword j = 0; j < uword(n); ++j) + { + eigvec.at(j, i) = std::complex<T>(z[n * i + j], T(0)); + } + } + } + + return (info == 0); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(sigma); + arma_ignore(opts); + + return false; + } + #endif + } + + + +//! immediate eigendecomposition of non-symmetric complex sparse object +template<typename T, typename T1> +inline +bool +sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpBase< std::complex<T>, T1>& X_expr, const uword n_eigvals, const form_type form_val, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X_expr.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must be square sized" ); + + if(arma_config::check_nonfinite && U.M.internal_has_nonfinite()) + { + arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements"); + return false; + } + + constexpr std::complex<T> sigma = T(0); + + 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 +template<typename T, typename T1> +inline +bool +sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpBase< std::complex<T>, T1>& X, const uword n_eigvals, const std::complex<T> sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat<T1> U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "eigs_gen(): given matrix must be square sized" ); + + if(arma_config::check_nonfinite && U.M.internal_has_nonfinite()) + { + arma_debug_warn_level(3, "eigs_gen(): detected non-finite elements"); + return false; + } + + #if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) + { + constexpr form_type form_val = form_sigma; + + return sp_auxlib::eigs_gen<T, true>(eigval, eigvec, U.M, n_eigvals, form_val, sigma, opts); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(n_eigvals); + arma_ignore(sigma); + arma_ignore(opts); + + arma_stop_logic_error("eigs_gen(): use of ARPACK and SuperLU must be enabled to use 'sigma'"); + return false; + } + #endif + } + + + +template<typename T, bool use_sigma> +inline +bool +sp_auxlib::eigs_gen(Col< std::complex<T> >& eigval, Mat< std::complex<T> >& eigvec, const SpMat< std::complex<T> >& X, const uword n_eigvals, const form_type form_val, const std::complex<T> sigma, const eigs_opts& opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_ARPACK) + { + // typedef typename std::complex<T> eT; + + arma_debug_check( (form_val != form_lm) && (form_val != form_sm) && (form_val != form_lr) && (form_val != form_sr) && (form_val != form_li) && (form_val != form_si) && (form_val != form_sigma), "eigs_gen(): unknown form specified" ); + + if(X.is_square() == false) { return false; } + + char which_lm[3] = "LM"; + char which_sm[3] = "SM"; + char which_lr[3] = "LR"; + char which_sr[3] = "SR"; + char which_li[3] = "LI"; + char which_si[3] = "SI"; + + char* which; + + switch(form_val) + { + case form_lm: which = which_lm; break; + case form_sm: which = which_sm; break; + case form_lr: which = which_lr; break; + case form_sr: which = which_sr; break; + case form_li: which = which_li; break; + case form_si: which = which_si; break; + + default: which = which_lm; + } + + // Make sure we aren't asking for every eigenvalue. + arma_debug_check( (n_eigvals + 1 >= X.n_rows), "eigs_gen(): n_eigvals + 1 must be less than the number of rows in the matrix" ); + + // If the matrix is empty, the case is trivial. + if( (X.n_cols == 0) || (n_eigvals == 0) ) // We already know n_cols == n_rows. + { + eigval.reset(); + eigvec.reset(); + return true; + } + + // Set up variables that get used for neupd(). + blas_int n, ncv, ncv_default, ldv, lworkl, info, maxiter; + + T tol = T(opts.tol); + maxiter = blas_int(opts.maxiter); + + podarray< std::complex<T> > resid, v, workd, workl; + podarray<blas_int> iparam, ipntr; + podarray<T> rwork; + + n = blas_int(X.n_rows); // The size of the matrix. + + // Use max(2*k+1, 20) as default subspace dimension for the gen case; same as MATLAB. + ncv_default = blas_int( ((2*n_eigvals+1)>(20)) ? (2*n_eigvals+1) : (20) ); + + // Use opts.subdim only if it's within the limits + ncv = ncv_default; + + if(opts.subdim != 0) + { + if(opts.subdim < (n_eigvals + 3)) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim must be greater than k+2; using k+3 instead of ", opts.subdim); + ncv = blas_int(n_eigvals + 3); + } + else + if(blas_int(opts.subdim) > n) + { + arma_debug_warn_level(1, "eigs_gen(): opts.subdim cannot be greater than n_rows; using n_rows instead of ", opts.subdim); + ncv = n; + } + else + { + ncv = blas_int(opts.subdim); + } + } + + if(use_sigma) + //if(form_val == form_sigma) + { + run_aupd_shiftinvert(n_eigvals, sigma, X, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + else + { + const SpMat< std::complex<T> > Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + } + + if(info != 0) { return false; } + + // The process has converged, and now we need to recover the actual eigenvectors using neupd(). + blas_int rvec = 1; // .TRUE + blas_int nev = blas_int(n_eigvals); + + char howmny = 'A'; + char bmat = 'I'; // We are considering the standard eigenvalue problem. + + podarray<blas_int> select(ncv, arma_zeros_indicator()); // logical array of dimension NCV + podarray<std::complex<T>> d(nev + 1, arma_zeros_indicator()); // complex array of dimension NEV + 1 + podarray<std::complex<T>> z(n * nev, arma_zeros_indicator()); // complex N by NEV array if HOWMNY = 'A' + podarray<std::complex<T>> workev(2 * ncv, arma_zeros_indicator()); + + blas_int ldz = n; + + // Prepare the outputs; neupd() will write directly to them. + eigval.zeros(n_eigvals); + eigvec.zeros(n, n_eigvals); + + 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(), workl.memptr(), &lworkl, rwork.memptr(), &info); + + // Check for errors. + if(info != 0) { arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in neupd()"); return false; } + + return (info == 0); + } + #else + { + arma_ignore(eigval); + arma_ignore(eigvec); + arma_ignore(X); + arma_ignore(n_eigvals); + arma_ignore(form_val); + arma_ignore(sigma); + arma_ignore(opts); + + arma_stop_logic_error("eigs_gen(): use of ARPACK must be enabled for decomposition of complex matrices"); + return false; + } + #endif + } + + + +template<typename T1, typename T2> +inline +bool +sp_auxlib::spsolve_simple(Mat<typename T1::elem_type>& X, const SpBase<typename T1::elem_type, T1>& A_expr, const Base<typename T1::elem_type, T2>& B_expr, const superlu_opts& user_opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_SUPERLU) + { + typedef typename T1::elem_type eT; + + superlu::superlu_options_t options; + sp_auxlib::set_superlu_opts(options, user_opts); + + const unwrap_spmat<T1> tmp1(A_expr.get_ref()); + const SpMat<eT>& A = tmp1.M; + + X = B_expr.get_ref(); // superlu::gssv() uses X as input (the B matrix) and as output (the solution) + + if(A.is_square() == false) + { + X.soft_reset(); + arma_stop_logic_error("spsolve(): solving under-determined / over-determined systems is currently not supported"); + return false; + } + + arma_debug_check( (A.n_rows != X.n_rows), "spsolve(): number of rows in the given objects must be the same", [&](){ X.soft_reset(); } ); + + if(A.is_empty() || X.is_empty()) + { + X.zeros(A.n_cols, X.n_cols); + return true; + } + + if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; } + + if(arma_config::check_nonfinite && (A.internal_has_nonfinite() || X.internal_has_nonfinite())) + { + arma_debug_warn_level(3, "spsolve(): detected non-finite elements"); + return false; + } + + if(arma_config::debug) + { + bool overflow = false; + + overflow = (A.n_nonzero > INT_MAX); + overflow = (A.n_rows > INT_MAX) || overflow; + overflow = (A.n_cols > INT_MAX) || overflow; + overflow = (X.n_rows > INT_MAX) || overflow; + overflow = (X.n_cols > INT_MAX) || overflow; + + if(overflow) + { + arma_stop_runtime_error("spsolve(): integer overflow: matrix dimensions are too large for integer type used by SuperLU"); + return false; + } + } + + superlu_supermatrix_wrangler x; + superlu_supermatrix_wrangler a; + + const bool status_x = wrap_to_supermatrix(x.get_ref(), X); + const bool status_a = copy_to_supermatrix(a.get_ref(), A); + + if( (status_x == false) || (status_a == false) ) { X.soft_reset(); return false; } + + superlu_supermatrix_wrangler l; + superlu_supermatrix_wrangler u; + + // paranoia: use SuperLU's memory allocation, in case it reallocs + + superlu_array_wrangler<int> perm_c(A.n_cols+1); // extra paranoia: increase array length by 1 + superlu_array_wrangler<int> perm_r(A.n_rows+1); + + superlu_stat_wrangler stat; + + int info = 0; // Return code. + + arma_extra_debug_print("superlu::gssv()"); + superlu::gssv<eT>(&options, a.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), l.get_ptr(), u.get_ptr(), x.get_ptr(), stat.get_ptr(), &info); + + + // Process the return code. + if( (info > 0) && (info <= int(A.n_cols)) ) + { + // std::ostringstream tmp; + // tmp << "spsolve(): could not solve system; LU factorisation completed, but detected zero in U(" << (info-1) << ',' << (info-1) << ')'; + // arma_debug_warn_level(1, tmp.str()); + } + else + if(info > int(A.n_cols)) + { + arma_debug_warn_level(1, "spsolve(): memory allocation failure"); + } + else + if(info < 0) + { + arma_debug_warn_level(1, "spsolve(): unknown SuperLU error code from gssv(): ", info); + } + + // No need to extract the data from x, since it's using the same memory as X + + return (info == 0); + } + #else + { + arma_ignore(X); + arma_ignore(A_expr); + arma_ignore(B_expr); + arma_ignore(user_opts); + arma_stop_logic_error("spsolve(): use of SuperLU must be enabled"); + return false; + } + #endif + } + + + +template<typename T1, typename T2> +inline +bool +sp_auxlib::spsolve_refine(Mat<typename T1::elem_type>& X, typename T1::pod_type& out_rcond, const SpBase<typename T1::elem_type, T1>& A_expr, const Base<typename T1::elem_type, T2>& B_expr, const superlu_opts& user_opts) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_SUPERLU) + { + typedef typename T1::pod_type T; + typedef typename T1::elem_type eT; + + superlu::superlu_options_t options; + sp_auxlib::set_superlu_opts(options, user_opts); + + const unwrap_spmat<T1> tmp1(A_expr.get_ref()); + const SpMat<eT>& A = tmp1.M; + + const unwrap<T2> tmp2(B_expr.get_ref()); + const Mat<eT>& B_unwrap = tmp2.M; + + const bool B_is_modified = ( (user_opts.equilibrate) || (&B_unwrap == &X) ); + + Mat<eT> B_copy; if(B_is_modified) { B_copy = B_unwrap; } + + const Mat<eT>& B = (B_is_modified) ? B_copy : B_unwrap; + + if(A.is_square() == false) + { + X.soft_reset(); + arma_stop_logic_error("spsolve(): solving under-determined / over-determined systems is currently not supported"); + return false; + } + + arma_debug_check( (A.n_rows != B.n_rows), "spsolve(): number of rows in the given objects must be the same", [&](){ X.soft_reset(); } ); + + 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()) { return true; } + + if(A.n_nonzero == uword(0)) { X.soft_reset(); return false; } + + if(arma_config::check_nonfinite && (A.internal_has_nonfinite() || B.internal_has_nonfinite())) + { + arma_debug_warn_level(3, "spsolve(): detected non-finite elements"); + return false; + } + + if(arma_config::debug) + { + bool overflow; + + overflow = (A.n_nonzero > INT_MAX); + overflow = (A.n_rows > INT_MAX) || overflow; + overflow = (A.n_cols > INT_MAX) || overflow; + overflow = (B.n_rows > INT_MAX) || overflow; + overflow = (B.n_cols > INT_MAX) || overflow; + overflow = (X.n_rows > INT_MAX) || overflow; + overflow = (X.n_cols > INT_MAX) || overflow; + + if(overflow) + { + arma_stop_runtime_error("spsolve(): integer overflow: matrix dimensions are too large for integer type used by SuperLU"); + return false; + } + } + + superlu_supermatrix_wrangler x; + superlu_supermatrix_wrangler a; + superlu_supermatrix_wrangler b; + + const bool status_x = wrap_to_supermatrix(x.get_ref(), X); + const bool status_a = copy_to_supermatrix(a.get_ref(), A); // NOTE: superlu::gssvx() modifies 'a' if equilibration is enabled + const bool status_b = wrap_to_supermatrix(b.get_ref(), B); // NOTE: superlu::gssvx() modifies 'b' if equilibration is enabled + + if( (status_x == false) || (status_a == false) || (status_b == false) ) { X.soft_reset(); return false; } + + superlu_supermatrix_wrangler l; + superlu_supermatrix_wrangler u; + + // paranoia: use SuperLU's memory allocation, in case it reallocs + + superlu_array_wrangler<int> perm_c(A.n_cols+1); // extra paranoia: increase array length by 1 + superlu_array_wrangler<int> perm_r(A.n_rows+1); + superlu_array_wrangler<int> etree(A.n_cols+1); + + superlu_array_wrangler<T> R(A.n_rows+1); + superlu_array_wrangler<T> C(A.n_cols+1); + superlu_array_wrangler<T> ferr(B.n_cols+1); + superlu_array_wrangler<T> berr(B.n_cols+1); + + superlu::GlobalLU_t glu; + arrayops::fill_zeros(reinterpret_cast<char*>(&glu), sizeof(superlu::GlobalLU_t)); + + superlu::mem_usage_t mu; + arrayops::fill_zeros(reinterpret_cast<char*>(&mu), sizeof(superlu::mem_usage_t)); + + superlu_stat_wrangler stat; + + char equed[8] = {}; // extra characters for paranoia + T rpg = T(0); + T rcond = T(0); + int info = int(0); // Return code. + + char work[8] = {}; + int lwork = int(0); // 0 means superlu will allocate memory + + 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(), &work[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; + + // Process the return code. + if(info == 0) + { + status = true; + } + if( (info > 0) && (info <= int(A.n_cols)) ) + { + // std::ostringstream tmp; + // tmp << "spsolve(): could not solve system; LU factorisation completed, but detected zero in U(" << (info-1) << ',' << (info-1) << ')'; + // arma_debug_warn_level(1, tmp.str()); + } + else + if( (info == int(A.n_cols+1)) && (user_opts.allow_ugly) ) + { + arma_debug_warn_level(2, "spsolve(): system is singular to working precision (rcond: ", rcond, ")"); + status = true; + } + else + if(info > int(A.n_cols+1)) + { + arma_debug_warn_level(1, "spsolve(): memory allocation failure"); + } + else + if(info < 0) + { + arma_debug_warn_level(1, "spsolve(): unknown SuperLU error code from gssvx(): ", info); + } + + // No need to extract the data from x, since it's using the same memory as X + + out_rcond = rcond; + + return status; + } + #else + { + arma_ignore(X); + arma_ignore(out_rcond); + arma_ignore(A_expr); + arma_ignore(B_expr); + arma_ignore(user_opts); + arma_stop_logic_error("spsolve(): use of SuperLU must be enabled"); + return false; + } + #endif + } + + + +#if defined(ARMA_USE_SUPERLU) + + template<typename eT> + inline + typename get_pod_type<eT>::result + sp_auxlib::norm1(superlu::SuperMatrix* A) + { + arma_extra_debug_sigprint(); + + char norm_id = '1'; + + arma_extra_debug_print("superlu::langs()"); + return superlu::langs<eT>(&norm_id, A); + } + + + + template<typename eT> + inline + typename get_pod_type<eT>::result + sp_auxlib::lu_rcond(superlu::SuperMatrix* L, superlu::SuperMatrix* U, typename get_pod_type<eT>::result norm_val) + { + arma_extra_debug_sigprint(); + + typedef typename get_pod_type<eT>::result T; + + char norm_id = '1'; + T rcond_out = T(0); + int info = int(0); + + superlu_stat_wrangler stat; + + arma_extra_debug_print("superlu::gscon()"); + superlu::gscon<eT>(&norm_id, L, U, norm_val, &rcond_out, stat.get_ptr(), &info); + + return (info == 0) ? T(rcond_out) : T(0); + } + + + + inline + void + sp_auxlib::set_superlu_opts(superlu::superlu_options_t& options, const superlu_opts& user_opts) + { + arma_extra_debug_sigprint(); + + // default options as the starting point + superlu::set_default_opts(&options); + + // our settings + options.Trans = superlu::NOTRANS; + options.ConditionNumber = superlu::YES; + + // process user_opts + + if(user_opts.equilibrate == true) { options.Equil = superlu::YES; } + if(user_opts.equilibrate == false) { options.Equil = superlu::NO; } + + if(user_opts.symmetric == true) { options.SymmetricMode = superlu::YES; } + if(user_opts.symmetric == false) { options.SymmetricMode = superlu::NO; } + + options.DiagPivotThresh = user_opts.pivot_thresh; + + if(user_opts.permutation == superlu_opts::NATURAL) { options.ColPerm = superlu::NATURAL; } + if(user_opts.permutation == superlu_opts::MMD_ATA) { options.ColPerm = superlu::MMD_ATA; } + if(user_opts.permutation == superlu_opts::MMD_AT_PLUS_A) { options.ColPerm = superlu::MMD_AT_PLUS_A; } + if(user_opts.permutation == superlu_opts::COLAMD) { options.ColPerm = superlu::COLAMD; } + + if(user_opts.refine == superlu_opts::REF_NONE) { options.IterRefine = superlu::NOREFINE; } + if(user_opts.refine == superlu_opts::REF_SINGLE) { options.IterRefine = superlu::SLU_SINGLE; } + if(user_opts.refine == superlu_opts::REF_DOUBLE) { options.IterRefine = superlu::SLU_DOUBLE; } + if(user_opts.refine == superlu_opts::REF_EXTRA) { options.IterRefine = superlu::SLU_EXTRA; } + } + + + + template<typename eT> + inline + bool + sp_auxlib::copy_to_supermatrix(superlu::SuperMatrix& out, const SpMat<eT>& A) + { + arma_extra_debug_sigprint(); + + // We store in column-major CSC. + out.Stype = superlu::SLU_NC; + + if( is_float<eT>::value) { out.Dtype = superlu::SLU_S; } + else if( is_double<eT>::value) { out.Dtype = superlu::SLU_D; } + else if( is_cx_float<eT>::value) { out.Dtype = superlu::SLU_C; } + else if(is_cx_double<eT>::value) { out.Dtype = superlu::SLU_Z; } + + out.Mtype = superlu::SLU_GE; // Just a general matrix. We don't know more now. + + // We have to actually create the object which stores the data. + // This gets cleaned by destroy_supermatrix(). + // We have to use SuperLU's problematic memory allocation routines since they are + // not guaranteed to be new and delete. See the comments in def_superlu.hpp + superlu::NCformat* nc = (superlu::NCformat*)superlu::malloc(sizeof(superlu::NCformat)); + + if(nc == nullptr) { return false; } + + A.sync(); + + nc->nnz = A.n_nonzero; + nc->nzval = (void*) superlu::malloc(sizeof(eT) * A.n_nonzero ); + nc->colptr = (superlu::int_t*)superlu::malloc(sizeof(superlu::int_t) * (A.n_cols + 1)); + nc->rowind = (superlu::int_t*)superlu::malloc(sizeof(superlu::int_t) * A.n_nonzero ); + + if( (nc->nzval == nullptr) || (nc->colptr == nullptr) || (nc->rowind == nullptr) ) { return false; } + + // Fill the matrix. + arrayops::copy((eT*) nc->nzval, A.values, A.n_nonzero); + + // // These have to be copied by hand, because the types may differ. + // for(uword i = 0; i <= A.n_cols; ++i) { nc->colptr[i] = (int_t) A.col_ptrs[i]; } + // for(uword i = 0; i < A.n_nonzero; ++i) { nc->rowind[i] = (int_t) A.row_indices[i]; } + + arrayops::convert(nc->colptr, A.col_ptrs, A.n_cols+1 ); + arrayops::convert(nc->rowind, A.row_indices, A.n_nonzero); + + out.nrow = superlu::int_t(A.n_rows); + out.ncol = superlu::int_t(A.n_cols); + out.Store = (void*) nc; + + return true; + } + + + + // memory efficient implementation of out = A - shift*I, where A is a square matrix + template<typename eT> + inline + bool + sp_auxlib::copy_to_supermatrix_with_shift(superlu::SuperMatrix& out, const SpMat<eT>& A, const eT shift) + { + arma_extra_debug_sigprint(); + + arma_debug_check( (A.is_square() == false), "sp_auxlib::copy_to_supermatrix_with_shift(): given matrix must be square sized" ); + + if(shift == eT(0)) + { + arma_extra_debug_print("sp_auxlib::copy_to_supermatrix_with_shift(): shift is zero; redirecting to sp_auxlib::copy_to_supermatrix()"); + return sp_auxlib::copy_to_supermatrix(out, A); + } + + // We store in column-major CSC. + out.Stype = superlu::SLU_NC; + + if( is_float<eT>::value) { out.Dtype = superlu::SLU_S; } + else if( is_double<eT>::value) { out.Dtype = superlu::SLU_D; } + else if( is_cx_float<eT>::value) { out.Dtype = superlu::SLU_C; } + else if(is_cx_double<eT>::value) { out.Dtype = superlu::SLU_Z; } + + out.Mtype = superlu::SLU_GE; // Just a general matrix. We don't know more now. + + // We have to actually create the object which stores the data. + // This gets cleaned by destroy_supermatrix(). + superlu::NCformat* nc = (superlu::NCformat*)superlu::malloc(sizeof(superlu::NCformat)); + + if(nc == nullptr) { return false; } + + A.sync(); + + uword n_nonzero_diag_old = 0; + uword n_nonzero_diag_new = 0; + + const uword n_search_cols = (std::min)(A.n_rows, A.n_cols); + + for(uword j=0; j < n_search_cols; ++j) + { + const uword col_offset = A.col_ptrs[j ]; + const uword next_col_offset = A.col_ptrs[j + 1]; + + const uword* start_ptr = &(A.row_indices[ col_offset]); + const uword* end_ptr = &(A.row_indices[next_col_offset]); + + const uword wanted_row = j; + + const uword* pos_ptr = std::lower_bound(start_ptr, end_ptr, wanted_row); // binary search + + if( (pos_ptr != end_ptr) && ((*pos_ptr) == wanted_row) ) + { + // element on the main diagonal is non-zero + ++n_nonzero_diag_old; + + const uword offset = uword(pos_ptr - start_ptr); + const uword index = offset + col_offset; + + const eT new_val = A.values[index] - shift; + + if(new_val != eT(0)) { ++n_nonzero_diag_new; } + } + else + { + // element on the main diagonal is zero, but sigma is non-zero, + // so the number of new non-zero elments on the diagonal is increased + ++n_nonzero_diag_new; + } + } + + const uword out_n_nonzero = A.n_nonzero - n_nonzero_diag_old + n_nonzero_diag_new; + + arma_extra_debug_print( arma_str::format("A.n_nonzero: %u") % A.n_nonzero ); + arma_extra_debug_print( arma_str::format("n_nonzero_diag_old: %u") % n_nonzero_diag_old ); + arma_extra_debug_print( arma_str::format("n_nonzero_diag_new: %u") % n_nonzero_diag_new ); + arma_extra_debug_print( arma_str::format("out_n_nonzero: %u") % out_n_nonzero ); + + nc->nnz = out_n_nonzero; + nc->nzval = (void*) superlu::malloc(sizeof(eT) * out_n_nonzero ); + nc->colptr = (superlu::int_t*)superlu::malloc(sizeof(superlu::int_t) * (A.n_cols + 1)); + nc->rowind = (superlu::int_t*)superlu::malloc(sizeof(superlu::int_t) * out_n_nonzero ); + + if( (nc->nzval == nullptr) || (nc->colptr == nullptr) || (nc->rowind == nullptr) ) { return false; } + + // fill the matrix column by column, and insert diagonal elements when necessary + + nc->colptr[0] = 0; + + eT* values_current = (eT*) nc->nzval; + superlu::int_t* rowind_current = nc->rowind; + + uword count = 0; + + for(uword j=0; j < A.n_cols; ++j) + { + const uword idx_start = A.col_ptrs[j ]; + const uword idx_end = A.col_ptrs[j + 1]; + + const eT* values_start = values_current; + + uword i = idx_start; + + // elements in the upper triangular part, excluding the main diagonal + for(; (i < idx_end) && (A.row_indices[i] < j); ++i) + { + (*values_current) = A.values[i]; + (*rowind_current) = superlu::int_t(A.row_indices[i]); + + ++values_current; + ++rowind_current; + + ++count; + } + + // elements on the main diagonal + if( (i < idx_end) && (A.row_indices[i] == j) ) + { + // A(j,j) is non-zero + + const eT new_diag_val = A.values[i] - shift; + + if(new_diag_val != eT(0)) + { + (*values_current) = new_diag_val; + (*rowind_current) = superlu::int_t(j); + + ++values_current; + ++rowind_current; + + ++count; + } + + ++i; + } + else + { + // A(j,j) is zero, so insert a new element + + if(j < n_search_cols) + { + (*values_current) = -shift; + (*rowind_current) = superlu::int_t(j); + + ++values_current; + ++rowind_current; + + ++count; + } + } + + // elements in the lower triangular part, excluding the main diagonal + for(; i < idx_end; ++i) + { + (*values_current) = A.values[i]; + (*rowind_current) = superlu::int_t(A.row_indices[i]); + + ++values_current; + ++rowind_current; + + ++count; + } + + // number of non-zero elements in the j-th column of out + const uword nnz_col = values_current - values_start; + nc->colptr[j + 1] = superlu::int_t(nc->colptr[j] + nnz_col); + } + + arma_extra_debug_print( arma_str::format("count: %u") % count ); + + arma_check( (count != out_n_nonzero), "internal error: sp_auxlib::copy_to_supermatrix_with_shift(): count != out_n_nonzero" ); + + out.nrow = superlu::int_t(A.n_rows); + out.ncol = superlu::int_t(A.n_cols); + out.Store = (void*) nc; + + return true; + } + + + +// // for debugging only +// template<typename eT> +// inline +// void +// sp_auxlib::copy_to_spmat(SpMat<eT>& out, const superlu::SuperMatrix& A) +// { +// arma_extra_debug_sigprint(); +// +// bool type_matched = false; +// +// if( is_float<eT>::value) { type_matched = (A.Dtype == superlu::SLU_S); } +// else if( is_double<eT>::value) { type_matched = (A.Dtype == superlu::SLU_D); } +// else if( is_cx_float<eT>::value) { type_matched = (A.Dtype == superlu::SLU_C); } +// else if(is_cx_double<eT>::value) { type_matched = (A.Dtype == superlu::SLU_Z); } +// +// arma_debug_check( (type_matched == false), "copy_to_spmat(): type mismatch" ); +// arma_debug_check( (A.Mtype != superlu::SLU_GE), "copy_to_spmat(): unknown layout" ); +// +// // NOTE: the l and u instances of SuperMatrix resulting from superlu::gstrf() +// // NOTE: do not have the superlu::SLU_GE layout +// +// const superlu::NCformat* nc = (const superlu::NCformat*)(A.Store); +// +// if(nc == nullptr) { out.reset(); return; } +// +// if( (nc->nzval == nullptr) || (nc->colptr == nullptr) || (nc->rowind == nullptr) ) { out.reset(); return; } +// +// const uword A_n_rows = uword(A.nrow ); +// const uword A_n_cols = uword(A.ncol ); +// const uword A_n_nonzero = uword(nc->nnz); +// +// if(A_n_nonzero == 0) { out.zeros(A_n_rows, A_n_cols); return; } +// +// out.reserve(A_n_rows, A_n_cols, A_n_nonzero); +// +// arrayops::copy(access::rwp(out.values), (const eT*)(nc->nzval), A_n_nonzero); +// +// arrayops::convert(access::rwp(out.col_ptrs), nc->colptr, A_n_cols+1 ); +// arrayops::convert(access::rwp(out.row_indices), nc->rowind, A_n_nonzero); +// +// out.remove_zeros(); // in case SuperLU has bugs and stores zeros in sparse matrices +// } + + + + template<typename eT> + inline + bool + sp_auxlib::wrap_to_supermatrix(superlu::SuperMatrix& out, const Mat<eT>& A) + { + arma_extra_debug_sigprint(); + + // NOTE: this function re-uses memory from matrix A + + // This is being stored as a dense matrix. + out.Stype = superlu::SLU_DN; + + if( is_float<eT>::value) { out.Dtype = superlu::SLU_S; } + else if( is_double<eT>::value) { out.Dtype = superlu::SLU_D; } + else if( is_cx_float<eT>::value) { out.Dtype = superlu::SLU_C; } + else if(is_cx_double<eT>::value) { out.Dtype = superlu::SLU_Z; } + + out.Mtype = superlu::SLU_GE; + + // We have to create the object that stores the data. + superlu::DNformat* dn = (superlu::DNformat*)superlu::malloc(sizeof(superlu::DNformat)); + + if(dn == nullptr) { return false; } + + dn->lda = A.n_rows; + dn->nzval = (void*) A.memptr(); // re-use memory instead of copying + + out.nrow = A.n_rows; + out.ncol = A.n_cols; + out.Store = (void*) dn; + + return true; + } + + + + inline + void + sp_auxlib::destroy_supermatrix(superlu::SuperMatrix& out) + { + arma_extra_debug_sigprint(); + + // Clean up. + if(out.Stype == superlu::SLU_NC) + { + superlu::destroy_compcol_mat(&out); + } + else + if(out.Stype == superlu::SLU_NCP) + { + superlu::destroy_compcolperm_mat(&out); + } + else + if(out.Stype == superlu::SLU_DN) + { + // superlu::destroy_dense_mat(&out); + + // since dn->nzval is set to re-use memory from a Mat object (which manages its own memory), + // we cannot simply call superlu::destroy_dense_mat(). + // Only the out.Store structure can be freed. + + superlu::DNformat* dn = (superlu::DNformat*) out.Store; + + if(dn != nullptr) { superlu::free(dn); } + } + else + if(out.Stype == superlu::SLU_SC) + { + superlu::destroy_supernode_mat(&out); + } + else + { + // Uh, crap. + + std::ostringstream tmp; + + tmp << "sp_auxlib::destroy_supermatrix(): unhandled Stype" << std::endl; + tmp << "Stype val: " << out.Stype << std::endl; + tmp << "Stype name: "; + + if(out.Stype == superlu::SLU_NC) { tmp << "SLU_NC"; } + if(out.Stype == superlu::SLU_NCP) { tmp << "SLU_NCP"; } + if(out.Stype == superlu::SLU_NR) { tmp << "SLU_NR"; } + if(out.Stype == superlu::SLU_SC) { tmp << "SLU_SC"; } + if(out.Stype == superlu::SLU_SCP) { tmp << "SLU_SCP"; } + if(out.Stype == superlu::SLU_SR) { tmp << "SLU_SR"; } + if(out.Stype == superlu::SLU_DN) { tmp << "SLU_DN"; } + if(out.Stype == superlu::SLU_NR_loc) { tmp << "SLU_NR_loc"; } + + arma_debug_warn_level(1, tmp.str()); + arma_stop_runtime_error("internal error: sp_auxlib::destroy_supermatrix()"); + } + } + +#endif + + + +template<typename eT, typename T> +inline +void +sp_auxlib::run_aupd_plain + ( + const uword n_eigvals, char* which, + const SpMat<T>& X, const SpMat<T>& Xst, const bool sym, + blas_int& n, eT& tol, blas_int& maxiter, + podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv, + podarray<blas_int>& iparam, podarray<blas_int>& ipntr, + podarray<T>& workd, podarray<T>& workl, blas_int& lworkl, podarray<eT>& rwork, + blas_int& info + ) + { + #if defined(ARMA_USE_ARPACK) + { + // ARPACK provides a "reverse communication interface" which is an + // entertainingly archaic FORTRAN software engineering technique that + // 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 + // then call it again. So this results in some type of iterative process + // where we call saupd()/naupd() many times. + + blas_int ido = 0; // This must be 0 for the first call. + char bmat = 'I'; // We are considering the standard eigenvalue problem. + n = X.n_rows; // The size of the matrix (should already be set outside). + blas_int nev = n_eigvals; + + // resid.zeros(n); + eigs_randu_filler<T> randu_filler; + randu_filler.fill(resid, n); // use deterministic starting point + + // Two contraints on NCV: (NCV > NEV) for sym problems or + // (NCV > NEV + 2) for gen problems and (NCV <= N) + // + // We're calling either arpack::saupd() or arpack::naupd(), + // 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/node138.html + + if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); } + if(ncv > n ) { ncv = n; } + + v.zeros(n * ncv); // Array N by NCV (output). + rwork.zeros(ncv); // Work array of size NCV for complex calls. + ldv = n; // "Leading dimension of V exactly as declared in the calling program." + + // IPARAM: integer array of length 11. + iparam.zeros(11); + iparam(0) = 1; // Exact shifts (not provided by us). + iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but they were written in the ancient times. + iparam(6) = 1; // Mode 1: A * x = lambda * x. + + // IPNTR: integer array of length 14 (output). + ipntr.zeros(14); + + // Real work array used in the basic Arnoldi iteration for reverse communication. + workd.zeros(3 * n); + + // lworkl must be at least 3 * NCV^2 + 6 * NCV. + lworkl = 3 * (ncv * ncv) + 6 * ncv; + + // Real work array of length lworkl. + workl.zeros(lworkl); + + // info = 0; // resid to be filled with random values by ARPACK (non-deterministic) + info = 1; // resid is already filled with random values (deterministic) + + // All the parameters have been set or created. Time to loop a lot. + while(ido != 99) + { + // Call saupd() or naupd() with the current parameters. + if(sym) + { + arma_extra_debug_print("arpack::saupd()"); + arpack::saupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, &info); + } + else + { + arma_extra_debug_print("arpack::naupd()"); + arpack::naupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); + } + + // What do we do now? + switch (ido) + { + case -1: + // fallthrough + case 1: + { + // We need to calculate the matrix-vector multiplication y = OP * x + // where x is of length n and starts at workd(ipntr(0)), and y is of + // length n and starts at workd(ipntr(1)). + + // // OLD METHOD + // + // // operator*(sp_mat, vec) doesn't properly put the result into the + // // right place so we'll just reimplement it here for now... + // + // // Set the output to point at the right memory. We have to subtract + // // one from FORTRAN pointers... + // Col<T> out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); + // // Set the input to point at the right memory. + // Col<T> in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); + // + // out.zeros(); + // + // T* out_mem = out.memptr(); + // const T* in_mem = in.memptr(); + // + // typename SpMat<T>::const_iterator X_it = X.begin(); + // + // const uword X_nnz = X.n_nonzero; + // + // for(uword count=0; count < X_nnz; ++count, ++X_it) + // { + // const eT X_it_val = (*X_it); + // const uword X_it_row = X_it.row(); + // const uword X_it_col = X_it.col(); + // + // out_mem[X_it_row] += X_it_val * in_mem[X_it_col]; + // } + // + // // No need to modify memory further since it was all done in-place. + + + // NEW METHOD + // + // both operator*(rowvec, sp_mat) and operator*(sp_mat, colvec) can now write to an existing object + + Row<T> out(workd.memptr() + ipntr(1) - 1, n, false, true); + Row<T> in(workd.memptr() + ipntr(0) - 1, n, false, true); + + out = in * Xst; + + break; + } + case 99: + // Nothing to do here, things have converged. + break; + default: + { + return; // Parent frame can look at the value of info. + } + } + } + + // The process has ended; check the return code. + if( (info != 0) && (info != 1) ) + { + // Print warnings if there was a failure. + + if(sym) + { + arma_debug_warn_level(1, "eigs_sym(): ARPACK error ", info, " in saupd()"); + } + else + { + arma_debug_warn_level(1, "eigs_gen(): ARPACK error ", info, " in naupd()"); + } + + return; // Parent frame can look at the value of info. + } + } + #else + { + arma_ignore(n_eigvals); + arma_ignore(which); + arma_ignore(X); + arma_ignore(sym); + arma_ignore(n); + arma_ignore(tol); + arma_ignore(maxiter); + arma_ignore(resid); + arma_ignore(ncv); + arma_ignore(v); + arma_ignore(ldv); + arma_ignore(iparam); + arma_ignore(ipntr); + arma_ignore(workd); + arma_ignore(workl); + arma_ignore(lworkl); + arma_ignore(rwork); + arma_ignore(info); + } + #endif + } + + + +// Here 'sigma' is 'T', but should be 'eT'. +// Applying complex shifts to real matrices is currently not directly implemented +template<typename eT, typename T> +inline +void +sp_auxlib::run_aupd_shiftinvert + ( + const uword n_eigvals, const T sigma, + const SpMat<T>& X, const bool sym, + blas_int& n, eT& tol, blas_int& maxiter, + podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv, + podarray<blas_int>& iparam, podarray<blas_int>& ipntr, + podarray<T>& workd, podarray<T>& workl, blas_int& lworkl, podarray<eT>& rwork, + blas_int& info + ) + { + // TODO: inconsistent use of type names: T can be complex while eT can be real + + #if (defined(ARMA_USE_ARPACK) && defined(ARMA_USE_SUPERLU)) + { + char which_lm[3] = "LM"; + + char* which = which_lm; // NOTE: which_lm is the assumed operation when using shift-invert + + blas_int ido = 0; // This must be 0 for the first call. + char bmat = 'I'; // We are considering the standard eigenvalue problem. + n = X.n_rows; // The size of the matrix (should already be set outside). + blas_int nev = n_eigvals; + + // resid.zeros(n); + eigs_randu_filler<T> randu_filler; + randu_filler.fill(resid, n); // use deterministic starting point + + // Two contraints on NCV: (NCV > NEV) for sym problems or + // (NCV > NEV + 2) for gen problems and (NCV <= N) + // + // We're calling either arpack::saupd() or arpack::naupd(), + // 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/node138.html + + if(ncv < (nev + (sym ? 1 : 3))) { ncv = (nev + (sym ? 1 : 3)); } + if(ncv > n ) { ncv = n; } + + v.zeros(n * ncv); // Array N by NCV (output). + rwork.zeros(ncv); // Work array of size NCV for complex calls. + ldv = n; // "Leading dimension of V exactly as declared in the calling program." + + // IPARAM: integer array of length 11. + iparam.zeros(11); + iparam(0) = 1; // Exact shifts (not provided by us). + iparam(2) = maxiter; // Maximum iterations; all the examples use 300, but they were written in the ancient times. + // iparam(6) = 1; // Mode 1: A * x = lambda * x. + + // 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. + + // IPNTR: integer array of length 14 (output). + ipntr.zeros(14); + + // Real work array used in the basic Arnoldi iteration for reverse communication. + workd.zeros(3 * n); + + // lworkl must be at least 3 * NCV^2 + 6 * NCV. + lworkl = 3 * (ncv * ncv) + 6 * ncv; + + // Real work array of length lworkl. + workl.zeros(lworkl); + + // info = 0; // resid to be filled with random values by ARPACK (non-deterministic) + info = 1; // resid is already filled with random values (deterministic) + + superlu_opts superlu_opts_default; + superlu::superlu_options_t options; + sp_auxlib::set_superlu_opts(options, superlu_opts_default); + int lwork = 0; + superlu::trans_t trans = superlu::NOTRANS; + + superlu::GlobalLU_t Glu; /* Not needed on return. */ + arrayops::fill_zeros(reinterpret_cast<char*>(&Glu), sizeof(superlu::GlobalLU_t)); + + superlu_supermatrix_wrangler x; + superlu_supermatrix_wrangler xC; + + const bool status_x = sp_auxlib::copy_to_supermatrix_with_shift(x.get_ref(), X, sigma); + + if(status_x == false) + { + arma_stop_runtime_error("run_aupd_shiftinvert(): could not construct SuperLU matrix"); + info = blas_int(-1); + return; + } + + // // for debugging only + // if(true) + // { + // cout << "*** testing output of copy_to_supermatrix_with_shift()" << endl; + // cout << "*** sigma: " << sigma << endl; + // + // SpMat<T> Y(X); + // Y.diag() -= sigma; + // + // SpMat<T> Z; + // + // sp_auxlib::copy_to_spmat(Z, x.get_ref()); + // + // cout << "*** size(Y): " << arma::size(Y) << endl; + // cout << "*** size(Z): " << arma::size(Z) << endl; + // cout << "*** accu(abs(Y)): " << accu(abs(Y)) << endl; + // cout << "*** accu(abs(Z)): " << accu(abs(Z)) << endl; + // + // if(arma::size(Y) == arma::size(Z)) + // { + // cout << "*** error: " << accu(abs(Y-Z)) << endl; + // } + // } + + superlu_supermatrix_wrangler l; + superlu_supermatrix_wrangler u; + + superlu_array_wrangler<int> perm_c(X.n_cols+1); // paranoia: increase array length by 1 + superlu_array_wrangler<int> perm_r(X.n_rows+1); + superlu_array_wrangler<int> etree(X.n_cols+1); + + superlu_stat_wrangler stat; + + int panel_size = superlu::sp_ispec_environ(1); + int relax = superlu::sp_ispec_environ(2); + int slu_info = 0; // Return code. + + arma_extra_debug_print("superlu::gstrf()"); + superlu::get_permutation_c(options.ColPerm, x.get_ptr(), perm_c.get_ptr()); + superlu::sp_preorder_mat(&options, x.get_ptr(), perm_c.get_ptr(), etree.get_ptr(), xC.get_ptr()); + superlu::gstrf<T>(&options, xC.get_ptr(), relax, panel_size, etree.get_ptr(), NULL, lwork, perm_c.get_ptr(), perm_r.get_ptr(), l.get_ptr(), u.get_ptr(), &Glu, stat.get_ptr(), &slu_info); + + if(slu_info != 0) + { + arma_debug_warn_level(2, "matrix is singular to working precision"); + info = blas_int(-1); + return; + } + + // NOTE: potential problem with inconsistent/mismatched use of eT and T types + eT x_norm_val = sp_auxlib::norm1<T>(x.get_ptr()); + eT x_rcond = sp_auxlib::lu_rcond<T>(l.get_ptr(), u.get_ptr(), x_norm_val); + + if( (x_rcond < std::numeric_limits<eT>::epsilon()) || arma_isnan(x_rcond) ) + { + arma_debug_warn_level(2, "matrix is singular to working precision (rcond: ", x_rcond, ")"); + info = blas_int(-1); + return; + } + + // All the parameters have been set or created. Time to loop a lot. + while(ido != 99) + { + // Call saupd() or naupd() with the current parameters. + if(sym) + { + arma_extra_debug_print("arpack::saupd()"); + arpack::saupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, &info); + } + else + { + arma_extra_debug_print("arpack::naupd()"); + arpack::naupd(&ido, &bmat, &n, which, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); + } + + // What do we do now? + switch (ido) + { + case -1: + // fallthrough + case 1: + { + // We need to calculate the matrix-vector multiplication y = OP * x + // where x is of length n and starts at workd(ipntr(0)), and y is of + // length n and starts at workd(ipntr(1)). + + // Set the output to point at the right memory. We have to subtract + // one from FORTRAN pointers... + Col<T> out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); + // Set the input to point at the right memory. + Col<T> in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); + + // Consider getting the LU factorization from ZGSTRF, and then + // solve the system L*U*out = in (possibly with permutation matrix?) + // Instead of "spsolve(out,X,in)" we call gstrf above and gstrs below + + out = in; + superlu_supermatrix_wrangler out_slu; + + const bool status_out_slu = sp_auxlib::wrap_to_supermatrix(out_slu.get_ref(), out); + + if(status_out_slu == false) { arma_stop_runtime_error("run_aupd_shiftinvert(): could not construct SuperLU matrix"); return; } + + arma_extra_debug_print("superlu::gstrs()"); + superlu::gstrs<T>(trans, l.get_ptr(), u.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), out_slu.get_ptr(), stat.get_ptr(), &info); + + // No need to modify memory further since it was all done in-place. + + break; + } + case 99: + // Nothing to do here, things have converged. + break; + default: + { + return; // Parent frame can look at the value of info. + } + } + } + + // The process has ended; check the return code. + if( (info != 0) && (info != 1) ) + { + // Print warnings if there was a failure. + + if(sym) + { + arma_debug_warn_level(2, "eigs_sym(): ARPACK error ", info, " in saupd()"); + } + else + { + arma_debug_warn_level(2, "eigs_gen(): ARPACK error ", info, " in naupd()"); + } + + return; // Parent frame can look at the value of info. + } + } + #else + { + arma_ignore(n_eigvals); + arma_ignore(sigma); + arma_ignore(X); + arma_ignore(sym); + arma_ignore(n); + arma_ignore(tol); + arma_ignore(maxiter); + arma_ignore(resid); + arma_ignore(ncv); + arma_ignore(v); + arma_ignore(ldv); + arma_ignore(iparam); + arma_ignore(ipntr); + arma_ignore(workd); + arma_ignore(workl); + arma_ignore(lworkl); + arma_ignore(rwork); + arma_ignore(info); + } + #endif + } + + + +template<typename eT> +inline +bool +sp_auxlib::rudimentary_sym_check(const SpMat<eT>& X) + { + arma_extra_debug_sigprint(); + + if(X.n_rows != X.n_cols) { return false; } + + const eT tol = eT(10000) * std::numeric_limits<eT>::epsilon(); // allow some leeway + + typename SpMat<eT>::const_iterator it = X.begin(); + typename SpMat<eT>::const_iterator it_end = X.end(); + + const uword n_check_limit = (std::max)( uword(2), uword(X.n_nonzero/100) ); + + uword n_check = 1; + + while( (it != it_end) && (n_check <= n_check_limit) ) + { + const uword it_row = it.row(); + const uword it_col = it.col(); + + if(it_row != it_col) + { + const eT A = (*it); + const eT B = X.at( it_col, it_row ); // deliberately swapped + + const eT C = (std::max)(std::abs(A), std::abs(B)); + + const eT delta = std::abs(A - B); + + if(( (delta <= tol) || (delta <= (C * tol)) ) == false) { return false; } + + ++n_check; + } + + ++it; + } + + return true; + } + + + +template<typename T> +inline +bool +sp_auxlib::rudimentary_sym_check(const SpMat< std::complex<T> >& X) + { + arma_extra_debug_sigprint(); + + // NOTE: the function name is a misnomer, as it checks for hermitian complex matrices; + // NOTE: for simplicity of use, the function name is the same as for real matrices + + typedef typename std::complex<T> eT; + + if(X.n_rows != X.n_cols) { return false; } + + const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway + + typename SpMat<eT>::const_iterator it = X.begin(); + typename SpMat<eT>::const_iterator it_end = X.end(); + + const uword n_check_limit = (std::max)( uword(2), uword(X.n_nonzero/100) ); + + uword n_check = 1; + + while( (it != it_end) && (n_check <= n_check_limit) ) + { + const uword it_row = it.row(); + const uword it_col = it.col(); + + if(it_row != it_col) + { + const eT A = (*it); + const eT B = X.at( it_col, it_row ); // deliberately swapped + + const T C_real = (std::max)(std::abs(A.real()), std::abs(B.real())); + const T C_imag = (std::max)(std::abs(A.imag()), std::abs(B.imag())); + + const T delta_real = std::abs(A.real() - B.real()); + const T delta_imag = std::abs(A.imag() + B.imag()); // take into account the conjugate + + const bool okay_real = ( (delta_real <= tol) || (delta_real <= (C_real * tol)) ); + const bool okay_imag = ( (delta_imag <= tol) || (delta_imag <= (C_imag * tol)) ); + + if( (okay_real == false) || (okay_imag == false) ) { return false; } + + ++n_check; + } + else + { + const eT A = (*it); + + if(std::abs(A.imag()) > tol) { return false; } + } + + ++it; + } + + return true; + } + + + +// + + + +template<typename eT> +inline +eigs_randu_filler<eT>::eigs_randu_filler() + { + arma_extra_debug_sigprint(); + + typedef typename std::mt19937_64::result_type local_seed_type; + + local_engine.seed(local_seed_type(123)); + + typedef typename std::uniform_real_distribution<eT>::param_type local_param_type; + + local_u_distr.param(local_param_type(-1.0, +1.0)); + } + + +template<typename eT> +inline +void +eigs_randu_filler<eT>::fill(podarray<eT>& X, const uword N) + { + arma_extra_debug_sigprint(); + + X.set_size(N); + + eT* X_mem = X.memptr(); + + for(uword i=0; i<N; ++i) { X_mem[i] = eT( local_u_distr(local_engine) ); } + } + + +template<typename T> +inline +eigs_randu_filler< std::complex<T> >::eigs_randu_filler() + { + arma_extra_debug_sigprint(); + + typedef typename std::mt19937_64::result_type local_seed_type; + + local_engine.seed(local_seed_type(123)); + + typedef typename std::uniform_real_distribution<T>::param_type local_param_type; + + local_u_distr.param(local_param_type(-1.0, +1.0)); + } + + +template<typename T> +inline +void +eigs_randu_filler< std::complex<T> >::fill(podarray< std::complex<T> >& X, const uword N) + { + arma_extra_debug_sigprint(); + + typedef typename std::complex<T> eT; + + X.set_size(N); + + eT* X_mem = X.memptr(); + + for(uword i=0; i<N; ++i) + { + eT& X_mem_i = X_mem[i]; + + X_mem_i.real( T(local_u_distr(local_engine)) ); + X_mem_i.imag( T(local_u_distr(local_engine)) ); + } + } + + + +// + + + +#if defined(ARMA_USE_SUPERLU) + +inline +superlu_supermatrix_wrangler::~superlu_supermatrix_wrangler() + { + arma_extra_debug_sigprint_this(this); + + if(used == false) { return; } + + char* m_char = reinterpret_cast<char*>(&m); + bool all_zero = true; + + for(size_t i=0; i < sizeof(superlu::SuperMatrix); ++i) + { + if(m_char[i] != char(0)) { all_zero = false; break; } + } + + if(all_zero == false) { sp_auxlib::destroy_supermatrix(m); } + } + +inline +superlu_supermatrix_wrangler::superlu_supermatrix_wrangler() + { + arma_extra_debug_sigprint_this(this); + + arrayops::fill_zeros(reinterpret_cast<char*>(&m), sizeof(superlu::SuperMatrix)); + } + +inline +superlu::SuperMatrix& +superlu_supermatrix_wrangler::get_ref() + { + used = true; + + return m; + } + +inline +superlu::SuperMatrix* +superlu_supermatrix_wrangler::get_ptr() + { + used = true; + + return &m; + } + + +// + + +inline +superlu_stat_wrangler::~superlu_stat_wrangler() + { + arma_extra_debug_sigprint_this(this); + + superlu::free_stat(&stat); + } + +inline +superlu_stat_wrangler::superlu_stat_wrangler() + { + arma_extra_debug_sigprint_this(this); + + arrayops::fill_zeros(reinterpret_cast<char*>(&stat), sizeof(superlu::SuperLUStat_t)); + + superlu::init_stat(&stat); + } + +inline +superlu::SuperLUStat_t* +superlu_stat_wrangler::get_ptr() + { + return &stat; + } + + +// + + +template<typename eT> +inline +superlu_array_wrangler<eT>::~superlu_array_wrangler() + { + arma_extra_debug_sigprint_this(this); + + (*this).reset(); + } + +template<typename eT> +inline +superlu_array_wrangler<eT>::superlu_array_wrangler() + : mem(nullptr) + { + arma_extra_debug_sigprint_this(this); + } + +template<typename eT> +inline +superlu_array_wrangler<eT>::superlu_array_wrangler(const uword n_elem) + : mem(nullptr) + { + arma_extra_debug_sigprint_this(this); + + (*this).set_size(n_elem); + } + +template<typename eT> +inline +void +superlu_array_wrangler<eT>::set_size(const uword n_elem) + { + arma_extra_debug_sigprint(); + + if(mem != nullptr) { (*this).reset(); } + + mem = (eT*)(superlu::malloc(n_elem * sizeof(eT))); + + arma_check_bad_alloc( (mem == nullptr), "superlu::malloc(): out of memory" ); + + arrayops::fill_zeros(mem, n_elem); + } + +template<typename eT> +inline +void +superlu_array_wrangler<eT>::reset() + { + arma_extra_debug_sigprint(); + + if(mem != nullptr) + { + superlu::free(mem); + mem = nullptr; + } + } + +template<typename eT> +inline +eT* +superlu_array_wrangler<eT>::get_ptr() + { + return mem; + } + + +// + + +template<typename eT> +inline +superlu_worker<eT>::~superlu_worker() + { + arma_extra_debug_sigprint_this(this); + + if(l != nullptr) { delete l; l = nullptr; } + if(u != nullptr) { delete u; u = nullptr; } + } + + +template<typename eT> +inline +superlu_worker<eT>::superlu_worker() + { + arma_extra_debug_sigprint_this(this); + } + + +template<typename eT> +inline +bool +superlu_worker<eT>::factorise(typename get_pod_type<eT>::result& out_rcond, const SpMat<eT>& A, const superlu_opts& user_opts) + { + arma_extra_debug_sigprint(); + + typedef typename get_pod_type<eT>::result T; + + factorisation_valid = false; + + if(l != nullptr) { delete l; l = nullptr; } + if(u != nullptr) { delete u; u = nullptr; } + + l = new(std::nothrow) superlu_supermatrix_wrangler; + u = new(std::nothrow) superlu_supermatrix_wrangler; + + if( (l == nullptr) || (u == nullptr) ) + { + arma_debug_warn_level(3, "superlu_worker()::factorise(): could not construct SuperLU matrix"); + return false; + } + + superlu_supermatrix_wrangler& l_ref = (*l); + superlu_supermatrix_wrangler& u_ref = (*u); + + superlu::superlu_options_t options; + sp_auxlib::set_superlu_opts(options, user_opts); + + superlu_supermatrix_wrangler AA; + superlu_supermatrix_wrangler AAc; + + const bool status_AA = sp_auxlib::copy_to_supermatrix(AA.get_ref(), A); + + if(status_AA == false) + { + arma_debug_warn_level(3, "superlu_worker()::factorise(): could not construct SuperLU matrix"); + return false; + } + + (*this).perm_c.set_size(A.n_cols+1); // paranoia: increase array length by 1 + (*this).perm_r.set_size(A.n_rows+1); + + superlu_array_wrangler<int> etree(A.n_cols+1); + + superlu::GlobalLU_t Glu; + arrayops::fill_zeros(reinterpret_cast<char*>(&Glu), sizeof(superlu::GlobalLU_t)); + + int panel_size = superlu::sp_ispec_environ(1); + int relax = superlu::sp_ispec_environ(2); + int lwork = 0; + int info = 0; + + arma_extra_debug_print("superlu::superlu::get_permutation_c()"); + superlu::get_permutation_c(options.ColPerm, AA.get_ptr(), perm_c.get_ptr()); + + arma_extra_debug_print("superlu::superlu::sp_preorder_mat()"); + superlu::sp_preorder_mat(&options, AA.get_ptr(), perm_c.get_ptr(), etree.get_ptr(), AAc.get_ptr()); + + arma_extra_debug_print("superlu::gstrf()"); + superlu::gstrf<eT>(&options, AAc.get_ptr(), relax, panel_size, etree.get_ptr(), NULL, lwork, perm_c.get_ptr(), perm_r.get_ptr(), l_ref.get_ptr(), u_ref.get_ptr(), &Glu, stat.get_ptr(), &info); + + if(info != 0) + { + arma_debug_warn_level(3, "superlu_worker()::factorise(): LU factorisation failed"); + return false; + } + + const T AA_norm = sp_auxlib::norm1<T>(AA.get_ptr()); + const T AA_rcond = sp_auxlib::lu_rcond<eT>(l_ref.get_ptr(), u_ref.get_ptr(), AA_norm); + + out_rcond = AA_rcond; + + if(arma_isnan(AA_rcond)) { return false; } + // if(AA_rcond == T(0)) { return false; } + + factorisation_valid = true; + + return true; + } + + +template<typename eT> +inline +bool +superlu_worker<eT>::solve(Mat<eT>& X, const Mat<eT>& B) + { + arma_extra_debug_sigprint(); + + if(factorisation_valid == false) { return false; } + if( (l == nullptr) || (u == nullptr) ) { return false; } + + superlu_supermatrix_wrangler& l_ref = (*l); + superlu_supermatrix_wrangler& u_ref = (*u); + + X = B; + + superlu_supermatrix_wrangler XX; + + const bool status_XX = sp_auxlib::wrap_to_supermatrix(XX.get_ref(), X); + + if(status_XX == false) + { + arma_debug_warn_level(3, "superlu_worker()::solve(): could not construct SuperLU matrix"); + return false; + } + + superlu::trans_t trans = superlu::NOTRANS; + int info = 0; + + arma_extra_debug_print("superlu::gstrs()"); + superlu::gstrs<eT>(trans, l_ref.get_ptr(), u_ref.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), XX.get_ptr(), stat.get_ptr(), &info); + + return (info == 0); + } + + +#endif + + +//! @} |