Skip to content

Commit 949f4e5

Browse files
pcarruscagpcarruscag
authored
Merge pull request #1228 from su2code/linsol_fixes
Linear solver changes to support hybrid parallel AD
2 parents 19826aa + 6f2230d commit 949f4e5

File tree

23 files changed

+350
-390
lines changed

23 files changed

+350
-390
lines changed

‎Common/include/linear_algebra/CMatrixVectorProduct.hpp‎

Lines changed: 0 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -99,43 +99,3 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
9999
matrix.MatrixVectorProduct(u, v, geometry, config);
100100
}
101101
};
102-
103-
104-
/*!
105-
* \class CSysMatrixVectorProductTransposed
106-
* \brief Specialization of matrix-vector product that uses CSysMatrix class for transposed products
107-
*/
108-
template<class ScalarType>
109-
class CSysMatrixVectorProductTransposed final : public CMatrixVectorProduct<ScalarType> {
110-
private:
111-
const CSysMatrix<ScalarType>& matrix; /*!< \brief pointer to matrix that defines the product. */
112-
CGeometry* geometry; /*!< \brief geometry associated with the matrix. */
113-
const CConfig *config; /*!< \brief config of the problem. */
114-
115-
public:
116-
/*!
117-
* \brief constructor of the class
118-
* \param[in] matrix_ref - matrix reference that will be used to define the products
119-
* \param[in] geometry_ref - geometry associated with the problem
120-
* \param[in] config_ref - config of the problem
121-
*/
122-
inline CSysMatrixVectorProductTransposed(const CSysMatrix<ScalarType> & matrix_ref,
123-
CGeometry *geometry_ref, const CConfig *config_ref) :
124-
matrix(matrix_ref),
125-
geometry(geometry_ref),
126-
config(config_ref) {}
127-
128-
/*!
129-
* \note This class cannot be default constructed as that would leave us with invalid pointers.
130-
*/
131-
CSysMatrixVectorProductTransposed() = delete;
132-
133-
/*!
134-
* \brief operator that defines the CSysMatrix-CSysVector product
135-
* \param[in] u - CSysVector that is being multiplied by the sparse matrix
136-
* \param[out] v - CSysVector that is the result of the product
137-
*/
138-
inline void operator()(const CSysVector<ScalarType> & u, CSysVector<ScalarType> & v) const override {
139-
matrix.MatrixVectorProductTransposed(u, v, geometry, config);
140-
}
141-
};

‎Common/include/linear_algebra/CPastixWrapper.hpp‎

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -167,9 +167,18 @@ class CPastixWrapper
167167
* \param[in] geometry - Geometrical definition of the problem.
168168
* \param[in] config - Definition of the particular problem.
169169
* \param[in] kind_fact - Type of factorization.
170-
* \param[in] transposed - Flag to use the transposed matrix during application of the preconditioner.
171170
*/
172-
void Factorize(CGeometry *geometry, const CConfig *config, unsigned short kind_fact, bool transposed);
171+
void Factorize(CGeometry *geometry, const CConfig *config, unsigned short kind_fact);
172+
173+
/*!
174+
* \brief Request solves with the transposed matrix.
175+
* \param[in] transposed - Yes or no.
176+
*/
177+
void SetTransposedSolve(bool transposed = true) {
178+
using namespace PaStiX;
179+
if (iparm[IPARM_SYM] == API_SYM_NO)
180+
iparm[IPARM_TRANSPOSE_SOLVE] = pastix_int_t(!transposed); // negated due to CSR to CSC copy
181+
}
173182

174183
/*!
175184
* \brief Runs the "solve" task for any rhs/sol with operator []

‎Common/include/linear_algebra/CPreconditioner.hpp‎

Lines changed: 47 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@
2626
* License along with SU2. If not, see <http://www.gnu.org/licenses/>.
2727
*/
2828

29-
3029
#pragma once
3130

3231
#include "../CConfig.hpp"
@@ -59,6 +58,17 @@ class CPreconditioner {
5958
* \brief Generic "preprocessing" hook derived classes may implement to build the preconditioner.
6059
*/
6160
virtual void Build() {}
61+
62+
/*!
63+
* \brief Return true to identify the identity preconditioner, may allow some solvers to be more efficient.
64+
*/
65+
virtual bool IsIdentity() const { return false; }
66+
67+
/*!
68+
* \brief Factory method.
69+
*/
70+
static CPreconditioner* Create(ENUM_LINEAR_SOLVER_PREC kind, CSysMatrix<ScalarType>& jacobian,
71+
CGeometry* geometry, const CConfig* config);
6272
};
6373
template<class ScalarType>
6474
CPreconditioner<ScalarType>::~CPreconditioner() {}
@@ -74,25 +84,22 @@ class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
7484
CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
7585
CGeometry* geometry; /*!< \brief Pointer to geometry associated with the matrix. */
7686
const CConfig *config; /*!< \brief Pointer to problem configuration. */
77-
bool transp; /*!< \brief If the transpose version of the preconditioner is required. */
7887

7988
public:
8089
/*!
8190
* \brief Constructor of the class.
8291
* \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
8392
* \param[in] geometry_ref - Geometry associated with the problem.
8493
* \param[in] config_ref - Config of the problem.
85-
* \param[in] transposed - If the transpose version of the preconditioner is required.
8694
*/
8795
inline CJacobiPreconditioner(CSysMatrix<ScalarType> & matrix_ref,
88-
CGeometry *geometry_ref, const CConfig *config_ref, bool transposed) :
96+
CGeometry *geometry_ref, const CConfig *config_ref) :
8997
sparse_matrix(matrix_ref)
9098
{
9199
if((geometry_ref == nullptr) || (config_ref == nullptr))
92100
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
93101
geometry = geometry_ref;
94102
config = config_ref;
95-
transp = transposed;
96103
}
97104

98105
/*!
@@ -113,7 +120,7 @@ class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
113120
* \note Request the associated matrix to build the preconditioner.
114121
*/
115122
inline void Build() override {
116-
sparse_matrix.BuildJacobiPreconditioner(transp);
123+
sparse_matrix.BuildJacobiPreconditioner();
117124
}
118125
};
119126

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

133139
public:
134140
/*!
135141
* \brief Constructor of the class.
136142
* \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
137143
* \param[in] geometry_ref - Geometry associated with the problem.
138144
* \param[in] config_ref - Config of the problem.
139-
* \param[in] transposed - If the transpose version of the preconditioner is required.
140145
*/
141146
inline CILUPreconditioner(CSysMatrix<ScalarType> & matrix_ref,
142-
CGeometry *geometry_ref, const CConfig *config_ref, bool transposed) :
147+
CGeometry *geometry_ref, const CConfig *config_ref) :
143148
sparse_matrix(matrix_ref)
144149
{
145150
if((geometry_ref == nullptr) || (config_ref == nullptr))
146151
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
147152
geometry = geometry_ref;
148153
config = config_ref;
149-
transp = transposed;
150154
}
151155

152156
/*!
@@ -167,7 +171,7 @@ class CILUPreconditioner final : public CPreconditioner<ScalarType> {
167171
* \note Request the associated matrix to build the preconditioner.
168172
*/
169173
inline void Build() override {
170-
sparse_matrix.BuildILUPreconditioner(transp);
174+
sparse_matrix.BuildILUPreconditioner();
171175
}
172176
};
173177

@@ -263,7 +267,7 @@ class CLineletPreconditioner final : public CPreconditioner<ScalarType> {
263267
* \note Request the associated matrix to build the preconditioner.
264268
*/
265269
inline void Build() override {
266-
sparse_matrix.BuildJacobiPreconditioner(false);
270+
sparse_matrix.BuildJacobiPreconditioner();
267271
}
268272
};
269273

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

284287
public:
285288
/*!
@@ -288,18 +291,16 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
288291
* \param[in] geometry_ref - Associated geometry.
289292
* \param[in] config_ref - Problem configuration.
290293
* \param[in] kind_factorization - Type of factorization required.
291-
* \param[in] transposed - If the transpose version of the preconditioner is required.
292294
*/
293295
inline CPastixPreconditioner(CSysMatrix<ScalarType> & matrix_ref, CGeometry *geometry_ref,
294-
const CConfig *config_ref, unsigned short kind_factorization, bool transposed) :
296+
const CConfig *config_ref, unsigned short kind_factorization) :
295297
sparse_matrix(matrix_ref)
296298
{
297299
if((geometry_ref == nullptr) || (config_ref == nullptr))
298300
SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
299301
geometry = geometry_ref;
300302
config = config_ref;
301303
kind_fact = kind_factorization;
302-
transp = transposed;
303304
}
304305

305306
/*!
@@ -320,6 +321,35 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
320321
* \note Request the associated matrix to build the preconditioner.
321322
*/
322323
inline void Build() override {
323-
sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact, transp);
324+
sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact);
324325
}
325326
};
327+
328+
329+
template<class ScalarType>
330+
CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOLVER_PREC kind,
331+
CSysMatrix<ScalarType>& jacobian,
332+
CGeometry* geometry,
333+
const CConfig* config) {
334+
CPreconditioner<ScalarType>* prec = nullptr;
335+
336+
switch (kind) {
337+
case JACOBI:
338+
prec = new CJacobiPreconditioner<ScalarType>(jacobian, geometry, config);
339+
break;
340+
case LINELET:
341+
prec = new CLineletPreconditioner<ScalarType>(jacobian, geometry, config);
342+
break;
343+
case LU_SGS:
344+
prec = new CLU_SGSPreconditioner<ScalarType>(jacobian, geometry, config);
345+
break;
346+
case ILU:
347+
prec = new CILUPreconditioner<ScalarType>(jacobian, geometry, config);
348+
break;
349+
case PASTIX_ILU: case PASTIX_LU_P: case PASTIX_LDLT_P:
350+
prec = new CPastixPreconditioner<ScalarType>(jacobian, geometry, config, kind);
351+
break;
352+
}
353+
354+
return prec;
355+
}

‎Common/include/linear_algebra/CSysMatrix.hpp‎

Lines changed: 17 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -164,8 +164,6 @@ class CSysMatrix {
164164
gemm_t MatrixVectorProductKernelBetaOne; /*!< \brief MKL JIT based GEMV kernel with BETA=1.0. */
165165
void * MatrixVectorProductJitterAlphaMinusOne; /*!< \brief Jitter handle for MKL JIT based GEMV with ALPHA=-1.0 and BETA=1.0. */
166166
gemm_t MatrixVectorProductKernelAlphaMinusOne; /*!< \brief MKL JIT based GEMV kernel with ALPHA=-1.0 and BETA=1.0. */
167-
void * MatrixVectorProductTranspJitterBetaOne; /*!< \brief Jitter handle for MKL JIT based GEMV (transposed) with BETA=1.0. */
168-
gemm_t MatrixVectorProductTranspKernelBetaOne; /*!< \brief MKL JIT based GEMV (transposed) kernel with BETA=1.0. */
169167
#endif
170168

171169
#ifdef HAVE_PASTIX
@@ -177,6 +175,9 @@ class CSysMatrix {
177175
*/
178176
struct {
179177
const unsigned long *ptr = nullptr;
178+
unsigned long nEdge = 0;
179+
180+
operator bool() { return nEdge != 0; }
180181

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

228-
/*!
229-
* \brief Calculates the matrix-vector product: product += matrix^T * vector
230-
* \param[in] matrix
231-
* \param[in] vector
232-
* \param[in,out] product
233-
*/
234-
void MatrixVectorProductTransp(const ScalarType *matrix, const ScalarType *vector, ScalarType *product) const;
235-
236229
/*!
237230
* \brief Calculates the matrix-matrix product
238231
*/
@@ -258,17 +251,10 @@ class CSysMatrix {
258251
/*!
259252
* \brief Copy matrix src into dst, transpose if required.
260253
*/
261-
FORCEINLINE void MatrixCopy(const ScalarType *src, ScalarType *dst, bool transposed = false) const {
262-
if (!transposed) {
263-
SU2_OMP_SIMD
264-
for(auto iVar = 0ul; iVar < nVar*nEqn; ++iVar)
265-
dst[iVar] = src[iVar];
266-
}
267-
else {
268-
for (auto iVar = 0ul; iVar < nVar; ++iVar)
269-
for (auto jVar = 0ul; jVar < nVar; ++jVar)
270-
dst[iVar*nVar+jVar] = src[jVar*nVar+iVar];
271-
}
254+
FORCEINLINE void MatrixCopy(const ScalarType *src, ScalarType *dst) const {
255+
SU2_OMP_SIMD
256+
for(auto iVar = 0ul; iVar < nVar*nEqn; ++iVar)
257+
dst[iVar] = src[iVar];
272258
}
273259

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

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

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

326-
/*!
327-
* \brief Set the transposed value of a block in the sparse matrix.
328-
* \param[in] block_i - Indexes of the block in the matrix-by-blocks structure.
329-
* \param[in] block_j - Indexes of the block in the matrix-by-blocks structure.
330-
* \param[in] **val_block - Block to set to A(i, j).
331-
*/
332-
inline void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block);
333-
334311
/*!
335312
* \brief Performs the product of i-th row of the upper part of a sparse matrix by a vector.
336313
* \param[in] vec - Vector to be multiplied by the upper part of the sparse matrix A.
@@ -805,6 +782,11 @@ class CSysMatrix {
805782
*/
806783
void SetDiagonalAsColumnSum();
807784

785+
/*!
786+
* \brief Transposes the matrix, any preconditioner that was computed may be invalid.
787+
*/
788+
void TransposeInPlace();
789+
808790
/*!
809791
* \brief Add a scaled sparse matrix to "this" (axpy-type operation, A = A+alpha*B).
810792
* \note Matrices must have the same sparse pattern.
@@ -823,20 +805,10 @@ class CSysMatrix {
823805
void MatrixVectorProduct(const CSysVector<ScalarType> & vec, CSysVector<ScalarType> & prod,
824806
CGeometry *geometry, const CConfig *config) const;
825807

826-
/*!
827-
* \brief Performs the product of a sparse matrix by a CSysVector.
828-
* \param[in] vec - CSysVector to be multiplied by the sparse matrix A.
829-
* \param[in] geometry - Geometrical definition of the problem.
830-
* \param[in] config - Definition of the particular problem.
831-
* \param[out] prod - Result of the product.
832-
*/
833-
void MatrixVectorProductTransposed(const CSysVector<ScalarType> & vec, CSysVector<ScalarType> & prod,
834-
CGeometry *geometry, const CConfig *config) const;
835-
836808
/*!
837809
* \brief Build the Jacobi preconditioner.
838810
*/
839-
void BuildJacobiPreconditioner(bool transpose = false);
811+
void BuildJacobiPreconditioner();
840812

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

851823
/*!
852824
* \brief Build the ILU preconditioner.
853-
* \param[in] transposed - Flag to use the transposed matrix to construct the preconditioner.
854825
*/
855-
void BuildILUPreconditioner(bool transposed = false);
826+
void BuildILUPreconditioner();
856827

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

909879
/*!
910880
* \brief Apply the PaStiX factorization to CSysVec.

0 commit comments

Comments
 (0)