StrumpackSparseSolver.hpp

/*
 * 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_HPP
#define STRUMPACK_SPARSE_SOLVER_HPP
 
#include <new>
#include <memory>
#include <vector>
#include <string>
 
#include "SparseSolverBase.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 SparseSolver :
    public SparseSolverBase<scalar_t,integer_t> {
 
    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:
 
    SparseSolver(int argc, char* argv[], bool verbose=true, bool root=true);
 
    SparseSolver(bool verbose=true, bool root=true);
 
    ~SparseSolver();
 
    void set_matrix(const CSRMatrix<scalar_t,integer_t>& A);
 
    void set_csr_matrix(integer_t N,
                        const integer_t* row_ptr, const integer_t* col_ind,
                        const scalar_t* values, bool symmetric_pattern=false);
 
    void update_matrix_values(integer_t N,
                              const integer_t* row_ptr,
                              const integer_t* col_ind,
                              const scalar_t* values,
                              bool symmetric_pattern=false);
 
    void update_matrix_values(const CSRMatrix<scalar_t,integer_t>& A);
 
  private:
    void setup_tree() override;
    void setup_reordering() override;
    int compute_reordering(const int* p, int base,
                           int nx, int ny, int nz,
                           int components, int width) override;
    void separator_reordering() override;
 
    SpMat_t* matrix() override { return mat_.get(); }
    std::unique_ptr<SpMat_t> matrix_nonzero_diag() override {
      return mat_->add_missing_diagonal(opts_.pivot_threshold());
    }
    Reord_t* reordering() override { return nd_.get(); }
    Tree_t* tree() override { return tree_.get(); }
    const SpMat_t* matrix() const override { return mat_.get(); }
    const Reord_t* reordering() const override { return nd_.get(); }
    const Tree_t* tree() const override { return tree_.get(); }
 
    void permute_matrix_values();
 
    ReturnCode solve_internal(const scalar_t* b, scalar_t* x,
                              bool use_initial_guess=false) override;
    ReturnCode solve_internal(const DenseM_t& b, DenseM_t& x,
                              bool use_initial_guess=false) override;
 
    void delete_factors_internal() override;
 
    void transform_x0(DenseM_t& x, DenseM_t& xtmp);
    void transform_b(const DenseM_t& b, DenseM_t& bloc);
    void transform_x(DenseM_t& x, DenseM_t& xtmp);
 
    std::unique_ptr<CSRMatrix<scalar_t,integer_t>> mat_;
    std::unique_ptr<MatrixReordering<scalar_t,integer_t>> nd_;
    std::unique_ptr<EliminationTree<scalar_t,integer_t>> tree_;
 
    using SPBase_t = SparseSolverBase<scalar_t,integer_t>;
    using SPBase_t::opts_;
    using SPBase_t::is_root_;
    using SPBase_t::matching_;
    using SPBase_t::equil_;
    using SPBase_t::factored_;
    using SPBase_t::reordered_;
    using SPBase_t::Krylov_its_;
  };
 
  template<typename scalar_t,typename integer_t>
  using StrumpackSparseSolver = SparseSolver<scalar_t,integer_t>;
 
} //end namespace strumpack
 
#endif // STRUMPACK_SPARSE_SOLVER_HPP