Skip to content

Commit

Permalink
Merge pull request #1228 from su2code/linsol_fixes
Browse files Browse the repository at this point in the history 
Linear solver changes to support hybrid parallel AD
  • Loading branch information
@pcarruscag
pcarruscag authored Mar 22, 2021
2 parents 19826aa + 6f2230d commit 949f4e5
Show file tree
Hide file tree
Showing 23 changed files with 350 additions and 390 deletions.
40 changes: 0 additions & 40 deletions Common/include/linear_algebra/CMatrixVectorProduct.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,43 +99,3 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
matrix.MatrixVectorProduct(u, v, geometry, config);
}
};


/*!
* \class CSysMatrixVectorProductTransposed
* \brief Specialization of matrix-vector product that uses CSysMatrix class for transposed products
*/
template<class ScalarType>
class CSysMatrixVectorProductTransposed final : public CMatrixVectorProduct<ScalarType> {
private:
const CSysMatrix<ScalarType>& matrix; /*!< \brief pointer to matrix that defines the product. */
CGeometry* geometry; /*!< \brief geometry associated with the matrix. */
const CConfig *config; /*!< \brief config of the problem. */

public:
/*!
* \brief constructor of the class
* \param[in] matrix_ref - matrix reference that will be used to define the products
* \param[in] geometry_ref - geometry associated with the problem
* \param[in] config_ref - config of the problem
*/
inline CSysMatrixVectorProductTransposed(const CSysMatrix<ScalarType> & matrix_ref,
CGeometry *geometry_ref, const CConfig *config_ref) :
matrix(matrix_ref),
geometry(geometry_ref),
config(config_ref) {}

/*!
* \note This class cannot be default constructed as that would leave us with invalid pointers.
*/
CSysMatrixVectorProductTransposed() = delete;

/*!
* \brief operator that defines the CSysMatrix-CSysVector product
* \param[in] u - CSysVector that is being multiplied by the sparse matrix
* \param[out] v - CSysVector that is the result of the product
*/
inline void operator()(const CSysVector<ScalarType> & u, CSysVector<ScalarType> & v) const override {
matrix.MatrixVectorProductTransposed(u, v, geometry, config);
}
};
13 changes: 11 additions & 2 deletions Common/include/linear_algebra/CPastixWrapper.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,9 +167,18 @@ class CPastixWrapper
* \param[in] geometry - Geometrical definition of the problem.
* \param[in] config - Definition of the particular problem.
* \param[in] kind_fact - Type of factorization.
* \param[in] transposed - Flag to use the transposed matrix during application of the preconditioner.
*/
void Factorize(CGeometry *geometry, const CConfig *config, unsigned short kind_fact, bool transposed);
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) {
using namespace PaStiX;
if (iparm[IPARM_SYM] == API_SYM_NO)
iparm[IPARM_TRANSPOSE_SOLVE] = pastix_int_t(!transposed); // negated due to CSR to CSC copy
}

/*!
* \brief Runs the "solve" task for any rhs/sol with operator []
Expand Down
64 changes: 47 additions & 17 deletions Common/include/linear_algebra/CPreconditioner.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@
* License along with SU2. If not, see <http://www.gnu.org/licenses/>.
*/


#pragma once

#include "../CConfig.hpp"
Expand Down Expand Up @@ -59,6 +58,17 @@ class CPreconditioner {
* \brief Generic "preprocessing" hook derived classes may implement to build the preconditioner.
*/
virtual void Build() {}

/*!
* \brief Return true to identify the identity preconditioner, may allow some solvers to be more efficient.
*/
virtual bool IsIdentity() const { return false; }

/*!
* \brief Factory method.
*/
static CPreconditioner* Create(ENUM_LINEAR_SOLVER_PREC kind, CSysMatrix<ScalarType>& jacobian,
CGeometry* geometry, const CConfig* config);
};
template<class ScalarType>
CPreconditioner<ScalarType>::~CPreconditioner() {}
Expand All @@ -74,25 +84,22 @@ class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
CGeometry* geometry; /*!< \brief Pointer to geometry associated with the matrix. */
const CConfig *config; /*!< \brief Pointer to problem configuration. */
bool transp; /*!< \brief If the transpose version of the preconditioner is required. */

public:
/*!
* \brief Constructor of the class.
* \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
* \param[in] geometry_ref - Geometry associated with the problem.
* \param[in] config_ref - Config of the problem.
* \param[in] transposed - If the transpose version of the preconditioner is required.
*/
inline CJacobiPreconditioner(CSysMatrix<ScalarType> & matrix_ref,
CGeometry *geometry_ref, const CConfig *config_ref, bool transposed) :
CGeometry *geometry_ref, const CConfig *config_ref) :
sparse_matrix(matrix_ref)
{
if((geometry_ref == nullptr) || (config_ref == nullptr))
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
geometry = geometry_ref;
config = config_ref;
transp = transposed;
}

/*!
Expand All @@ -113,7 +120,7 @@ class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
* \note Request the associated matrix to build the preconditioner.
*/
inline void Build() override {
sparse_matrix.BuildJacobiPreconditioner(transp);
sparse_matrix.BuildJacobiPreconditioner();
}
};

Expand All @@ -128,25 +135,22 @@ class CILUPreconditioner final : public CPreconditioner<ScalarType> {
CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
CGeometry* geometry; /*!< \brief Pointer to geometry associated with the matrix. */
const CConfig *config; /*!< \brief Pointer to problem configuration. */
bool transp; /*!< \brief If the transpose version of the preconditioner is required. */

public:
/*!
* \brief Constructor of the class.
* \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
* \param[in] geometry_ref - Geometry associated with the problem.
* \param[in] config_ref - Config of the problem.
* \param[in] transposed - If the transpose version of the preconditioner is required.
*/
inline CILUPreconditioner(CSysMatrix<ScalarType> & matrix_ref,
CGeometry *geometry_ref, const CConfig *config_ref, bool transposed) :
CGeometry *geometry_ref, const CConfig *config_ref) :
sparse_matrix(matrix_ref)
{
if((geometry_ref == nullptr) || (config_ref == nullptr))
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
geometry = geometry_ref;
config = config_ref;
transp = transposed;
}

/*!
Expand All @@ -167,7 +171,7 @@ class CILUPreconditioner final : public CPreconditioner<ScalarType> {
* \note Request the associated matrix to build the preconditioner.
*/
inline void Build() override {
sparse_matrix.BuildILUPreconditioner(transp);
sparse_matrix.BuildILUPreconditioner();
}
};

Expand Down Expand Up @@ -263,7 +267,7 @@ class CLineletPreconditioner final : public CPreconditioner<ScalarType> {
* \note Request the associated matrix to build the preconditioner.
*/
inline void Build() override {
sparse_matrix.BuildJacobiPreconditioner(false);
sparse_matrix.BuildJacobiPreconditioner();
}
};

Expand All @@ -279,7 +283,6 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
CGeometry* geometry; /*!< \brief Geometry associated with the problem. */
const CConfig *config; /*!< \brief Configuration of the problem. */
unsigned short kind_fact; /*!< \brief The type of factorization desired. */
bool transp; /*!< \brief If the transpose version of the preconditioner is required. */

public:
/*!
Expand All @@ -288,18 +291,16 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
* \param[in] geometry_ref - Associated geometry.
* \param[in] config_ref - Problem configuration.
* \param[in] kind_factorization - Type of factorization required.
* \param[in] transposed - If the transpose version of the preconditioner is required.
*/
inline CPastixPreconditioner(CSysMatrix<ScalarType> & matrix_ref, CGeometry *geometry_ref,
const CConfig *config_ref, unsigned short kind_factorization, bool transposed) :
const CConfig *config_ref, unsigned short kind_factorization) :
sparse_matrix(matrix_ref)
{
if((geometry_ref == nullptr) || (config_ref == nullptr))
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
geometry = geometry_ref;
config = config_ref;
kind_fact = kind_factorization;
transp = transposed;
}

/*!
Expand All @@ -320,6 +321,35 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
* \note Request the associated matrix to build the preconditioner.
*/
inline void Build() override {
sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact, transp);
sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact);
}
};


template<class ScalarType>
CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOLVER_PREC kind,
CSysMatrix<ScalarType>& jacobian,
CGeometry* geometry,
const CConfig* config) {
CPreconditioner<ScalarType>* prec = nullptr;

switch (kind) {
case JACOBI:
prec = new CJacobiPreconditioner<ScalarType>(jacobian, geometry, config);
break;
case LINELET:
prec = new CLineletPreconditioner<ScalarType>(jacobian, geometry, config);
break;
case LU_SGS:
prec = new CLU_SGSPreconditioner<ScalarType>(jacobian, geometry, config);
break;
case ILU:
prec = new CILUPreconditioner<ScalarType>(jacobian, geometry, config);
break;
case PASTIX_ILU: case PASTIX_LU_P: case PASTIX_LDLT_P:
prec = new CPastixPreconditioner<ScalarType>(jacobian, geometry, config, kind);
break;
}

return prec;
}
64 changes: 17 additions & 47 deletions Common/include/linear_algebra/CSysMatrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,8 +164,6 @@ class CSysMatrix {
gemm_t MatrixVectorProductKernelBetaOne; /*!< \brief MKL JIT based GEMV kernel with BETA=1.0. */
void * MatrixVectorProductJitterAlphaMinusOne; /*!< \brief Jitter handle for MKL JIT based GEMV with ALPHA=-1.0 and BETA=1.0. */
gemm_t MatrixVectorProductKernelAlphaMinusOne; /*!< \brief MKL JIT based GEMV kernel with ALPHA=-1.0 and BETA=1.0. */
void * MatrixVectorProductTranspJitterBetaOne; /*!< \brief Jitter handle for MKL JIT based GEMV (transposed) with BETA=1.0. */
gemm_t MatrixVectorProductTranspKernelBetaOne; /*!< \brief MKL JIT based GEMV (transposed) kernel with BETA=1.0. */
#endif

#ifdef HAVE_PASTIX
Expand All @@ -177,6 +175,9 @@ class CSysMatrix {
*/
struct {
const unsigned long *ptr = nullptr;
unsigned long nEdge = 0;

operator bool() { return nEdge != 0; }

inline unsigned long operator() (unsigned long edge, unsigned long node) const {
return ptr[2*edge+node];
Expand Down Expand Up @@ -225,14 +226,6 @@ class CSysMatrix {
*/
void MatrixVectorProductSub(const ScalarType *matrix, const ScalarType *vector, ScalarType *product) const;

/*!
* \brief Calculates the matrix-vector product: product += matrix^T * vector
* \param[in] matrix
* \param[in] vector
* \param[in,out] product
*/
void MatrixVectorProductTransp(const ScalarType *matrix, const ScalarType *vector, ScalarType *product) const;

/*!
* \brief Calculates the matrix-matrix product
*/
Expand All @@ -258,17 +251,10 @@ class CSysMatrix {
/*!
* \brief Copy matrix src into dst, transpose if required.
*/
FORCEINLINE void MatrixCopy(const ScalarType *src, ScalarType *dst, bool transposed = false) const {
if (!transposed) {
SU2_OMP_SIMD
for(auto iVar = 0ul; iVar < nVar*nEqn; ++iVar)
dst[iVar] = src[iVar];
}
else {
for (auto iVar = 0ul; iVar < nVar; ++iVar)
for (auto jVar = 0ul; jVar < nVar; ++jVar)
dst[iVar*nVar+jVar] = src[jVar*nVar+iVar];
}
FORCEINLINE void MatrixCopy(const ScalarType *src, ScalarType *dst) const {
SU2_OMP_SIMD
for(auto iVar = 0ul; iVar < nVar*nEqn; ++iVar)
dst[iVar] = src[iVar];
}

/*!
Expand All @@ -289,17 +275,16 @@ class CSysMatrix {
* \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i,i) subblock and rhs.
* \param[in] block_i - Index of the (i,i) diagonal block.
* \param[in] rhs - Right-hand-side of the linear system.
* \param[in] transposed - If true the transposed of the block is used (default = false).
* \return Solution of the linear system (overwritten on rhs).
*/
inline void Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed = false) const;
inline void Gauss_Elimination(unsigned long block_i, ScalarType* rhs) const;

/*!
* \brief Inverse diagonal block.
* \param[in] block_i - Indexes of the block in the matrix-by-blocks structure.
* \param[out] invBlock - Inverse block.
*/
inline void InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transposed = false) const;
inline void InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock) const;

/*!
* \brief Inverse diagonal block.
Expand All @@ -323,14 +308,6 @@ class CSysMatrix {
*/
inline void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block);

/*!
* \brief Set the transposed value of a block in the sparse matrix.
* \param[in] block_i - Indexes of the block in the matrix-by-blocks structure.
* \param[in] block_j - Indexes of the block in the matrix-by-blocks structure.
* \param[in] **val_block - Block to set to A(i, j).
*/
inline void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block);

/*!
* \brief Performs the product of i-th row of the upper part of a sparse matrix by a vector.
* \param[in] vec - Vector to be multiplied by the upper part of the sparse matrix A.
Expand Down Expand Up @@ -805,6 +782,11 @@ class CSysMatrix {
*/
void SetDiagonalAsColumnSum();

/*!
* \brief Transposes the matrix, any preconditioner that was computed may be invalid.
*/
void TransposeInPlace();

/*!
* \brief Add a scaled sparse matrix to "this" (axpy-type operation, A = A+alpha*B).
* \note Matrices must have the same sparse pattern.
Expand All @@ -823,20 +805,10 @@ class CSysMatrix {
void MatrixVectorProduct(const CSysVector<ScalarType> & vec, CSysVector<ScalarType> & prod,
CGeometry *geometry, const CConfig *config) const;

/*!
* \brief Performs the product of a sparse matrix by a CSysVector.
* \param[in] vec - CSysVector to be multiplied by the sparse matrix A.
* \param[in] geometry - Geometrical definition of the problem.
* \param[in] config - Definition of the particular problem.
* \param[out] prod - Result of the product.
*/
void MatrixVectorProductTransposed(const CSysVector<ScalarType> & vec, CSysVector<ScalarType> & prod,
CGeometry *geometry, const CConfig *config) const;

/*!
* \brief Build the Jacobi preconditioner.
*/
void BuildJacobiPreconditioner(bool transpose = false);
void BuildJacobiPreconditioner();

/*!
* \brief Multiply CSysVector by the preconditioner
Expand All @@ -850,9 +822,8 @@ class CSysMatrix {

/*!
* \brief Build the ILU preconditioner.
* \param[in] transposed - Flag to use the transposed matrix to construct the preconditioner.
*/
void BuildILUPreconditioner(bool transposed = false);
void BuildILUPreconditioner();

/*!
* \brief Multiply CSysVector by the preconditioner
Expand Down Expand Up @@ -902,9 +873,8 @@ class CSysMatrix {
* \param[in] geometry - Geometrical definition of the problem.
* \param[in] config - Definition of the particular problem.
* \param[in] kind_fact - Type of factorization.
* \param[in] transposed - Flag to use the transposed matrix during application of the preconditioner.
*/
void BuildPastixPreconditioner(CGeometry *geometry, const CConfig *config, unsigned short kind_fact, bool transposed = false);
void BuildPastixPreconditioner(CGeometry *geometry, const CConfig *config, unsigned short kind_fact);

/*!
* \brief Apply the PaStiX factorization to CSysVec.
Expand Down
Loading

0 comments on commit 949f4e5

Please sign in to comment.