CPastixWrapper.hpp

Go to the documentation of this file.

/*!
 * \file CPastixWrapper.hpp
 * \brief An interface to the INRIA solver PaStiX
 *        (http://pastix.gforge.inria.fr/files/README-txt.html)
 * \author P. Gomes
 * \version 8.5.0 "Harrier"
 *
 * SU2 Project Website: https://su2code.github.io
 *
 * The SU2 Project is maintained by the SU2 Foundation
 * (http://su2foundation.org)
 *
 * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
 *
 * SU2 is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * SU2 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
 */
 
#pragma once
 
#include "../code_config.hpp"
 
#ifdef HAVE_PASTIX
 
#ifdef CODI_FORWARD_TYPE
#error Cannot use PaStiX with forward mode AD
#endif
 
#include <pastix.h>
#include <spm.h>
#include <vector>
 
using namespace std;
 
class CConfig;
class CGeometry;
 
/*!
 * \class CPastixWrapper
 * \ingroup SpLinSys
 * \brief Wrapper class that converts between SU2 sparse format and PaStiX
 *        format and simplifies calls to the external solver.
 */
template <class ScalarType>
class CPastixWrapper {
 private:
  pastix_data_t* state{};        /*!< \brief Internal state of the solver. */
  spmatrix_t spm;                /*!< \brief Matrix format used by the solver. */
  pastix_int_t nCols;            /*!< \brief Local number of columns. */
  vector<pastix_int_t> colptr;   /*!< \brief Equiv. to our "row_ptr". */
  vector<pastix_int_t> rowidx;   /*!< \brief Equiv. to our "col_ind". */
  vector<su2mixedfloat> values;  /*!< \brief Equiv. to our "matrix". */
  vector<pastix_int_t> loc2glb;  /*!< \brief Global index of the columns held by this rank. */
  vector<pastix_int_t> perm;     /*!< \brief Ordering computed by PaStiX. */
  vector<su2mixedfloat> workvec; /*!< \brief RHS vector which then becomes the solution. */
 
  vector<unsigned long> csr_row_ptr; /*!< \brief Owned CSR row pointers (built from LDU). */
  vector<unsigned long> csr_col_ind; /*!< \brief Owned CSR column indices (built from LDU). */
 
  pastix_int_t iparm[IPARM_SIZE]; /*!< \brief Integer parameters for PaStiX. */
  double dparm[DPARM_SIZE];       /*!< \brief Floating point parameters for PaStiX. */
 
  struct {
    unsigned long nVar = 0;
    unsigned long nPoint = 0;
    unsigned long nPointDomain = 0;
    unsigned long blkSz = 0; /*!< \brief Block size (nVar * nVar) for value assembly. */
 
    const su2uint* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */
    const su2uint* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */
    const ScalarType* d = nullptr;      /*!< \brief Diagonal blocks (matrix-owned). */
    const ScalarType* l = nullptr;      /*!< \brief Lower blocks (matrix-owned). */
    const ScalarType* u = nullptr;      /*!< \brief Upper blocks (matrix-owned). */
 
    unsigned long size_rhs() const { return nPointDomain * nVar; }
  } matrix; /*!< \brief Dimensions and LDU pointers captured from the owning CSysMatrix. */
 
  bool issetup{};        /*!< \brief Signals that the structure has been provided. */
  bool isinitialized{};  /*!< \brief Signals that the sparsity pattern has been set. */
  bool isfactorized{};   /*!< \brief Signals that a factorization has been computed. */
  bool transpose{};      /*!< \brief Solve A^T x = b instead of A x = b. */
  unsigned long iter{};  /*!< \brief Number of times a factorization has been requested. */
  unsigned short verb{}; /*!< \brief Verbosity level. */
  const int mpi_size = SU2_MPI::GetSize();
  const int mpi_rank = SU2_MPI::GetRank();
 
  vector<unsigned long> sort_rows;          /*!< \brief List of rows with halo points. */
  vector<vector<unsigned long>> sort_order; /*!< \brief How each of those rows needs to be sorted. */
 
  /*!
   * \brief Run the "clean" task, releases all memory, leaves object in unusable state.
   */
  void Clean() {
    if (isinitialized) {
      iparm[IPARM_VERBOSE] = (verb > 0) ? PastixVerboseNo : PastixVerboseNot;
      pastixFinalize(&state);
      spm.colptr = nullptr;
      spm.rowptr = nullptr;
      spm.values = nullptr;
      if (mpi_size > 1) spm.loc2glob = nullptr;
      spmExit(&spm);
      isfactorized = false;
      isinitialized = false;
    }
  }
 
  /*!
   * \brief Initialize the matrix format that PaStiX requires.
   */
  void Initialize(CGeometry* geometry, const CConfig* config);
 
  /*!
   * \brief Assemble CSR values from the stored LDU pointers directly into the values buffer.
   */
  void AssembleValues();
 
 public:
  CPastixWrapper() = default;
 
  /*--- Move or copy is not allowed. ---*/
  CPastixWrapper(CPastixWrapper&&) = delete;
  CPastixWrapper(const CPastixWrapper&) = delete;
  CPastixWrapper& operator=(CPastixWrapper&&) = delete;
  CPastixWrapper& operator=(const CPastixWrapper&) = delete;
 
  /*!
   * \brief Class destructor.
   */
  ~CPastixWrapper() { Clean(); }
 
  /*!
   * \brief Returns true once SetLDU has been called.
   */
  bool IsSetup() const { return issetup; }
 
  /*!
   * \brief Set LDU structure and value pointers; builds and owns assembled CSR (called once).
   * \param[in] nVar - DOF per point (square blocks: nVar x nVar).
   * \param[in] nPoint - Total number of points including halos.
   * \param[in] nPointDomain - Number of internal points (domain rows).
   * \param[in] row_ptr_l/u - LDU lower/upper row pointers (geometry-owned, must outlive wrapper).
   * \param[in] col_ind_l/u - LDU lower/upper column indices (geometry-owned).
   * \param[in] d/l/u - LDU value blocks (matrix-owned, must outlive wrapper).
   */
  void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const su2uint* row_ptr_l,
              const su2uint* col_ind_l, const su2uint* row_ptr_u, const su2uint* col_ind_u, const ScalarType* d,
              const ScalarType* l, const ScalarType* u) {
    if (issetup) return;
    matrix.nVar = nVar;
    matrix.nPoint = nPoint;
    matrix.nPointDomain = nPointDomain;
    matrix.row_ptr_l = row_ptr_l;
    matrix.row_ptr_u = row_ptr_u;
    matrix.d = d;
    matrix.l = l;
    matrix.u = u;
    matrix.blkSz = nVar * nVar;
 
    const unsigned long nnz_domain = row_ptr_l[nPointDomain] + nPointDomain + row_ptr_u[nPointDomain];
    csr_row_ptr.resize(nPointDomain + 1);
    csr_col_ind.reserve(nnz_domain);
    for (auto i = 0ul; i < nPointDomain; ++i) {
      csr_row_ptr[i] = static_cast<unsigned long>(csr_col_ind.size());
      for (auto k = row_ptr_l[i]; k < row_ptr_l[i + 1]; ++k) csr_col_ind.push_back(col_ind_l[k]);
      csr_col_ind.push_back(i);
      for (auto k = row_ptr_u[i]; k < row_ptr_u[i + 1]; ++k) csr_col_ind.push_back(col_ind_u[k]);
    }
    csr_row_ptr[nPointDomain] = static_cast<unsigned long>(csr_col_ind.size());
    issetup = true;
  }
 
  /*!
   * \brief Factorize matrix.
   * \param[in] geometry - Geometrical definition of the problem.
   * \param[in] config   - Definition of the particular problem.
   * \param[in] kind_fact - Type of factorization.
   */
  void Factorize(CGeometry* geometry, const CConfig* config, unsigned short kind_fact);
 
  /*!
   * \brief Request solves with the transposed matrix.
   * \param[in] transposed - Yes or no.
   */
  void SetTransposedSolve(bool transposed = true) { transpose = transposed; }
 
  /*!
   * \brief Runs the "solve" task for any rhs/sol with operator []
   * \param[in]  rhs - Right hand side of the linear system.
   * \param[out] sol - Solution of the system.
   */
  template <class T>
  void Solve(const T& rhs, T& sol) {
    if (!isfactorized) SU2_MPI::Error("The factorization has not been computed yet.", CURRENT_FUNCTION);
 
    /*--- Inverted logic due to CSR to CSC direct copy (PaStiX's A is our A^T). ---*/
    if (iparm[IPARM_FACTORIZATION] == PastixFactLDLT || transpose) {
      iparm[IPARM_TRANSPOSE_SOLVE] = PastixNoTrans;
    } else {
      iparm[IPARM_TRANSPOSE_SOLVE] = PastixTrans;
    }
    iparm[IPARM_VERBOSE] = PastixVerboseNot;
 
    for (auto i = 0ul; i < matrix.size_rhs(); ++i) workvec[i] = rhs[i];
    if (pastix_task_solve(state, matrix.size_rhs(), 1, workvec.data(), matrix.size_rhs()) != PASTIX_SUCCESS) {
      SU2_MPI::Error("Error solving linear system.", CURRENT_FUNCTION);
    }
    for (auto i = 0ul; i < matrix.size_rhs(); ++i) sol[i] = workvec[i];
  }
};
#endif