SparseSolverBase.hpp

Go to the documentation of this file.

/*
 * STRUMPACK -- STRUctured Matrices PACKage, Copyright (c) 2014, The
 * Regents of the University of California, through Lawrence Berkeley
 * National Laboratory (subject to receipt of any required approvals
 * from the U.S. Dept. of Energy).  All rights reserved.
 *
 * If you have questions about your rights to use or distribute this
 * software, please contact Berkeley Lab's Technology Transfer
 * Department at TTD@lbl.gov.
 *
 * NOTICE. This software is owned by the U.S. Department of Energy. As
 * such, the U.S. Government has been granted for itself and others
 * acting on its behalf a paid-up, nonexclusive, irrevocable,
 * worldwide license in the Software to reproduce, prepare derivative
 * works, and perform publicly and display publicly.  Beginning five
 * (5) years after the date permission to assert copyright is obtained
 * from the U.S. Department of Energy, and subject to any subsequent
 * five (5) year renewals, the U.S. Government is granted for itself
 * and others acting on its behalf a paid-up, nonexclusive,
 * irrevocable, worldwide license in the Software to reproduce,
 * prepare derivative works, distribute copies to the public, perform
 * publicly and display publicly, and to permit others to do so.
 *
 * Developers: Pieter Ghysels, Francois-Henry Rouet, Xiaoye S. Li.
 *             (Lawrence Berkeley National Lab, Computational Research
 *             Division).
 */
#ifndef STRUMPACK_SPARSE_SOLVER_BASE_HPP
#define STRUMPACK_SPARSE_SOLVER_BASE_HPP
 
#include <new>
#include <memory>
#include <vector>
#include <string>
 
#include "StrumpackConfig.hpp"
#include "StrumpackOptions.hpp"
#include "sparse/CSRMatrix.hpp"
#include "dense/DenseMatrix.hpp"
 
namespace strumpack {
 
  // forward declarations
  template<typename scalar_t,typename integer_t> class MatrixReordering;
  template<typename scalar_t,typename integer_t> class EliminationTree;
  class TaskTimer;
 
  template<typename scalar_t,typename integer_t=int>
  class SparseSolverBase {
 
    using SpMat_t = CompressedSparseMatrix<scalar_t,integer_t>;
    using Tree_t = EliminationTree<scalar_t,integer_t>;
    using Reord_t = MatrixReordering<scalar_t,integer_t>;
    using DenseM_t = DenseMatrix<scalar_t>;
    using DenseMW_t = DenseMatrixWrapper<scalar_t>;
 
  public:
 
    SparseSolverBase(int argc, char* argv[],
                     bool verbose=true, bool root=true);
 
    SparseSolverBase(bool verbose=true, bool root=true);
 
    virtual ~SparseSolverBase();
 
    ReturnCode reorder(int nx=1, int ny=1, int nz=1,
                       int components=1, int width=1);
 
    ReturnCode reorder(const int* p, int base=0);
 
    ReturnCode factor();
 
    void move_to_gpu();
    void remove_from_gpu();
 
    ReturnCode solve(const scalar_t* b, scalar_t* x,
                     bool use_initial_guess=false);
 
    ReturnCode solve(const DenseM_t& b, DenseM_t& x,
                     bool use_initial_guess=false);
 
    ReturnCode solve(int nrhs, const scalar_t* b, int ldb,
                     scalar_t* x, int ldx,
                     bool use_initial_guess=false);
 
    SPOptions<scalar_t>& options();
 
    const SPOptions<scalar_t>& options() const;
 
    void set_from_options();
 
    void set_from_options(int argc, char* argv[]);
 
    int maximum_rank() const;
 
    std::size_t factor_nonzeros() const;
 
    std::size_t factor_memory() const;
 
    int Krylov_iterations() const;
 
 
    ReturnCode inertia(integer_t& neg, integer_t& zero, integer_t& pos);
 
    ReturnCode subnormals(std::size_t& ns, std::size_t& nz);
 
    ReturnCode pivot_growth(scalar_t& pg);
 
    void draw(const std::string& name) const;
 
    void delete_factors();
 
  protected:
    virtual void setup_tree() = 0;
    virtual void setup_reordering() = 0;
    virtual
    int compute_reordering(const int* p, int base,
                           int nx, int ny, int nz,
                           int components, int width) = 0;
    virtual void separator_reordering() = 0;
 
    virtual SpMat_t* matrix() = 0;
    virtual Reord_t* reordering() = 0;
    virtual Tree_t* tree() = 0;
    virtual const SpMat_t* matrix() const = 0;
    virtual const Reord_t* reordering() const = 0;
    virtual const Tree_t* tree() const = 0;
 
    virtual void perf_counters_start();
    virtual void perf_counters_stop(const std::string& s);
 
    virtual void synchronize() {}
    virtual void communicate_ordering() {}
    virtual double max_peak_memory() const
    { return double(params::peak_memory); }
    virtual double min_peak_memory() const
    { return double(params::peak_memory); }
 
    void papi_initialize();
    long long dense_factor_nonzeros() const;
    void print_solve_stats(TaskTimer& t) const;
 
    virtual void reduce_flop_counters() const {}
    void print_flop_breakdown_HSS() const;
    void print_flop_breakdown_HODLR() const;
    void flop_breakdown_reset() const;
 
    void print_wrong_sparsity_error();
 
    // TODO do these all need to be virtual, can some be private?
    virtual
    ReturnCode solve_internal(const scalar_t* b, scalar_t* x,
                              bool use_initial_guess=false) = 0;
    virtual
    ReturnCode solve_internal(const DenseM_t& b, DenseM_t& x,
                              bool use_initial_guess=false) = 0;
    virtual
    ReturnCode solve_internal(int nrhs, const scalar_t* b, int ldb,
                              scalar_t* x, int ldx,
                              bool use_initial_guess=false);
 
    SPOptions<scalar_t> opts_;
    bool is_root_;
 
    MatchingData<scalar_t,integer_t> matching_;
    Equilibration<scalar_t> equil_;
 
    std::new_handler old_handler_;
    std::ostream* rank_out_ = nullptr;
    bool factored_ = false;
    bool reordered_ = false;
    int Krylov_its_ = 0;
 
#if defined(STRUMPACK_USE_PAPI)
    float rtime_ = 0., ptime_ = 0.;
    long_long _flpops = 0;
#endif
#if defined(STRUMPACK_COUNT_FLOPS)
    long long int f0_ = 0, ftot_ = 0, fmin_ = 0, fmax_ = 0;
    long long int b0_ = 0, btot_ = 0, bmin_ = 0, bmax_ = 0;
    long long int m0_ = 0, mtot_ = 0, mmin_ = 0, mmax_ = 0;
    long long int ptot_ = 0, pmin_ = 0, pmax_ = 0;
    long long int dm0_ = 0, dmtot_ = 0, dmmin_ = 0, dmmax_ = 0;
    long long int dptot_ = 0, dpmin_ = 0, dpmax_ = 0;
#endif
 
  private:
    ReturnCode reorder_internal(const int* p, int base,
                                int nx, int ny, int nz,
                                int components, int width);
 
    virtual void delete_factors_internal() = 0;
  };
 
  template<typename scalar_t,typename integer_t>
  using StrumpackSparseSolverBase =
    SparseSolverBase<scalar_t,integer_t>;
 
} //end namespace strumpack
 
#endif // STRUMPACK_SPARSE_SOLVER_BASE_HPP