diff --git a/Common/include/grid_movement_structure.hpp b/Common/include/grid_movement_structure.hpp index 7d4fa1ffe508..ec631e5a7cb0 100644 --- a/Common/include/grid_movement_structure.hpp +++ b/Common/include/grid_movement_structure.hpp @@ -968,10 +968,10 @@ class CVolumetricMovement : public CGridMovement { unsigned long nIterMesh; /*!< \brief Number of iterations in the mesh update. +*/ - CSysSolve System; - CSysMatrix StiffMatrix; /*!< \brief Matrix to store the point-to-point stiffness. */ - CSysVector LinSysSol; - CSysVector LinSysRes; + CSysSolve System; + CSysMatrix StiffMatrix; /*!< \brief Matrix to store the point-to-point stiffness. */ + CSysVector LinSysSol; + CSysVector LinSysRes; public: @@ -1331,10 +1331,15 @@ class CElasticityMovement : public CVolumetricMovement { su2double MinVolume; su2double MaxVolume; - CSysSolve System; - CSysMatrix StiffMatrix; /*!< \brief Matrix to store the point-to-point stiffness. */ - CSysVector LinSysSol; - CSysVector LinSysRes; +#ifndef CODI_FORWARD_TYPE + CSysSolve System; + CSysMatrix StiffMatrix; /*!< \brief Matrix to store the point-to-point stiffness. */ +#else + CSysSolve System; + CSysMatrix StiffMatrix; +#endif + CSysVector LinSysSol; + CSysVector LinSysRes; su2double E; /*!< \brief Young's modulus of elasticity. */ su2double Nu; /*!< \brief Poisson's ratio. */ diff --git a/Common/include/linear_solvers_structure.hpp b/Common/include/linear_solvers_structure.hpp index 5a2273af99a6..2133a6bf26c0 100644 --- a/Common/include/linear_solvers_structure.hpp +++ b/Common/include/linear_solvers_structure.hpp @@ -67,27 +67,40 @@ using namespace std; * matrix-vector products and preconditioners to different problems * that may arise in a hierarchical solver (i.e. multigrid). */ +template class CSysSolve { - + +public: + /*--- Some typedefs for simplicity ---*/ + typedef CSysVector VectorType; + typedef CSysMatrix MatrixType; + typedef CMatrixVectorProduct ProductType; + typedef CPreconditioner PrecondType; + private: - bool mesh_deform; /*!< \brief Operate in mesh deformation mode, changes the source of solver options. */ - su2double Residual;/*!< \brief Residual at the end of a call to Solve. */ + bool mesh_deform; /*!< \brief Operate in mesh deformation mode, changes the source of solver options. */ + ScalarType Residual; /*!< \brief Residual at the end of a call to Solve. */ bool cg_ready; /*!< \brief Indicate if memory used by CG is allocated. */ bool bcg_ready; /*!< \brief Indicate if memory used by BCGSTAB is allocated. */ bool gmres_ready; /*!< \brief Indicate if memory used by FGMRES is allocated. */ - CSysVector r; /*!< \brief Residual in CG and BCGSTAB. */ - CSysVector A_x; /*!< \brief Result of matrix-vector product in CG and BCGSTAB. */ - CSysVector p; /*!< \brief Direction in CG and BCGSTAB. */ - CSysVector z; /*!< \brief Preconditioned residual/direction in CG/BCGSTAB. */ + VectorType r; /*!< \brief Residual in CG and BCGSTAB. */ + VectorType A_x; /*!< \brief Result of matrix-vector product in CG and BCGSTAB. */ + VectorType p; /*!< \brief Direction in CG and BCGSTAB. */ + VectorType z; /*!< \brief Preconditioned residual/direction in CG/BCGSTAB. */ - CSysVector r_0; /*!< \brief The "arbitrary" vector in BCGSTAB. */ - CSysVector v; /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */ + VectorType r_0; /*!< \brief The "arbitrary" vector in BCGSTAB. */ + VectorType v; /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */ - vector W; /*!< \brief Large matrix used by FGMRES, w^i+1 = A * z^i. */ - vector Z; /*!< \brief Large matrix used by FGMRES, preconditioned W. */ + vector W; /*!< \brief Large matrix used by FGMRES, w^i+1 = A * z^i. */ + vector Z; /*!< \brief Large matrix used by FGMRES, preconditioned W. */ + + VectorType LinSysRes_tmp; /*!< \brief Temporary used when it is necessary to interface between active and passive types. */ + VectorType LinSysSol_tmp; /*!< \brief Temporary used when it is necessary to interface between active and passive types. */ + VectorType* LinSysRes_ptr; /*!< \brief Pointer to appropriate LinSysRes (set to original or temporary in call to Solve). */ + VectorType* LinSysSol_ptr; /*!< \brief Pointer to appropriate LinSysSol (set to original or temporary in call to Solve). */ /*! * \brief sign transfer function @@ -98,28 +111,28 @@ class CSysSolve { * so, feel free to delete this and replace it as needed with the * appropriate global function */ - su2double Sign(const su2double & x, const su2double & y) const; + ScalarType Sign(const ScalarType & x, const ScalarType & y) const; /*! * \brief applys a Givens rotation to a 2-vector * \param[in] s - sine of the Givens rotation angle * \param[in] c - cosine of the Givens rotation angle - * \param[in, out] h1 - first element of 2x1 vector being transformed - * \param[in, out] h2 - second element of 2x1 vector being transformed + * \param[in,out] h1 - first element of 2x1 vector being transformed + * \param[in,out] h2 - second element of 2x1 vector being transformed */ - void ApplyGivens(const su2double & s, const su2double & c, su2double & h1, su2double & h2); + void ApplyGivens(const ScalarType & s, const ScalarType & c, ScalarType & h1, ScalarType & h2); /*! * \brief generates the Givens rotation matrix for a given 2-vector - * \param[in, out] dx - element of 2x1 vector being transformed - * \param[in, out] dy - element of 2x1 vector being set to zero - * \param[in, out] s - sine of the Givens rotation angle - * \param[in, out] c - cosine of the Givens rotation angle + * \param[in,out] dx - element of 2x1 vector being transformed + * \param[in,out] dy - element of 2x1 vector being set to zero + * \param[in,out] s - sine of the Givens rotation angle + * \param[in,out] c - cosine of the Givens rotation angle * * Based on givens() of SPARSKIT, which is based on p.202 of * "Matrix Computations" by Golub and van Loan. */ - void GenerateGivens(su2double & dx, su2double & dy, su2double & s, su2double & c); + void GenerateGivens(ScalarType & dx, ScalarType & dy, ScalarType & s, ScalarType & c); /*! * \brief finds the solution of the upper triangular system Hsbg*x = rhs @@ -132,8 +145,8 @@ class CSysSolve { * \pre the upper Hessenberg matrix has been transformed into a * triangular matrix. */ - void SolveReduced(const int & n, const vector > & Hsbg, - const vector & rhs, vector & x); + void SolveReduced(const int & n, const vector > & Hsbg, + const vector & rhs, vector & x); /*! * \brief Modified Gram-Schmidt orthogonalization @@ -153,7 +166,7 @@ class CSysSolve { * vector is kept in nrm0 and updated after operating with each vector * */ - void ModGramSchmidt(int i, vector > & Hsbg, vector & w); + void ModGramSchmidt(int i, vector > & Hsbg, vector & w); /*! * \brief writes header information for a CSysSolve residual history @@ -163,7 +176,7 @@ class CSysSolve { * * \pre the ostream object os should be open */ - void WriteHeader(const string & solver, const su2double & restol, const su2double & resinit); + void WriteHeader(const string & solver, const ScalarType & restol, const ScalarType & resinit); /*! * \brief writes residual convergence data for one iteration to a stream @@ -173,8 +186,21 @@ class CSysSolve { * * \pre the ostream object os should be open */ - void WriteHistory(const int & iter, const su2double & res, const su2double & resinit); + void WriteHistory(const int & iter, const ScalarType & res, const ScalarType & resinit); + + /*! + * \brief Used by Solve for compatibility between passive and active CSysVector, see specializations. + * \param[in] LinSysRes - Linear system residual + * \param[in,out] LinSysSol - Linear system solution + */ + void HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol); + /*! + * \brief Used by Solve for compatibility between passive and active CSysVector, see specializations. + * \param[out] LinSysSol - Linear system solution + */ + void HandleTemporariesOut(CSysVector & LinSysSol); + public: /*! @@ -193,9 +219,9 @@ class CSysSolve { * \param[in] monitoring - turn on priting residuals from solver to screen. * \param[in] config - Definition of the particular problem. */ - unsigned long CG_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, - unsigned long m, su2double *residual, bool monitoring, CConfig *config); + unsigned long CG_LinSolver(const VectorType & b, VectorType & x, ProductType & mat_vec, + PrecondType & precond, ScalarType tol, unsigned long m, + ScalarType *residual, bool monitoring, CConfig *config); /*! * \brief Flexible Generalized Minimal Residual method @@ -209,9 +235,9 @@ class CSysSolve { * \param[in] monitoring - turn on priting residuals from solver to screen. * \param[in] config - Definition of the particular problem. */ - unsigned long FGMRES_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, - unsigned long m, su2double *residual, bool monitoring, CConfig *config); + unsigned long FGMRES_LinSolver(const VectorType & b, VectorType & x, ProductType & mat_vec, + PrecondType & precond, ScalarType tol, unsigned long m, + ScalarType *residual, bool monitoring, CConfig *config); /*! * \brief Biconjugate Gradient Stabilized Method (BCGSTAB) @@ -225,9 +251,9 @@ class CSysSolve { * \param[in] monitoring - turn on priting residuals from solver to screen. * \param[in] config - Definition of the particular problem. */ - unsigned long BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, - unsigned long m, su2double *residual, bool monitoring, CConfig *config); + unsigned long BCGSTAB_LinSolver(const VectorType & b, VectorType & x, ProductType & mat_vec, + PrecondType & precond, ScalarType tol, unsigned long m, + ScalarType *residual, bool monitoring, CConfig *config); /*! * \brief Solve the linear system using a Krylov subspace method @@ -237,7 +263,8 @@ class CSysSolve { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - unsigned long Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config); + unsigned long Solve(MatrixType & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, + CGeometry *geometry, CConfig *config); /*! * \brief Solve the adjoint linear system using a Krylov subspace method @@ -247,13 +274,14 @@ class CSysSolve { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - unsigned long Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config); - + unsigned long Solve_b(MatrixType & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, + CGeometry *geometry, CConfig *config); + /*! * \brief Get the final residual. * \return The residual at the end of Solve */ - su2double GetResidual(void) const; + ScalarType GetResidual(void) const; }; diff --git a/Common/include/linear_solvers_structure.inl b/Common/include/linear_solvers_structure.inl index 0597ec571f20..b229828027df 100644 --- a/Common/include/linear_solvers_structure.inl +++ b/Common/include/linear_solvers_structure.inl @@ -37,7 +37,8 @@ #pragma once -inline su2double CSysSolve::Sign(const su2double & x, const su2double & y) const { +template +inline ScalarType CSysSolve::Sign(const ScalarType & x, const ScalarType & y) const { if (y == 0.0) return 0.0; else { @@ -47,4 +48,11 @@ inline su2double CSysSolve::Sign(const su2double & x, const su2double & y) const } } -inline su2double CSysSolve::GetResidual(void) const { return Residual; } +template +inline ScalarType CSysSolve::GetResidual(void) const { return Residual; } + +template +void CSysSolve::HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol) {} + +template +void CSysSolve::HandleTemporariesOut(CSysVector & LinSysSol) {} diff --git a/Common/include/linear_solvers_structure_b.hpp b/Common/include/linear_solvers_structure_b.hpp index 877b1f758645..eb2a645f2391 100644 --- a/Common/include/linear_solvers_structure_b.hpp +++ b/Common/include/linear_solvers_structure_b.hpp @@ -40,9 +40,12 @@ #include "config_structure.hpp" #ifdef CODI_REVERSE_TYPE +template class CSysSolve_b{ public: - static void Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m, const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n, codi::DataStore* d); + static void Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m, + const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n, + codi::DataStore* d); }; #endif diff --git a/Common/include/matrix_structure.hpp b/Common/include/matrix_structure.hpp index 756ec439daa2..b538e0970b97 100644 --- a/Common/include/matrix_structure.hpp +++ b/Common/include/matrix_structure.hpp @@ -63,6 +63,7 @@ const su2double eps = numeric_limits::epsilon(); /*!< \brief mach with compressed row format. * \author A. Bueno, F. Palacios */ +template class CSysMatrix { private: int rank, /*!< \brief MPI Rank. */ @@ -71,8 +72,8 @@ class CSysMatrix { nPointDomain, /*!< \brief Number of points in the grid. */ nVar, /*!< \brief Number of variables. */ nEqn; /*!< \brief Number of equations. */ - su2double *matrix; /*!< \brief Entries of the sparse matrix. */ - su2double *ILU_matrix; /*!< \brief Entries of the ILU sparse matrix. */ + ScalarType *matrix; /*!< \brief Entries of the sparse matrix. */ + ScalarType *ILU_matrix; /*!< \brief Entries of the ILU sparse matrix. */ unsigned long nnz; /*!< \brief Number of possible nonzero entries in the matrix. */ unsigned long *row_ptr; /*!< \brief Pointers to the first element in each row. */ unsigned long *col_ind; /*!< \brief Column index for each of the elements in val(). */ @@ -81,19 +82,19 @@ class CSysMatrix { unsigned long *col_ind_ilu; /*!< \brief Column index for each of the elements in val() (ILU). */ unsigned short ilu_fill_in; /*!< \brief Fill in level for the ILU preconditioner. */ - su2double *block; /*!< \brief Internal array to store a subblock of the matrix. */ - su2double *block_inverse; /*!< \brief Internal array to store a subblock of the matrix. */ - su2double *block_weight; /*!< \brief Internal array to store a subblock of the matrix. */ - su2double *prod_block_vector; /*!< \brief Internal array to store the product of a subblock with a vector. */ - su2double *prod_row_vector; /*!< \brief Internal array to store the product of a matrix-by-blocks "row" with a vector. */ - su2double *aux_vector; /*!< \brief Auxiliary array to store intermediate results. */ - su2double *sum_vector; /*!< \brief Auxiliary array to store intermediate results. */ - su2double *invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ + ScalarType *block; /*!< \brief Internal array to store a subblock of the matrix. */ + ScalarType *block_inverse; /*!< \brief Internal array to store a subblock of the matrix. */ + ScalarType *block_weight; /*!< \brief Internal array to store a subblock of the matrix. */ + ScalarType *prod_block_vector; /*!< \brief Internal array to store the product of a subblock with a vector. */ + ScalarType *prod_row_vector; /*!< \brief Internal array to store the product of a matrix-by-blocks "row" with a vector. */ + ScalarType *aux_vector; /*!< \brief Auxiliary array to store intermediate results. */ + ScalarType *sum_vector; /*!< \brief Auxiliary array to store intermediate results. */ + ScalarType *invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ bool *LineletBool; /*!< \brief Identify if a point belong to a linelet. */ vector *LineletPoint; /*!< \brief Linelet structure. */ unsigned long nLinelet; /*!< \brief Number of Linelets in the system. */ - su2double **UBlock, **invUBlock, **LBlock, + ScalarType **UBlock, **invUBlock, **LBlock, **yVector, **zVector, **rVector, *LFBlock, *LyVector, *FzVector; /*!< \brief Arrays of the Linelet preconditioner methodology. */ unsigned long max_nElem; @@ -108,6 +109,18 @@ class CSysMatrix { bool useMKL; #endif + /*! + * \brief Handle type conversion for when we Set, Add, etc. blocks, preserving derivative information (if supported by types). + */ + template + DstType ActiveAssign(const SrcType & val) const; + + /*! + * \brief Handle type conversion for when we Set, Add, etc. blocks, discarding derivative information. + */ + template + DstType PassiveAssign(const SrcType & val) const; + public: /*! @@ -165,7 +178,8 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[in] commType - Enumerated type for the quantity to be communicated. */ - void InitiateComms(CSysVector & x, + template + void InitiateComms(CSysVector & x, CGeometry *geometry, CConfig *config, unsigned short commType); @@ -177,7 +191,8 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[in] commType - Enumerated type for the quantity to be unpacked. */ - void CompleteComms(CSysVector & x, + template + void CompleteComms(CSysVector & x, CGeometry *geometry, CConfig *config, unsigned short commType); @@ -187,14 +202,14 @@ class CSysMatrix { * \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. */ - su2double *GetBlock(unsigned long block_i, unsigned long block_j); + ScalarType *GetBlock(unsigned long block_i, unsigned long block_j); /*! * \brief Copies the block (i, j) of the matrix-by-blocks structure in the internal variable *block. * \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. */ - su2double GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar); + ScalarType GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar); /*! * \brief Set the value of a block in the sparse matrix. @@ -202,7 +217,8 @@ class CSysMatrix { * \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). */ - void SetBlock(unsigned long block_i, unsigned long block_j, su2double **val_block); + template + void SetBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); /*! * \brief Set the value of a block in the sparse matrix. @@ -210,7 +226,8 @@ class CSysMatrix { * \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). */ - void SetBlock(unsigned long block_i, unsigned long block_j, su2double *val_block); + template + void SetBlock(unsigned long block_i, unsigned long block_j, OtherType *val_block); /*! * \brief Adds the specified block to the sparse matrix. @@ -218,7 +235,8 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to add to A(i, j). */ - void AddBlock(unsigned long block_i, unsigned long block_j, su2double **val_block); + template + void AddBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); /*! * \brief Subtracts the specified block to the sparse matrix. @@ -226,14 +244,15 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to subtract to A(i, j). */ - void SubtractBlock(unsigned long block_i, unsigned long block_j, su2double **val_block); + template + void SubtractBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); /*! * \brief Copies the block (i, j) of the matrix-by-blocks structure in the internal variable *block. * \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. */ - su2double *GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j); + ScalarType *GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j); /*! * \brief Set the value of a block in the sparse matrix. @@ -241,7 +260,7 @@ class CSysMatrix { * \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). */ - void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block); + void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! @@ -250,7 +269,7 @@ class CSysMatrix { * \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). */ - void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block); + void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! * \brief Subtracts the specified block to the sparse matrix. @@ -258,7 +277,7 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to subtract to A(i, j). */ - void SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block); + void SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! * \brief Adds the specified value to the diagonal of the (i, i) subblock @@ -266,7 +285,8 @@ class CSysMatrix { * \param[in] block_i - Index of the block in the matrix-by-blocks structure. * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). */ - void AddVal2Diag(unsigned long block_i, su2double val_matrix); + template + void AddVal2Diag(unsigned long block_i, OtherType val_matrix); /*! * \brief Sets the specified value to the diagonal of the (i, i) subblock @@ -274,7 +294,8 @@ class CSysMatrix { * \param[in] block_i - Index of the block in the matrix-by-blocks structure. * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). */ - void SetVal2Diag(unsigned long block_i, su2double val_matrix); + template + void SetVal2Diag(unsigned long block_i, OtherType val_matrix); /*! * \brief Calculates the matrix-vector product @@ -282,7 +303,7 @@ class CSysMatrix { * \param[in] vector * \param[out] product */ - void MatrixVectorProduct(su2double *matrix, su2double *vector, su2double *product); + void MatrixVectorProduct(ScalarType *matrix, ScalarType *vector, ScalarType *product); /*! * \brief Calculates the matrix-matrix product @@ -290,7 +311,7 @@ class CSysMatrix { * \param[in] matrix_b * \param[out] product */ - void MatrixMatrixProduct(su2double *matrix_a, su2double *matrix_b, su2double *product); + void MatrixMatrixProduct(ScalarType *matrix_a, ScalarType *matrix_b, ScalarType *product); /*! * \brief Deletes the values of the row i of the sparse matrix. @@ -304,7 +325,7 @@ class CSysMatrix { * \param[in] n - Size of the quare matrix. * \return Value of the determinant. */ - su2double MatrixDeterminant(su2double **a, unsigned long n); + ScalarType MatrixDeterminant(ScalarType **a, unsigned long n); /*! * \brief Find the cofactor matrix of a square matrix. Written by Paul Bourke @@ -312,14 +333,14 @@ class CSysMatrix { * \param[in] n - Size of the quare matrix. * \param[out] b - cofactor matrix */ - void MatrixCoFactor(su2double **a, unsigned long n, su2double **b) ; + void MatrixCoFactor(ScalarType **a, unsigned long n, ScalarType **b) ; /*! * \brief Transpose of a square matrix, do it in place. Written by Paul Bourke * \param[in] a - Matrix to compute the determinant. * \param[in] n - Size of the quare matrix. */ - void MatrixTranspose(su2double **a, unsigned long n) ; + void MatrixTranspose(ScalarType **a, unsigned long n) ; /*! * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. @@ -328,7 +349,7 @@ class CSysMatrix { * \param[in] transposed - If true the transposed of the block is used (default = false). * \return Solution of the linear system (overwritten on rhs). */ - void Gauss_Elimination(unsigned long block_i, su2double* rhs, bool transposed = false); + void Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed = false); /*! * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. @@ -336,7 +357,7 @@ class CSysMatrix { * \param[in] rhs - Right-hand-side of the linear system. * \return Solution of the linear system (overwritten on rhs). */ - void Gauss_Elimination(su2double* Block, su2double* rhs); + void Gauss_Elimination(ScalarType* Block, ScalarType* rhs); /*! * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. @@ -344,7 +365,7 @@ class CSysMatrix { * \param[in] rhs - Right-hand-side of the linear system. * \return Solution of the linear system (overwritten on rhs). */ - void Gauss_Elimination_ILUMatrix(unsigned long block_i, su2double* rhs); + void Gauss_Elimination_ILUMatrix(unsigned long block_i, ScalarType* rhs); /*! * \fn void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, su2double* vec); @@ -354,7 +375,7 @@ class CSysMatrix { * \param[in] vec - Vector to be multiplied by the block (i, j) of the sparse matrix A. * \return Product of A(i, j) by vector *vec (stored at *prod_block_vector). */ - void ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec); + void ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec); /*! * \brief Performs the product of i-th row of the upper part of a sparse matrix by a vector. @@ -362,7 +383,7 @@ class CSysMatrix { * \param[in] row_i - Row of the matrix to be multiplied by vector vec. * \return prod Result of the product U(A)*vec (stored at *prod_row_vector). */ - void UpperProduct(CSysVector & vec, unsigned long row_i); + void UpperProduct(CSysVector & vec, unsigned long row_i); /*! * \brief Performs the product of i-th row of the lower part of a sparse matrix by a vector. @@ -370,7 +391,7 @@ class CSysMatrix { * \param[in] row_i - Row of the matrix to be multiplied by vector vec. * \return prod Result of the product L(A)*vec (stored at *prod_row_vector). */ - void LowerProduct(CSysVector & vec, unsigned long row_i); + void LowerProduct(CSysVector & vec, unsigned long row_i); /*! * \brief Performs the product of i-th row of the diagonal part of a sparse matrix by a vector. @@ -378,7 +399,7 @@ class CSysMatrix { * \param[in] row_i - Row of the matrix to be multiplied by vector vec. * \return prod Result of the product D(A)*vec (stored at *prod_row_vector). */ - void DiagonalProduct(CSysVector & vec, unsigned long row_i); + void DiagonalProduct(CSysVector & vec, unsigned long row_i); /*! * \brief Send receive the solution using MPI. @@ -386,7 +407,8 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config); + template + void SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config); /*! * \brief Send receive the solution using MPI and the transposed structure of the matrix. @@ -394,7 +416,7 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config); + void SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config); /*! * \brief Performs the product of i-th row of a sparse matrix by a vector. @@ -402,7 +424,7 @@ class CSysMatrix { * \param[in] row_i - Row of the matrix to be multiplied by vector vec. * \return Result of the product (stored at *prod_row_vector). */ - void RowProduct(const CSysVector & vec, unsigned long row_i); + void RowProduct(const CSysVector & vec, unsigned long row_i); /*! * \brief Performs the product of a sparse matrix by a vector. @@ -410,7 +432,7 @@ class CSysMatrix { * \param[out] prod - Result of the product. * \return Result of the product A*vec. */ - void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod); + void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod); /*! * \brief Performs the product of a sparse matrix by a CSysVector. @@ -419,7 +441,7 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[out] prod - Result of the product. */ - void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Performs the product of a sparse matrix by a CSysVector. @@ -428,48 +450,48 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[out] prod - Result of the product. */ - void MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Performs the product of two block matrices. */ - void GetMultBlockBlock(su2double *c, su2double *a, su2double *b); + void GetMultBlockBlock(ScalarType *c, ScalarType *a, ScalarType *b); /*! * \brief Performs the product of a block matrices by a vector. */ - void GetMultBlockVector(su2double *c, su2double *a, su2double *b); + void GetMultBlockVector(ScalarType *c, ScalarType *a, ScalarType *b); /*! * \brief Performs the subtraction of two matrices. */ - void GetSubsBlock(su2double *c, su2double *a, su2double *b); + void GetSubsBlock(ScalarType *c, ScalarType *a, ScalarType *b); /*! * \brief Performs the subtraction of two vectors. */ - void GetSubsVector(su2double *c, su2double *a, su2double *b); + void GetSubsVector(ScalarType *c, ScalarType *a, ScalarType *b); /*! * \brief Inverse diagonal block. * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. * \param[out] invBlock - Inverse block. */ - void InverseDiagonalBlock(unsigned long block_i, su2double *invBlock, bool transpose = false); + void InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose = false); /*! * \brief Inverse diagonal block. * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. * \param[out] invBlock - Inverse block. */ - void InverseDiagonalBlock_ILUMatrix(unsigned long block_i, su2double *invBlock); + void InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock); /*! * \brief Inverse a block. * \param[in] Block - block matrix. * \param[out] invBlock - Inverse block. */ - void InverseBlock(su2double *Block, su2double *invBlock); + void InverseBlock(ScalarType *Block, ScalarType *invBlock); /*! * \brief Build the Jacobi preconditioner. @@ -483,7 +505,7 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Apply Jacobi as a classical iterative smoother @@ -498,7 +520,7 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). */ - unsigned long Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config); + unsigned long Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); /*! * \brief Build the ILU preconditioner. @@ -513,7 +535,7 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Apply ILU as a classical iterative smoother @@ -528,14 +550,14 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). */ - unsigned long ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config); + unsigned long ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); /*! * \brief Multiply CSysVector by the preconditioner * \param[in] vec - CSysVector to be multiplied by the preconditioner. * \param[out] prod - Result of the product A*vec. */ - void ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Apply LU_SGS as a classical iterative smoother @@ -550,7 +572,7 @@ class CSysMatrix { * \param[in] config - Definition of the particular problem. * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). */ - unsigned long LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config); + unsigned long LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); /*! * \brief Build the Linelet preconditioner. @@ -564,7 +586,7 @@ class CSysMatrix { * \param[in] vec - CSysVector to be multiplied by the preconditioner. * \param[out] prod - Result of the product A*vec. */ - void ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + void ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Compute the residual Ax-b @@ -572,7 +594,7 @@ class CSysMatrix { * \param[in] f - Result of the product A*vec. * \param[out] res - Result of the product A*vec. */ - void ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res); + void ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res); }; @@ -580,9 +602,10 @@ class CSysMatrix { * \class CSysMatrixVectorProduct * \brief specialization of matrix-vector product that uses CSysMatrix class */ -class CSysMatrixVectorProduct : public CMatrixVectorProduct { +template +class CSysMatrixVectorProduct : public CMatrixVectorProduct { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -594,7 +617,7 @@ class CSysMatrixVectorProduct : public CMatrixVectorProduct { * \param[in] geometry_ref - * \param[in] config_ref - */ - CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -606,16 +629,17 @@ class CSysMatrixVectorProduct : public CMatrixVectorProduct { * \param[in] u - CSysVector that is being multiplied by the sparse matrix * \param[out] v - CSysVector that is the result of the product */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CSysMatrixVectorProduct * \brief specialization of matrix-vector product that uses CSysMatrix class */ -class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { +template +class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -627,7 +651,7 @@ class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { * \param[in] geometry_ref - * \param[in] config_ref - */ - CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -639,16 +663,17 @@ class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { * \param[in] u - CSysVector that is being multiplied by the sparse matrix * \param[out] v - CSysVector that is the result of the product */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CJacobiPreconditioner * \brief specialization of preconditioner that uses CSysMatrix class */ -class CJacobiPreconditioner : public CPreconditioner { +template +class CJacobiPreconditioner : public CPreconditioner { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -660,7 +685,7 @@ class CJacobiPreconditioner : public CPreconditioner { * \param[in] geometry_ref - * \param[in] config_ref - */ - CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -672,16 +697,17 @@ class CJacobiPreconditioner : public CPreconditioner { * \param[in] u - CSysVector that is being preconditioned * \param[out] v - CSysVector that is the result of the preconditioning */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CJacobiTransposedPreconditioner * \brief specialization of preconditioner that uses CSysMatrix class */ -class CJacobiTransposedPreconditioner : public CPreconditioner { +template +class CJacobiTransposedPreconditioner : public CPreconditioner { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -693,7 +719,7 @@ class CJacobiTransposedPreconditioner : public CPreconditioner { * \param[in] geometry_ref - * \param[in] config_ref - */ - CJacobiTransposedPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CJacobiTransposedPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -705,16 +731,17 @@ class CJacobiTransposedPreconditioner : public CPreconditioner { * \param[in] u - CSysVector that is being preconditioned * \param[out] v - CSysVector that is the result of the preconditioning */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CILUPreconditioner * \brief specialization of preconditioner that uses CSysMatrix class */ -class CILUPreconditioner : public CPreconditioner { +template +class CILUPreconditioner : public CPreconditioner { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -726,7 +753,7 @@ class CILUPreconditioner : public CPreconditioner { * \param[in] geometry_ref - * \param[in] config_ref - */ - CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -738,16 +765,17 @@ class CILUPreconditioner : public CPreconditioner { * \param[in] u - CSysVector that is being preconditioned * \param[out] v - CSysVector that is the result of the preconditioning */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CLU_SGSPreconditioner * \brief specialization of preconditioner that uses CSysMatrix class */ -class CLU_SGSPreconditioner : public CPreconditioner { +template +class CLU_SGSPreconditioner : public CPreconditioner { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -759,7 +787,7 @@ class CLU_SGSPreconditioner : public CPreconditioner { * \param[in] geometry_ref - * \param[in] config_ref - */ - CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -771,16 +799,17 @@ class CLU_SGSPreconditioner : public CPreconditioner { * \param[in] u - CSysVector that is being preconditioned * \param[out] v - CSysVector that is the result of the preconditioning */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; /*! * \class CLineletPreconditioner * \brief specialization of preconditioner that uses CSysMatrix class */ -class CLineletPreconditioner : public CPreconditioner { +template +class CLineletPreconditioner : public CPreconditioner { private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ CConfig* config; /*!< \brief pointer to matrix that defines the config. */ @@ -792,7 +821,7 @@ class CLineletPreconditioner : public CPreconditioner { * \param[in] geometry_ref - * \param[in] config_ref - */ - CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); + CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); /*! * \brief destructor of the class @@ -804,7 +833,7 @@ class CLineletPreconditioner : public CPreconditioner { * \param[in] u - CSysVector that is being preconditioned * \param[out] v - CSysVector that is the result of the preconditioning */ - void operator()(const CSysVector & u, CSysVector & v) const; + void operator()(const CSysVector & u, CSysVector & v) const; }; #include "matrix_structure.inl" diff --git a/Common/include/matrix_structure.inl b/Common/include/matrix_structure.inl index 5e4ca11d2029..1d02bad80b6d 100644 --- a/Common/include/matrix_structure.inl +++ b/Common/include/matrix_structure.inl @@ -34,23 +34,160 @@ * You should have received a copy of the GNU Lesser General Public * License along with SU2. If not, see . */ - + #pragma once -inline void CSysMatrix::SetValZero(void) { +template +inline void CSysMatrix::SetValZero(void) { if(NULL != matrix) { for (unsigned long index = 0; index < nnz*nVar*nEqn; index++) matrix[index] = 0.0; } } -inline CSysMatrixVectorProduct::CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +template +inline DstType CSysMatrix::ActiveAssign(const SrcType & val) const { return val; } + +#ifdef CODI_REVERSE_TYPE +template<> template<> +inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } + +template<> template<> +inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } +#endif + +template +template +inline DstType CSysMatrix::PassiveAssign(const SrcType & val) const { +#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) + return SU2_TYPE::GetValue(val); +#else + return val; +#endif +} + +template +template +inline void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index, step = 0; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = PassiveAssign(val_block[iVar][jVar]); + break; + } + } + +} + +template +template +inline void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, OtherType *val_block) { + + unsigned long iVar, jVar, index, step = 0; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = PassiveAssign(val_block[iVar*nVar+jVar]); + break; + } + } + +} + +template +template +inline void CSysMatrix::AddBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index, step = 0; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] += PassiveAssign(val_block[iVar][jVar]); + break; + } + } + +} + +template +template +inline void CSysMatrix::SubtractBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index, step = 0; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] -= PassiveAssign(val_block[iVar][jVar]); + break; + } + } + +} + +template +template +inline void CSysMatrix::AddVal2Diag(unsigned long block_i, OtherType val_matrix) { + + unsigned long step = 0, iVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_i) { // Only elements on the diagonal + for (iVar = 0; iVar < nVar; iVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] += PassiveAssign(val_matrix); + break; + } + } + +} + +template +template +inline void CSysMatrix::SetVal2Diag(unsigned long block_i, OtherType val_matrix) { + + unsigned long step = 0, iVar, jVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + step++; + if (col_ind[index] == block_i) { // Only elements on the diagonal + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+jVar] = 0.0; + + for (iVar = 0; iVar < nVar; iVar++) + matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] = PassiveAssign(val_matrix); + + break; + } + } + +} + +template +inline CSysMatrixVectorProduct::CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CSysMatrixVectorProduct::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CSysMatrixVectorProduct::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CSysMatrixVectorProduct::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; @@ -59,13 +196,15 @@ inline void CSysMatrixVectorProduct::operator()(const CSysVector & u, CSysVector sparse_matrix->MatrixVectorProduct(u, v, geometry, config); } -inline CSysMatrixVectorProductTransposed::CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +inline CSysMatrixVectorProductTransposed::CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CSysMatrixVectorProductTransposed::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CSysMatrixVectorProductTransposed::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CSysMatrixVectorProduct::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; @@ -74,14 +213,15 @@ inline void CSysMatrixVectorProductTransposed::operator()(const CSysVector & u, sparse_matrix->MatrixVectorProductTransposed(u, v, geometry, config); } - -inline CJacobiPreconditioner::CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +inline CJacobiPreconditioner::CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CJacobiPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CJacobiPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CJacobiPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; @@ -90,13 +230,15 @@ inline void CJacobiPreconditioner::operator()(const CSysVector & u, CSysVector & sparse_matrix->ComputeJacobiPreconditioner(u, v, geometry, config); } -inline CILUPreconditioner::CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +inline CILUPreconditioner::CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CILUPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CILUPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CILUPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; @@ -105,13 +247,15 @@ inline void CILUPreconditioner::operator()(const CSysVector & u, CSysVector & v) sparse_matrix->ComputeILUPreconditioner(u, v, geometry, config); } -inline CLU_SGSPreconditioner::CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +inline CLU_SGSPreconditioner::CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CLU_SGSPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CLU_SGSPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CLU_SGSPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; @@ -120,13 +264,15 @@ inline void CLU_SGSPreconditioner::operator()(const CSysVector & u, CSysVector & sparse_matrix->ComputeLU_SGSPreconditioner(u, v, geometry, config); } -inline CLineletPreconditioner::CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { +template +inline CLineletPreconditioner::CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { sparse_matrix = &matrix_ref; geometry = geometry_ref; config = config_ref; } -inline void CLineletPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { +template +inline void CLineletPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { if (sparse_matrix == NULL) { cerr << "CLineletPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; cerr << "pointer to sparse matrix is NULL." << endl; diff --git a/Common/include/mpi_structure.hpp b/Common/include/mpi_structure.hpp index af8a68cd5978..b9a2b89b2702 100644 --- a/Common/include/mpi_structure.hpp +++ b/Common/include/mpi_structure.hpp @@ -77,6 +77,18 @@ class CBaseMPIWrapper; typedef CBaseMPIWrapper SU2_MPI; #endif // defined CODI_REVERSE_TYPE || defined CODI_FORWARD_TYPE +/*--- Select the appropriate MPI wrapper based on datatype, to use in templated classes. ---*/ +template struct SelectMPIWrapper {}; + +/*--- This one is the default in both direct and AD. ---*/ +template<> struct SelectMPIWrapper { typedef SU2_MPI W; }; + +/*--- In AD we overload also for the passive wrapper. ---*/ +#if defined CODI_REVERSE_TYPE +class CBaseMPIWrapper; +template<> struct SelectMPIWrapper { typedef CBaseMPIWrapper W; }; +#endif + /*! * \class CMPIWrapper * \brief Class for defining the MPI wrapper routines; this class features as a base class for diff --git a/Common/include/vector_structure.hpp b/Common/include/vector_structure.hpp index cfea994acce6..edb26337ca08 100644 --- a/Common/include/vector_structure.hpp +++ b/Common/include/vector_structure.hpp @@ -50,6 +50,10 @@ using namespace std; +/*--- Forward declaration of template friend functions. ---*/ +template class CSysVector; +template T dotProd(const CSysVector & u, const CSysVector & v); + /*! * \class CSysVector * \brief Class for holding and manipulating vectors needed by linear solvers @@ -60,6 +64,7 @@ using namespace std; * use a block storage scheme rather than a continuous storage * scheme). */ +template class CSysVector { private: @@ -71,7 +76,7 @@ class CSysVector { unsigned short nVar; /*!< \brief number of elements in a block */ unsigned long nBlk; /*!< \brief number of blocks (or number of blocks on this processor) */ unsigned long nBlkDomain; /*!< \brief number of blocks (or number of blocks on this processor without Ghost cells) */ - su2double* vec_val; /*!< \brief storage for the element values */ + ScalarType* vec_val; /*!< \brief storage for the element values */ public: @@ -85,7 +90,7 @@ class CSysVector { * \param[in] size - number of elements locally * \param[in] val - default value for elements */ - CSysVector(const unsigned long & size, const su2double & val = 0.0); + CSysVector(const unsigned long & size, const ScalarType & val = 0.0); /*! * \brief constructor of the class. @@ -94,7 +99,7 @@ class CSysVector { * \param[in] numVar - number of variables in each block * \param[in] val - default value for elements */ - CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const su2double & val = 0.0); + CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const ScalarType & val = 0.0); /*! * \brief copy constructor of the class. @@ -112,7 +117,7 @@ class CSysVector { * \param[in] size - number of elements locally * \param[in] u_array - vector stored as array being copied */ - explicit CSysVector(const unsigned long & size, const su2double* u_array); + explicit CSysVector(const unsigned long & size, const ScalarType* u_array); /*! * \brief constructor from array @@ -122,7 +127,7 @@ class CSysVector { * \param[in] u_array - vector stored as array being copied */ explicit CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, - const su2double* u_array); + const ScalarType* u_array); /*! * \brief class destructor @@ -136,13 +141,18 @@ class CSysVector { * \param[in] numVar - number of variables in each block * \param[in] val - default value for elements */ - void Initialize(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const su2double & val = 0.0); + void Initialize(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const ScalarType & val = 0.0); /*! * \brief return the number of local elements in the CSysVector */ unsigned long GetLocSize() const; + /*! + * \brief return the number of local elements in the CSysVector without ghost cells + */ + unsigned long GetNElmDomain() const; + /*! * \brief return the size of the CSysVector (over all processors) */ @@ -168,14 +178,14 @@ class CSysVector { * \param[in] a - scalar factor for x * \param[in] x - CSysVector that is being scaled */ - void Equals_AX(const su2double & a, CSysVector & x); + void Equals_AX(const ScalarType & a, CSysVector & x); /*! * \brief adds a scaled CSysVector to calling CSysVector * \param[in] a - scalar factor for x * \param[in] x - CSysVector that is being scaled */ - void Plus_AX(const su2double & a, CSysVector & x); + void Plus_AX(const ScalarType & a, CSysVector & x); /*! * \brief general linear combination of two CSysVectors @@ -184,7 +194,7 @@ class CSysVector { * \param[in] b - scalar factor for y * \param[in] y - second CSysVector in linear combination */ - void Equals_AX_Plus_BY(const su2double & a, CSysVector & x, const su2double & b, CSysVector & y); + void Equals_AX_Plus_BY(const ScalarType & a, CSysVector & x, const ScalarType & b, CSysVector & y); /*! * \brief assignment operator with deep copy @@ -196,7 +206,7 @@ class CSysVector { * \brief CSysVector=su2double assignment operator * \param[in] val - value assigned to each element of CSysVector */ - CSysVector & operator=(const su2double & val); + CSysVector & operator=(const ScalarType & val); /*! * \brief addition operator @@ -226,71 +236,64 @@ class CSysVector { * \brief vector * scalar multiplication operator * \param[in] val - value to multiply *this by */ - CSysVector operator*(const su2double & val) const; - - /*! - * \brief scalar * vector multiplication operator - * \param[in] val - scalar value to multiply by - * \param[in] u - CSysVector having its elements scaled - */ - friend CSysVector operator*(const su2double & val, const CSysVector & u); + CSysVector operator*(const ScalarType & val) const; /*! * \brief compound scalar multiplication-assignment operator * \param[in] val - value to multiply calling object by */ - CSysVector & operator*=(const su2double & val); + CSysVector & operator*=(const ScalarType & val); /*! * \brief vector-scalar division operator (no scalar/vector operator) * \param[in] val - value to divide elements of *this by */ - CSysVector operator/(const su2double & val) const; + CSysVector operator/(const ScalarType & val) const; /*! * \brief compound scalar division-assignment operator * \param[in] val - value to divide elements of calling object by */ - CSysVector & operator/=(const su2double & val); + CSysVector & operator/=(const ScalarType & val); /*! * \brief indexing operator with assignment permitted * \param[in] i = local index to access */ - su2double & operator[](const unsigned long & i); + ScalarType & operator[](const unsigned long & i); /*! * \brief indexing operator with assignment not permitted * \param[in] i = local index to access */ - const su2double & operator[](const unsigned long & i) const; + const ScalarType & operator[](const unsigned long & i) const; /*! * \brief the L2 norm of the CSysVector * \result the L2 norm */ - su2double norm() const; + ScalarType norm() const; /*! * \brief copies the contents of the calling CSysVector into an array * \param[out] u_array - array into which information is being copied * \pre u_array must be allocated and have the same size as CSysVector */ - void CopyToArray(su2double* u_array); + void CopyToArray(ScalarType* u_array); /*! * \brief Subtract val_residual to the residual. * \param[in] val_ipoint - index of the point where subtract the residual. * \param[in] val_residual - Value to subtract to the residual. */ - void SubtractBlock(unsigned long val_ipoint, su2double *val_residual); + void SubtractBlock(unsigned long val_ipoint, ScalarType *val_residual); /*! * \brief Add val_residual to the residual. * \param[in] val_ipoint - index of the point where add the residual. * \param[in] val_residual - Value to add to the residual. */ - void AddBlock(unsigned long val_ipoint, su2double *val_residual); + void AddBlock(unsigned long val_ipoint, ScalarType *val_residual); /*! * \brief Set val_residual to the residual. @@ -298,14 +301,14 @@ class CSysVector { * \param[in] val_var - inde of the residual to be set. * \param[in] val_residual - Value to set to the residual. */ - void SetBlock(unsigned long val_ipoint, unsigned short val_var, su2double val_residual); + void SetBlock(unsigned long val_ipoint, unsigned short val_var, ScalarType val_residual); /*! * \brief Set val_residual to the residual. * \param[in] val_ipoint - index of the point where set the residual. * \param[in] val_residual - Value to set to the residual. */ - void SetBlock(unsigned long val_ipoint, su2double *val_residual); + void SetBlock(unsigned long val_ipoint, ScalarType *val_residual); /*! * \brief Set the residual to zero. @@ -325,7 +328,7 @@ class CSysVector { * \param[in] val_ipoint - index of the point where set the residual. * \return Pointer to the residual. */ - su2double *GetBlock(unsigned long val_ipoint); + ScalarType *GetBlock(unsigned long val_ipoint); /*! * \brief Get the value of the residual. @@ -333,18 +336,31 @@ class CSysVector { * \param[in] val_var - inde of the residual to be set. * \return Value of the residual. */ - su2double GetBlock(unsigned long val_ipoint, unsigned short val_var); - + ScalarType GetBlock(unsigned long val_ipoint, unsigned short val_var); /*! * \brief dot-product between two CSysVectors * \param[in] u - first CSysVector in dot product * \param[in] v - second CSysVector in dot product */ - friend su2double dotProd(const CSysVector & u, const CSysVector & v); + friend ScalarType dotProd(const CSysVector & u, const CSysVector & v); + /*! + * \brief Set our values (resizing if required) by copying from other, the derivative information is lost. + * \param[in] other - source CSysVector + */ + template + void PassiveCopy(const CSysVector& other); }; +/*! + * \brief scalar * vector multiplication operator + * \param[in] val - scalar value to multiply by + * \param[in] u - CSysVector having its elements scaled + */ +template +CSysVector operator*(const ScalarType & val, const CSysVector & u); + /*! * \class CMatrixVectorProduct * \brief abstract base class for defining matrix-vector products @@ -361,13 +377,15 @@ class CSysVector { * handle the different types of matrix-vector products and still be * passed to a single implementation of the Krylov solvers. */ +template class CMatrixVectorProduct { public: virtual ~CMatrixVectorProduct() = 0; ///< class destructor - virtual void operator()(const CSysVector & u, CSysVector & v) + virtual void operator()(const CSysVector & u, CSysVector & v) const = 0; ///< matrix-vector product operation }; -inline CMatrixVectorProduct::~CMatrixVectorProduct() {} +template +inline CMatrixVectorProduct::~CMatrixVectorProduct() {} /*! * \class CPreconditioner @@ -377,12 +395,14 @@ inline CMatrixVectorProduct::~CMatrixVectorProduct() {} * See the remarks regarding the CMatrixVectorProduct class. The same * idea applies here to the preconditioning operation. */ +template class CPreconditioner { public: virtual ~CPreconditioner() = 0; ///< class destructor - virtual void operator()(const CSysVector & u, CSysVector & v) + virtual void operator()(const CSysVector & u, CSysVector & v) const = 0; ///< preconditioning operation }; -inline CPreconditioner::~CPreconditioner() {} +template +inline CPreconditioner::~CPreconditioner() {} #include "vector_structure.inl" diff --git a/Common/include/vector_structure.inl b/Common/include/vector_structure.inl index 074be5da3ddb..6d724333ee16 100644 --- a/Common/include/vector_structure.inl +++ b/Common/include/vector_structure.inl @@ -37,14 +37,20 @@ #pragma once -inline void CSysVector::SetValZero(void) { +template +inline void CSysVector::SetValZero(void) { for (unsigned long i = 0; i < nElm; i++) vec_val[i] = 0.0; } -inline unsigned long CSysVector::GetLocSize() const { return nElm; } +template +inline unsigned long CSysVector::GetLocSize() const { return nElm; } -inline unsigned long CSysVector::GetSize() const { +template +unsigned long CSysVector::GetNElmDomain() const { return nElmDomain; } + +template +inline unsigned long CSysVector::GetSize() const { #ifdef HAVE_MPI return nElmGlobal; #else @@ -52,12 +58,51 @@ inline unsigned long CSysVector::GetSize() const { #endif } -inline unsigned short CSysVector::GetNVar() const { return nVar; } +template +inline unsigned short CSysVector::GetNVar() const { return nVar; } + +template +inline unsigned long CSysVector::GetNBlk() const { return nBlk; } + +template +inline unsigned long CSysVector::GetNBlkDomain() const { return nBlkDomain; } + +template +inline ScalarType & CSysVector::operator[](const unsigned long & i) { return vec_val[i]; } + +template +inline const ScalarType & CSysVector::operator[](const unsigned long & i) const { return vec_val[i]; } + +template +template +void CSysVector::PassiveCopy(const CSysVector& other) { + + /*--- This is a method and not the overload of an operator to make sure who + calls it knows the consequence to the derivative information (lost) ---*/ -inline unsigned long CSysVector::GetNBlk() const { return nBlk; } + /*--- check if self-assignment, otherwise perform deep copy ---*/ + if ((const void*)this == (const void*)&other) return; -inline unsigned long CSysVector::GetNBlkDomain() const { return nBlkDomain; } + /*--- determine if (re-)allocation is needed ---*/ + if (nElm != other.GetLocSize() && vec_val != NULL) { + delete [] vec_val; + vec_val = NULL; + } -inline su2double & CSysVector::operator[](const unsigned long & i) { return vec_val[i]; } + /*--- copy ---*/ + nElm = other.GetLocSize(); + nElmDomain = other.GetNElmDomain(); + nBlk = other.GetNBlk(); + nBlkDomain = other.GetNBlkDomain(); + nVar = other.GetNVar(); -inline const su2double & CSysVector::operator[](const unsigned long & i) const { return vec_val[i]; } + if (vec_val == NULL) + vec_val = new ScalarType[nElm]; + + for (unsigned long i = 0; i < nElm; i++) + vec_val[i] = SU2_TYPE::GetValue(other[i]); + +#ifdef HAVE_MPI + nElmGlobal = other.GetSize(); +#endif +} diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index 073dee7ecc57..35ac618d0bd7 100644 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -186,8 +186,8 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co if (Derivative) { SetBoundaryDerivatives(geometry, config); } - CMatrixVectorProduct* mat_vec = NULL; - CPreconditioner* precond = NULL; + CMatrixVectorProduct* mat_vec = NULL; + CPreconditioner* precond = NULL; /*--- Communicate any prescribed boundary displacements via MPI, so that all nodes have the same solution and r.h.s. entries @@ -209,20 +209,20 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co if (config->GetKind_Deform_Linear_Solver_Prec() == LU_SGS) { if ((rank == MASTER_NODE) && Screen_Output) cout << "\n# LU_SGS preconditioner." << endl; - mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); - precond = new CLU_SGSPreconditioner(StiffMatrix, geometry, config); + mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); + precond = new CLU_SGSPreconditioner(StiffMatrix, geometry, config); } if (config->GetKind_Deform_Linear_Solver_Prec() == ILU) { if ((rank == MASTER_NODE) && Screen_Output) cout << "\n# ILU preconditioner." << endl; StiffMatrix.BuildILUPreconditioner(); - mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); - precond = new CILUPreconditioner(StiffMatrix, geometry, config); + mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); + precond = new CILUPreconditioner(StiffMatrix, geometry, config); } if (config->GetKind_Deform_Linear_Solver_Prec() == JACOBI) { if ((rank == MASTER_NODE) && Screen_Output) cout << "\n# Jacobi preconditioner." << endl; StiffMatrix.BuildJacobiPreconditioner(); - mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); - precond = new CJacobiPreconditioner(StiffMatrix, geometry, config); + mat_vec = new CSysMatrixVectorProduct(StiffMatrix, geometry, config); + precond = new CJacobiPreconditioner(StiffMatrix, geometry, config); } } else if (Derivative && (config->GetKind_SU2() == SU2_DOT)) { @@ -233,14 +233,14 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co (config->GetKind_Deform_Linear_Solver_Prec() == LU_SGS)) { if ((rank == MASTER_NODE) && Screen_Output) cout << "\n# ILU preconditioner." << endl; StiffMatrix.BuildILUPreconditioner(true); - mat_vec = new CSysMatrixVectorProductTransposed(StiffMatrix, geometry, config); - precond = new CILUPreconditioner(StiffMatrix, geometry, config); + mat_vec = new CSysMatrixVectorProductTransposed(StiffMatrix, geometry, config); + precond = new CILUPreconditioner(StiffMatrix, geometry, config); } if (config->GetKind_Deform_Linear_Solver_Prec() == JACOBI) { if ((rank == MASTER_NODE) && Screen_Output) cout << "\n# Jacobi preconditioner." << endl; StiffMatrix.BuildJacobiPreconditioner(true); - mat_vec = new CSysMatrixVectorProductTransposed(StiffMatrix, geometry, config); - precond = new CJacobiPreconditioner(StiffMatrix, geometry, config); + mat_vec = new CSysMatrixVectorProductTransposed(StiffMatrix, geometry, config); + precond = new CJacobiPreconditioner(StiffMatrix, geometry, config); } } diff --git a/Common/src/linear_solvers_structure.cpp b/Common/src/linear_solvers_structure.cpp index 32ee01388549..166cfec38c9e 100644 --- a/Common/src/linear_solvers_structure.cpp +++ b/Common/src/linear_solvers_structure.cpp @@ -38,32 +38,37 @@ #include "../include/linear_solvers_structure.hpp" #include "../include/linear_solvers_structure_b.hpp" -CSysSolve::CSysSolve(const bool mesh_deform_mode) : cg_ready(false), bcg_ready(false), gmres_ready(false) { +template +CSysSolve::CSysSolve(const bool mesh_deform_mode) : cg_ready(false), bcg_ready(false), gmres_ready(false) { mesh_deform = mesh_deform_mode; + LinSysRes_ptr = NULL; + LinSysSol_ptr = NULL; } -void CSysSolve::ApplyGivens(const su2double & s, const su2double & c, su2double & h1, su2double & h2) { +template +void CSysSolve::ApplyGivens(const ScalarType & s, const ScalarType & c, ScalarType & h1, ScalarType & h2) { - su2double temp = c*h1 + s*h2; + ScalarType temp = c*h1 + s*h2; h2 = c*h2 - s*h1; h1 = temp; } -void CSysSolve::GenerateGivens(su2double & dx, su2double & dy, su2double & s, su2double & c) { +template +void CSysSolve::GenerateGivens(ScalarType & dx, ScalarType & dy, ScalarType & s, ScalarType & c) { if ( (dx == 0.0) && (dy == 0.0) ) { c = 1.0; s = 0.0; } else if ( fabs(dy) > fabs(dx) ) { - su2double tmp = dx/dy; + ScalarType tmp = dx/dy; dx = sqrt(1.0 + tmp*tmp); s = Sign(1.0/dx, dy); c = tmp*s; } else if ( fabs(dy) <= fabs(dx) ) { - su2double tmp = dy/dx; + ScalarType tmp = dy/dx; dy = sqrt(1.0 + tmp*tmp); c = Sign(1.0/dy, dx); s = tmp*c; @@ -79,8 +84,9 @@ void CSysSolve::GenerateGivens(su2double & dx, su2double & dy, su2double & s, su dy = 0.0; } -void CSysSolve::SolveReduced(const int & n, const vector > & Hsbg, - const vector & rhs, vector & x) { +template +void CSysSolve::SolveReduced(const int & n, const vector > & Hsbg, + const vector & rhs, vector & x) { // initialize... for (int i = 0; i < n; i++) x[i] = rhs[i]; @@ -93,19 +99,20 @@ void CSysSolve::SolveReduced(const int & n, const vector > & H } } -void CSysSolve::ModGramSchmidt(int i, vector > & Hsbg, vector & w) { +template +void CSysSolve::ModGramSchmidt(int i, vector > & Hsbg, vector > & w) { bool Convergence = true; /*--- Parameter for reorthonormalization ---*/ - static const su2double reorth = 0.98; + static const ScalarType reorth = 0.98; /*--- Get the norm of the vector being orthogonalized, and find the threshold for re-orthogonalization ---*/ - su2double nrm = dotProd(w[i+1], w[i+1]); - su2double thr = nrm*reorth; + ScalarType nrm = dotProd(w[i+1], w[i+1]); + ScalarType thr = nrm*reorth; /*--- The norm of w[i+1] < 0.0 or w[i+1] = NaN ---*/ @@ -152,7 +159,7 @@ void CSysSolve::ModGramSchmidt(int i, vector > & Hsbg, vector< /*--- Begin main Gram-Schmidt loop ---*/ for (int k = 0; k < i+1; k++) { - su2double prod = dotProd(w[i+1], w[k]); + ScalarType prod = dotProd(w[i+1], w[k]); Hsbg[k][i] = prod; w[i+1].Plus_AX(-prod, w[k]); @@ -182,7 +189,8 @@ void CSysSolve::ModGramSchmidt(int i, vector > & Hsbg, vector< } -void CSysSolve::WriteHeader(const string & solver, const su2double & restol, const su2double & resinit) { +template +void CSysSolve::WriteHeader(const string & solver, const ScalarType & restol, const ScalarType & resinit) { cout << "\n# " << solver << " residual history" << endl; cout << "# Residual tolerance target = " << restol << endl; @@ -190,17 +198,20 @@ void CSysSolve::WriteHeader(const string & solver, const su2double & restol, con } -void CSysSolve::WriteHistory(const int & iter, const su2double & res, const su2double & resinit) { +template +void CSysSolve::WriteHistory(const int & iter, const ScalarType & res, const ScalarType & resinit) { cout << " " << iter << " " << res/resinit << endl; } -unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, unsigned long m, su2double *residual, bool monitoring, CConfig *config) { +template +unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, + CMatrixVectorProduct & mat_vec, CPreconditioner & precond, + ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CConfig *config) { int rank = SU2_MPI::GetRank(); - su2double norm_r = 0.0, norm0 = 0.0; + ScalarType norm_r = 0.0, norm0 = 0.0; int i = 0; /*--- Check the subspace size ---*/ @@ -248,7 +259,7 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, CMat } - su2double alpha, beta, r_dot_z; + ScalarType alpha, beta, r_dot_z; precond(r, z); p = z; @@ -309,7 +320,7 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, CMat mat_vec(x, A_x); r = b; r -= A_x; - su2double true_res = r.norm(); + ScalarType true_res = r.norm(); if (fabs(true_res - norm_r) > tol*10.0) { if (rank == MASTER_NODE) { @@ -327,8 +338,10 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, CMat } -unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, unsigned long m, su2double *residual, bool monitoring, CConfig *config) { +template +unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, + CMatrixVectorProduct & mat_vec, CPreconditioner & precond, + ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CConfig *config) { int rank = SU2_MPI::GetRank(); @@ -360,15 +373,15 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, /*--- Define various arrays ---*/ - vector g(m+1, 0.0); - vector sn(m+1, 0.0); - vector cs(m+1, 0.0); - vector y(m, 0.0); - vector > H(m+1, vector(m, 0.0)); + vector g(m+1, 0.0); + vector sn(m+1, 0.0); + vector cs(m+1, 0.0); + vector y(m, 0.0); + vector > H(m+1, vector(m, 0.0)); /*--- Calculate the norm of the rhs vector ---*/ - su2double norm0 = b.norm(); + ScalarType norm0 = b.norm(); /*--- Calculate the initial residual (actually the negative residual) and compute its norm ---*/ @@ -376,7 +389,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, mat_vec(x, W[0]); W[0] -= b; - su2double beta = W[0].norm(); + ScalarType beta = W[0].norm(); if ( (beta < tol*norm0) || (beta < eps) ) { @@ -453,17 +466,18 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, x.Plus_AX(y[k], Z[k]); } - if ((monitoring) && (rank == MASTER_NODE)) { - cout << "# FGMRES final (true) residual:" << endl; - cout << "# Iteration = " << i << ": |res|/|res0| = " << beta/norm0 << ".\n" << endl; - } - /*--- Recalculate final (neg.) residual (this should be optional) ---*/ - if (monitoring) { + if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) { + + if (rank == MASTER_NODE) { + cout << "# FGMRES final (true) residual:" << endl; + cout << "# Iteration = " << i << ": |res|/|res0| = " << beta/norm0 << ".\n" << endl; + } + mat_vec(x, W[0]); W[0] -= b; - su2double res = W[0].norm(); + ScalarType res = W[0].norm(); if (fabs(res - beta) > tol*10) { if (rank == MASTER_NODE) { @@ -473,6 +487,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, cout << "# res - beta = " << res - beta << endl << endl; } } + } (*residual) = beta; @@ -480,11 +495,13 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector & b, CSysVector & x, } -unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, - CPreconditioner & precond, su2double tol, unsigned long m, su2double *residual, bool monitoring, CConfig *config) { +template +unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, + CMatrixVectorProduct & mat_vec, CPreconditioner & precond, + ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CConfig *config) { int rank = SU2_MPI::GetRank(); - su2double norm_r = 0.0, norm0 = 0.0; + ScalarType norm_r = 0.0, norm0 = 0.0; int i = 0; /*--- Check the subspace size ---*/ @@ -536,8 +553,8 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, /*--- Initialization ---*/ - su2double alpha = 1.0, beta = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0; - p = su2double(0.0); v = su2double(0.0); r_0 = r; + ScalarType alpha = 1.0, beta = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0; + p = ScalarType(0.0); v = ScalarType(0.0); r_0 = r; /*--- Loop over all search directions ---*/ @@ -557,7 +574,7 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, /*--- p_{i} = r_{i-1} + beta * p_{i-1} - beta * omega * v_{i-1} ---*/ - su2double beta_omega = -beta*omega; + ScalarType beta_omega = -beta*omega; p.Equals_AX_Plus_BY(beta, p, beta_omega, v); p.Plus_AX(1.0, r); @@ -568,7 +585,7 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, /*--- Calculate step-length alpha ---*/ - su2double r_0_v = dotProd(r_0, v); + ScalarType r_0_v = dotProd(r_0, v); alpha = rho / r_0_v; /*--- Update solution and residual: ---*/ @@ -619,7 +636,7 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, mat_vec(x, A_x); r = b; r -= A_x; - su2double true_res = r.norm(); + ScalarType true_res = r.norm(); if ((fabs(true_res - norm_r) > tol*10.0) && (rank == MASTER_NODE)) { cout << "# WARNING in CSysSolve::BCGSTAB_LinSolver(): " << endl; @@ -634,11 +651,67 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVector & b, CSysVector & x, return (unsigned long) i; } -unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) { +template<> +void CSysSolve::HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol) { + + /*--- When the type is the same the temporaties are not required ---*/ + /*--- Set the pointers ---*/ + LinSysRes_ptr = &LinSysRes; + LinSysSol_ptr = &LinSysSol; +} + +template<> +void CSysSolve::HandleTemporariesOut(CSysVector & LinSysSol) { + + /*--- When the type is the same the temporaties are not required ---*/ + /*--- Reset the pointers ---*/ + LinSysRes_ptr = NULL; + LinSysSol_ptr = NULL; +} + +#ifdef CODI_REVERSE_TYPE +template<> +void CSysSolve::HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol) { + + /*--- When the type is different we need to copy data to the temporaries ---*/ + /*--- Copy data, the solution is also copied because it serves as initial conditions ---*/ + LinSysRes_tmp.PassiveCopy(LinSysRes); + LinSysSol_tmp.PassiveCopy(LinSysSol); + + /*--- Set the pointers ---*/ + LinSysRes_ptr = &LinSysRes_tmp; + LinSysSol_ptr = &LinSysSol_tmp; +} + +template<> +void CSysSolve::HandleTemporariesOut(CSysVector & LinSysSol) { + + /*--- When the type is different we need to copy data from the temporaries ---*/ + /*--- Copy data, only the solution needs to be copied ---*/ + LinSysSol.PassiveCopy(LinSysSol_tmp); + + /*--- Reset the pointers ---*/ + LinSysRes_ptr = NULL; + LinSysSol_ptr = NULL; +} +#endif + +template +unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, + CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) { + /*--- + A word about the templated types. It is assumed that the residual and solution vectors are always of su2doubles, + meaning that they are active in the discrete adjoint. The same assumption is made in SetExternalSolve. + When the Jacobian is passive (and therefore not compatible with the vectors) we go through the "HandleTemporaries" + mechanisms. Note that CG, BCGSTAB, and FGMRES, all expect the vector to be compatible with the Product and + Preconditioner (and therefore with the Matrix). Likewise for Solve_b (which is used by CSysSolve_b). + There are no provisions here for active Matrix and passive Vectors as that makes no sense since we only handle the + derivatives of the residual in CSysSolve_b. + ---*/ unsigned short KindSolver, KindPrecond; unsigned long MaxIter, RestartIter, IterLinSol = 0; - su2double SolverTol, Norm0 = 0.0; + ScalarType SolverTol, Norm0 = 0.0; bool ScreenOutput; /*--- Normal mode ---*/ @@ -649,7 +722,7 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS KindPrecond = config->GetKind_Linear_Solver_Prec(); MaxIter = config->GetLinear_Solver_Iter(); RestartIter = config->GetLinear_Solver_Restart_Frequency(); - SolverTol = config->GetLinear_Solver_Error(); + SolverTol = SU2_TYPE::GetValue(config->GetLinear_Solver_Error()); ScreenOutput = false; } @@ -661,11 +734,11 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS KindPrecond = config->GetKind_Deform_Linear_Solver_Prec(); MaxIter = config->GetDeform_Linear_Solver_Iter(); RestartIter = config->GetLinear_Solver_Restart_Frequency(); - SolverTol = config->GetDeform_Linear_Solver_Error(); + SolverTol = SU2_TYPE::GetValue(config->GetDeform_Linear_Solver_Error()); ScreenOutput = config->GetDeform_Output(); } - CMatrixVectorProduct *mat_vec = NULL; + CMatrixVectorProduct *mat_vec = NULL; bool TapeActive = NO; @@ -685,52 +758,54 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS } /*--- Solve the linear system using a Krylov subspace method ---*/ - + + HandleTemporariesIn(LinSysRes, LinSysSol); + if (KindSolver == BCGSTAB || KindSolver == CONJUGATE_GRADIENT || KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES ) { - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - CPreconditioner* precond = NULL; + mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); + CPreconditioner* precond = NULL; switch (KindPrecond) { case JACOBI: Jacobian.BuildJacobiPreconditioner(); - precond = new CJacobiPreconditioner(Jacobian, geometry, config); + precond = new CJacobiPreconditioner(Jacobian, geometry, config); break; case ILU: Jacobian.BuildILUPreconditioner(); - precond = new CILUPreconditioner(Jacobian, geometry, config); + precond = new CILUPreconditioner(Jacobian, geometry, config); break; case LU_SGS: - precond = new CLU_SGSPreconditioner(Jacobian, geometry, config); + precond = new CLU_SGSPreconditioner(Jacobian, geometry, config); break; case LINELET: Jacobian.BuildJacobiPreconditioner(); - precond = new CLineletPreconditioner(Jacobian, geometry, config); + precond = new CLineletPreconditioner(Jacobian, geometry, config); break; default: Jacobian.BuildJacobiPreconditioner(); - precond = new CJacobiPreconditioner(Jacobian, geometry, config); + precond = new CJacobiPreconditioner(Jacobian, geometry, config); break; } switch (KindSolver) { case BCGSTAB: - IterLinSol = BCGSTAB_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); break; case FGMRES: - IterLinSol = FGMRES_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); break; case CONJUGATE_GRADIENT: - IterLinSol = CG_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); break; case RESTARTED_FGMRES: IterLinSol = 0; - Norm0 = LinSysRes.norm(); + Norm0 = LinSysRes_ptr->norm(); while (IterLinSol < MaxIter) { /*--- Enforce a hard limit on total number of iterations ---*/ unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol); - IterLinSol += FGMRES_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol, IterLimit, &Residual, ScreenOutput, config); + IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, IterLimit, &Residual, ScreenOutput, config); if ( Residual < SolverTol*Norm0 ) break; } break; @@ -748,30 +823,31 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS else { switch (KindSolver) { case SMOOTHER_LUSGS: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - IterLinSol = Jacobian.LU_SGS_Smoother(LinSysRes, LinSysSol, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); + mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); + IterLinSol = Jacobian.LU_SGS_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); delete mat_vec; break; case SMOOTHER_JACOBI: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); + mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); Jacobian.BuildJacobiPreconditioner(); - IterLinSol = Jacobian.Jacobi_Smoother(LinSysRes, LinSysSol, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); + IterLinSol = Jacobian.Jacobi_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); delete mat_vec; break; case SMOOTHER_ILU: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); + mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); Jacobian.BuildILUPreconditioner(); - IterLinSol = Jacobian.ILU_Smoother(LinSysRes, LinSysSol, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); + IterLinSol = Jacobian.ILU_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); delete mat_vec; break; case SMOOTHER_LINELET: Jacobian.BuildJacobiPreconditioner(); - Jacobian.ComputeLineletPreconditioner(LinSysRes, LinSysSol, geometry, config); + Jacobian.ComputeLineletPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config); IterLinSol = 1; break; } } + HandleTemporariesOut(LinSysSol); if(TapeActive) { @@ -784,7 +860,7 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS AD::StartRecording(); AD::SetExtFuncOut(&LinSysSol[0], (int)LinSysSol.GetLocSize()); - + #ifdef CODI_REVERSE_TYPE AD::FuncHelper->addUserData(&LinSysRes); AD::FuncHelper->addUserData(&LinSysSol); @@ -792,38 +868,37 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS AD::FuncHelper->addUserData(geometry); AD::FuncHelper->addUserData(config); AD::FuncHelper->addUserData(this); - AD::FuncHelper->addToTape(CSysSolve_b::Solve_b); + AD::FuncHelper->addToTape(CSysSolve_b::Solve_b); #endif - /*--- Build preconditioner for the transposed Jacobian ---*/ + /*--- Build preconditioner for the transposed Jacobian ---*/ - switch(KindPrecond) { - case ILU: - Jacobian.BuildILUPreconditioner(RequiresTranspose); - break; - case JACOBI: - Jacobian.BuildJacobiPreconditioner(RequiresTranspose); - break; - default: - SU2_MPI::Error("The specified preconditioner is not yet implemented for the discrete adjoint method.", CURRENT_FUNCTION); - break; - } - - AD::EndExtFunc(); + switch(KindPrecond) { + case ILU: + Jacobian.BuildILUPreconditioner(RequiresTranspose); + break; + case JACOBI: + Jacobian.BuildJacobiPreconditioner(RequiresTranspose); + break; + default: + SU2_MPI::Error("The specified preconditioner is not yet implemented for the discrete adjoint method.", CURRENT_FUNCTION); + break; + } + AD::EndExtFunc(); } return IterLinSol; - } - -unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) { +template +unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, + CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) { #ifdef CODI_REVERSE_TYPE unsigned short KindSolver, KindPrecond; unsigned long MaxIter, RestartIter, IterLinSol = 0; - su2double SolverTol, Norm0 = 0.0; + ScalarType SolverTol, Norm0 = 0.0; bool ScreenOutput; /*--- Normal mode ---*/ @@ -834,7 +909,7 @@ unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, KindPrecond = config->GetKind_DiscAdj_Linear_Prec(); MaxIter = config->GetLinear_Solver_Iter(); RestartIter = config->GetLinear_Solver_Restart_Frequency(); - SolverTol = config->GetLinear_Solver_Error(); + SolverTol = SU2_TYPE::GetValue(config->GetLinear_Solver_Error()); ScreenOutput = false; } @@ -846,44 +921,46 @@ unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, KindPrecond = config->GetKind_Deform_Linear_Solver_Prec(); MaxIter = config->GetDeform_Linear_Solver_Iter(); RestartIter = config->GetLinear_Solver_Restart_Frequency(); - SolverTol = config->GetDeform_Linear_Solver_Error(); + SolverTol = SU2_TYPE::GetValue(config->GetDeform_Linear_Solver_Error()); ScreenOutput = config->GetDeform_Output(); } /*--- Set up preconditioner and matrix-vector product ---*/ - CPreconditioner* precond = NULL; + CPreconditioner* precond = NULL; switch(KindPrecond) { case ILU: - precond = new CILUPreconditioner(Jacobian, geometry, config); + precond = new CILUPreconditioner(Jacobian, geometry, config); break; case JACOBI: - precond = new CJacobiPreconditioner(Jacobian, geometry, config); + precond = new CJacobiPreconditioner(Jacobian, geometry, config); break; } - CMatrixVectorProduct* mat_vec = new CSysMatrixVectorProductTransposed(Jacobian, geometry, config); + CMatrixVectorProduct* mat_vec = new CSysMatrixVectorProductTransposed(Jacobian, geometry, config); /*--- Solve the system ---*/ + HandleTemporariesIn(LinSysRes, LinSysSol); + switch(KindSolver) { case FGMRES: - IterLinSol = FGMRES_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol , MaxIter, &Residual, ScreenOutput, config); + IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol , MaxIter, &Residual, ScreenOutput, config); break; case BCGSTAB: - IterLinSol = BCGSTAB_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol , MaxIter, &Residual, ScreenOutput, config); + IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol , MaxIter, &Residual, ScreenOutput, config); break; case CONJUGATE_GRADIENT: - IterLinSol = CG_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); break; case RESTARTED_FGMRES: IterLinSol = 0; - Norm0 = LinSysRes.norm(); + Norm0 = LinSysRes_ptr->norm(); while (IterLinSol < MaxIter) { /*--- Enforce a hard limit on total number of iterations ---*/ unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol); - IterLinSol += FGMRES_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, SolverTol , IterLimit, &Residual, ScreenOutput, config); + IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol , IterLimit, &Residual, ScreenOutput, config); if ( Residual < SolverTol*Norm0 ) break; } break; @@ -892,6 +969,8 @@ unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, break; } + HandleTemporariesOut(LinSysSol); + delete mat_vec; delete precond; @@ -901,3 +980,9 @@ unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, #endif } +/*--- Explicit instantiations ---*/ +template class CSysSolve; + +#ifdef CODI_REVERSE_TYPE +template class CSysSolve; +#endif diff --git a/Common/src/linear_solvers_structure_b.cpp b/Common/src/linear_solvers_structure_b.cpp index be986badefa2..a3092ed01ff8 100644 --- a/Common/src/linear_solvers_structure_b.cpp +++ b/Common/src/linear_solvers_structure_b.cpp @@ -41,26 +41,28 @@ #include "../include/matrix_structure.hpp" #ifdef CODI_REVERSE_TYPE -void CSysSolve_b::Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m, const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n, codi::DataStore* d){ +template +void CSysSolve_b::Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m, + const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n, + codi::DataStore* d) { - CSysVector* LinSysRes_b = NULL; + CSysVector* LinSysRes_b = NULL; d->getData(LinSysRes_b); - CSysVector* LinSysSol_b = NULL; + CSysVector* LinSysSol_b = NULL; d->getData(LinSysSol_b); - CSysMatrix* Jacobian = NULL; + CSysMatrix* Jacobian = NULL; d->getData(Jacobian); - + CGeometry* geometry = NULL; d->getData(geometry); - + CConfig* config = NULL; d->getData(config); - CSysSolve* solver; + CSysSolve* solver = NULL; d->getData(solver); - /*--- Initialize the right-hand side with the gradient of the solution of the primal linear system ---*/ @@ -70,10 +72,14 @@ void CSysSolve_b::Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::R } solver->Solve_b(*Jacobian, *LinSysRes_b, *LinSysSol_b, geometry, config); - + for (unsigned long i = 0; i < n; i ++) { x_b[i] = SU2_TYPE::GetValue(LinSysSol_b->operator [](i)); } } + +template class CSysSolve_b; +template class CSysSolve_b; + #endif diff --git a/Common/src/matrix_structure.cpp b/Common/src/matrix_structure.cpp index 82efa98295a6..840da34e44ef 100644 --- a/Common/src/matrix_structure.cpp +++ b/Common/src/matrix_structure.cpp @@ -37,7 +37,8 @@ #include "../include/matrix_structure.hpp" -CSysMatrix::CSysMatrix(void) { +template +CSysMatrix::CSysMatrix(void) { size = SU2_MPI::GetSize(); rank = SU2_MPI::GetRank(); @@ -85,7 +86,8 @@ CSysMatrix::CSysMatrix(void) { } -CSysMatrix::~CSysMatrix(void) { +template +CSysMatrix::~CSysMatrix(void) { unsigned long iElem; @@ -140,7 +142,8 @@ CSysMatrix::~CSysMatrix(void) { } -void CSysMatrix::Initialize(unsigned long nPoint, unsigned long nPointDomain, +template +void CSysMatrix::Initialize(unsigned long nPoint, unsigned long nPointDomain, unsigned short nVar, unsigned short nEqn, bool EdgeConnect, CGeometry *geometry, CConfig *config) { @@ -293,7 +296,7 @@ void CSysMatrix::Initialize(unsigned long nPoint, unsigned long nPointDomain, } - ILU_matrix = new su2double [nnz_ilu*nVar*nEqn]; + ILU_matrix = new ScalarType [nnz_ilu*nVar*nEqn]; for (iVar = 0; iVar < nnz_ilu*nVar*nEqn; iVar++) ILU_matrix[iVar] = 0.0; delete [] nNeigh_ilu; @@ -302,7 +305,8 @@ void CSysMatrix::Initialize(unsigned long nPoint, unsigned long nPointDomain, } -void CSysMatrix::SetNeighbours(CGeometry *geometry, unsigned long iPoint, unsigned short deep_level, unsigned short fill_level, +template +void CSysMatrix::SetNeighbours(CGeometry *geometry, unsigned long iPoint, unsigned short deep_level, unsigned short fill_level, bool EdgeConnect, vector & vneighs) { unsigned long Point, iElem, Elem; unsigned short iNode; @@ -329,7 +333,8 @@ void CSysMatrix::SetNeighbours(CGeometry *geometry, unsigned long iPoint, unsign } -void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDomain, unsigned short val_nVar, unsigned short val_nEq, unsigned long* val_row_ptr, unsigned long* val_col_ind, unsigned long val_nnz, CConfig *config) { +template +void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDomain, unsigned short val_nVar, unsigned short val_nEq, unsigned long* val_row_ptr, unsigned long* val_col_ind, unsigned long val_nnz, CConfig *config) { unsigned long iVar; @@ -348,15 +353,15 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDo nnz_ilu = val_nnz; // Assign number of possible non zero blocks in the spare system structure (ILU structure) } - matrix = new su2double [nnz*nVar*nEqn]; // Reserve memory for the values of the matrix - block = new su2double [nVar*nEqn]; - block_weight = new su2double [nVar*nEqn]; - block_inverse = new su2double [nVar*nEqn]; + matrix = new ScalarType [nnz*nVar*nEqn]; // Reserve memory for the values of the matrix + block = new ScalarType [nVar*nEqn]; + block_weight = new ScalarType [nVar*nEqn]; + block_inverse = new ScalarType [nVar*nEqn]; - prod_block_vector = new su2double [nEqn]; - prod_row_vector = new su2double [nVar]; - aux_vector = new su2double [nVar]; - sum_vector = new su2double [nVar]; + prod_block_vector = new ScalarType [nEqn]; + prod_row_vector = new ScalarType [nVar]; + aux_vector = new ScalarType [nVar]; + sum_vector = new ScalarType [nVar]; /*--- Memory initialization ---*/ @@ -383,7 +388,7 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDo /*--- Reserve memory for the ILU matrix. ---*/ - ILU_matrix = new su2double [nnz_ilu*nVar*nEqn]; + ILU_matrix = new ScalarType [nnz_ilu*nVar*nEqn]; for (iVar = 0; iVar < nnz_ilu*nVar*nEqn; iVar++) ILU_matrix[iVar] = 0.0; } @@ -403,17 +408,19 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDo /*--- Reserve memory for the values of the inverse of the preconditioner. ---*/ - invM = new su2double [nPoint*nVar*nEqn]; + invM = new ScalarType [nPoint*nVar*nEqn]; for (iVar = 0; iVar < nPoint*nVar*nEqn; iVar++) invM[iVar] = 0.0; } } -void CSysMatrix::InitiateComms(CSysVector & x, - CGeometry *geometry, - CConfig *config, - unsigned short commType) { +template +template +void CSysMatrix::InitiateComms(CSysVector & x, + CGeometry *geometry, + CConfig *config, + unsigned short commType) { /*--- Local variables ---*/ @@ -558,10 +565,12 @@ void CSysMatrix::InitiateComms(CSysVector & x, } -void CSysMatrix::CompleteComms(CSysVector & x, - CGeometry *geometry, - CConfig *config, - unsigned short commType) { +template +template +void CSysMatrix::CompleteComms(CSysVector & x, + CGeometry *geometry, + CConfig *config, + unsigned short commType) { /*--- Local variables ---*/ @@ -621,7 +630,7 @@ void CSysMatrix::CompleteComms(CSysVector & x, /*--- Store the data correctly depending on the quantity. ---*/ for (iVar = 0; iVar < nVar; iVar++) - x[iPoint*nVar+iVar] = bufDRecv[buf_offset+iVar]; + x[iPoint*nVar+iVar] = ActiveAssign(bufDRecv[buf_offset+iVar]); } break; @@ -659,7 +668,7 @@ void CSysMatrix::CompleteComms(CSysVector & x, for (iVar = 0; iVar < nVar; iVar++) - x[iPoint*nVar+iVar] += bufDRecv[buf_offset+iVar]; + x[iPoint*nVar+iVar] += ActiveAssign(bufDRecv[buf_offset+iVar]); } @@ -683,7 +692,8 @@ void CSysMatrix::CompleteComms(CSysVector & x, } -su2double *CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j) { +template +ScalarType *CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j) { unsigned long step = 0, index; @@ -695,7 +705,8 @@ su2double *CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j) { } -su2double CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar) { +template +ScalarType CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar) { unsigned long step = 0, index; @@ -707,75 +718,8 @@ su2double CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j, uns } -void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, su2double **val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = val_block[iVar][jVar]; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = SU2_TYPE::GetValue(val_block[iVar][jVar]); - break; - } - } - -} - -void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, su2double *val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = val_block[iVar*nVar+jVar]; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = SU2_TYPE::GetValue(val_block[iVar*nVar+jVar]); - break; - } - } - -} - -void CSysMatrix::AddBlock(unsigned long block_i, unsigned long block_j, su2double **val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] += val_block[iVar][jVar]; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] += SU2_TYPE::GetValue(val_block[iVar][jVar]); - break; - } - } - -} - -void CSysMatrix::SubtractBlock(unsigned long block_i, unsigned long block_j, su2double **val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] -= val_block[iVar][jVar]; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] -= SU2_TYPE::GetValue(val_block[iVar][jVar]); - break; - } - } - -} - -su2double *CSysMatrix::GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j) { +template +ScalarType *CSysMatrix::GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j) { unsigned long step = 0, index; @@ -787,7 +731,8 @@ su2double *CSysMatrix::GetBlock_ILUMatrix(unsigned long block_i, unsigned long b } -void CSysMatrix::SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block) { +template +void CSysMatrix::SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { unsigned long iVar, jVar, index, step = 0; @@ -803,7 +748,8 @@ void CSysMatrix::SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j } -void CSysMatrix::SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block) { +template +void CSysMatrix::SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { unsigned long iVar, jVar, index, step = 0; @@ -819,7 +765,8 @@ void CSysMatrix::SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned lo } -void CSysMatrix::SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block) { +template +void CSysMatrix::SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { unsigned long iVar, jVar, index, step = 0; @@ -835,7 +782,8 @@ void CSysMatrix::SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long bl } -void CSysMatrix::MatrixVectorProduct(su2double *matrix, su2double *vector, su2double *product) { +template +void CSysMatrix::MatrixVectorProduct(ScalarType *matrix, ScalarType *vector, ScalarType *product) { #if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) // NOTE: matrix/vector swapped due to column major kernel -- manual "CBLAS" setup. @@ -857,7 +805,8 @@ void CSysMatrix::MatrixVectorProduct(su2double *matrix, su2double *vector, su2do } -void CSysMatrix::MatrixMatrixProduct(su2double *matrix_a, su2double *matrix_b, su2double *product) { +template +void CSysMatrix::MatrixMatrixProduct(ScalarType *matrix_a, ScalarType *matrix_b, ScalarType *product) { #if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) if (useMKL) @@ -880,45 +829,8 @@ void CSysMatrix::MatrixMatrixProduct(su2double *matrix_a, su2double *matrix_b, s } -void CSysMatrix::AddVal2Diag(unsigned long block_i, su2double val_matrix) { - - unsigned long step = 0, iVar, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_i) { // Only elements on the diagonal - for (iVar = 0; iVar < nVar; iVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] += val_matrix; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] += SU2_TYPE::GetValue(val_matrix); - break; - } - } - -} - -void CSysMatrix::SetVal2Diag(unsigned long block_i, su2double val_matrix) { - - unsigned long step = 0, iVar, jVar, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_i) { // Only elements on the diagonal - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+jVar] = 0.0; - - for (iVar = 0; iVar < nVar; iVar++) -// matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] = val_matrix; // Allow AD in Matrix Structure (disabled temporarily to avoid conflicts) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] = SU2_TYPE::GetValue(val_matrix); - - break; - } - } - -} - -void CSysMatrix::DeleteValsRowi(unsigned long i) { +template +void CSysMatrix::DeleteValsRowi(unsigned long i) { unsigned long block_i = i/nVar; unsigned long row = i - block_i*nVar; @@ -933,11 +845,12 @@ void CSysMatrix::DeleteValsRowi(unsigned long i) { } -su2double CSysMatrix::MatrixDeterminant(su2double **a, unsigned long n) { +template +ScalarType CSysMatrix::MatrixDeterminant(ScalarType **a, unsigned long n) { unsigned long i, j, j1, j2; - su2double det = 0; - su2double **m = NULL; + ScalarType det = 0; + ScalarType **m = NULL; if (n < 1) { } else if (n == 1) { det = a[0][0]; } @@ -946,9 +859,9 @@ su2double CSysMatrix::MatrixDeterminant(su2double **a, unsigned long n) { det = 0.0; for (j1=0;j1 +void CSysMatrix::MatrixCoFactor(ScalarType **a, unsigned long n, ScalarType **b) { unsigned long i,j,ii,jj,i1,j1; - su2double det; - su2double **c; + ScalarType det; + ScalarType **c; - c = new su2double*[n-1]; + c = new ScalarType*[n-1]; for (i=0;i +void CSysMatrix::MatrixTranspose(ScalarType **a, unsigned long n) { unsigned long i, j; - su2double tmp; + ScalarType tmp; for (i=1;i +void CSysMatrix::Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed) { short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - su2double weight, aux; + ScalarType weight, aux; - su2double *Block = GetBlock(block_i, block_i); + ScalarType *Block = GetBlock(block_i, block_i); /*--- Copy block matrix, note that the original matrix is modified by the algorithm---*/ @@ -1079,19 +995,20 @@ void CSysMatrix::Gauss_Elimination(unsigned long block_i, su2double* rhs, bool t } -void CSysMatrix::Gauss_Elimination_ILUMatrix(unsigned long block_i, su2double* rhs) { +template +void CSysMatrix::Gauss_Elimination_ILUMatrix(unsigned long block_i, ScalarType* rhs) { short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - su2double weight, aux; + ScalarType weight, aux; - su2double *Block = GetBlock_ILUMatrix(block_i, block_i); + ScalarType *Block = GetBlock_ILUMatrix(block_i, block_i); /*--- Copy block matrix, note that the original matrix is modified by the algorithm---*/ // If source and dest overlap higher level problems occur, so memcpy is safe. And it is faster. - memcpy( block, Block, (nVar * nVar * sizeof(su2double)) ); + memcpy( block, Block, (nVar * nVar * sizeof(ScalarType)) ); //for (iVar = 0; iVar < (short)nVar; iVar++) // for (jVar = 0; jVar < (short)nVar; jVar++) @@ -1141,10 +1058,11 @@ void CSysMatrix::Gauss_Elimination_ILUMatrix(unsigned long block_i, su2double* r } -void CSysMatrix::Gauss_Elimination(su2double* Block, su2double* rhs) { +template +void CSysMatrix::Gauss_Elimination(ScalarType* Block, ScalarType* rhs) { short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - su2double weight, aux; + ScalarType weight, aux; /*--- Copy block matrix, note that the original matrix is modified by the algorithm---*/ @@ -1181,12 +1099,13 @@ void CSysMatrix::Gauss_Elimination(su2double* Block, su2double* rhs) { } -void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec) { +template +void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec) { unsigned long j = block_j*nVar; unsigned short iVar, jVar; - su2double *block = GetBlock(block_i, block_j); + ScalarType *block = GetBlock(block_i, block_j); for (iVar = 0; iVar < nVar; iVar++) { prod_block_vector[iVar] = 0; @@ -1196,7 +1115,8 @@ void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, c } -void CSysMatrix::UpperProduct(CSysVector & vec, unsigned long row_i) { +template +void CSysMatrix::UpperProduct(CSysVector & vec, unsigned long row_i) { unsigned long iVar, index; @@ -1213,7 +1133,8 @@ void CSysMatrix::UpperProduct(CSysVector & vec, unsigned long row_i) { } -void CSysMatrix::LowerProduct(CSysVector & vec, unsigned long row_i) { +template +void CSysMatrix::LowerProduct(CSysVector & vec, unsigned long row_i) { unsigned long iVar, index; @@ -1230,7 +1151,8 @@ void CSysMatrix::LowerProduct(CSysVector & vec, unsigned long row_i) { } -void CSysMatrix::DiagonalProduct(CSysVector & vec, unsigned long row_i) { +template +void CSysMatrix::DiagonalProduct(CSysVector & vec, unsigned long row_i) { unsigned long iVar, index; @@ -1247,11 +1169,13 @@ void CSysMatrix::DiagonalProduct(CSysVector & vec, unsigned long row_i) { } -void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config) { +template +template +void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config) { unsigned short iVar, iMarker, MarkerS, MarkerR; unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double *Buffer_Receive = NULL, *Buffer_Send = NULL; + OtherType *Buffer_Receive = NULL, *Buffer_Send = NULL; #ifdef HAVE_MPI int send_to, receive_from; @@ -1277,8 +1201,8 @@ void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConf /*--- Allocate Receive and send buffers ---*/ - Buffer_Receive = new su2double [nBufferR_Vector]; - Buffer_Send = new su2double[nBufferS_Vector]; + Buffer_Receive = new OtherType [nBufferR_Vector]; + Buffer_Send = new OtherType[nBufferS_Vector]; /*--- Copy the solution that should be sended ---*/ @@ -1292,7 +1216,7 @@ void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConf /*--- Send/Receive information using Sendrecv ---*/ - SU2_MPI::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, + SelectMPIWrapper::W::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, Buffer_Receive, nBufferR_Vector, MPI_DOUBLE, receive_from, 0, MPI_COMM_WORLD, &status); #else @@ -1335,11 +1259,12 @@ void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConf } -void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config) { unsigned short iVar, iMarker, MarkerS, MarkerR; unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double *Buffer_Receive = NULL, *Buffer_Send = NULL; + ScalarType *Buffer_Receive = NULL, *Buffer_Send = NULL; #ifdef HAVE_MPI int send_to, receive_from; @@ -1365,8 +1290,8 @@ void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geome /*--- Allocate Receive and send buffers ---*/ - Buffer_Receive = new su2double [nBufferR_Vector]; - Buffer_Send = new su2double[nBufferS_Vector]; + Buffer_Receive = new ScalarType [nBufferR_Vector]; + Buffer_Send = new ScalarType[nBufferS_Vector]; /*--- Copy the solution that should be sended ---*/ @@ -1380,7 +1305,7 @@ void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geome /*--- Send/Receive information using Sendrecv ---*/ - SU2_MPI::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, + SelectMPIWrapper::W::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, Buffer_Receive, nBufferR_Vector, MPI_DOUBLE, receive_from, 0, MPI_COMM_WORLD, &status); #else @@ -1423,7 +1348,8 @@ void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geome } -void CSysMatrix::RowProduct(const CSysVector & vec, unsigned long row_i) { +template +void CSysMatrix::RowProduct(const CSysVector & vec, unsigned long row_i) { unsigned long iVar, index; @@ -1438,7 +1364,8 @@ void CSysMatrix::RowProduct(const CSysVector & vec, unsigned long row_i) { } -void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod) { +template +void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod) { unsigned long iPoint, iVar; @@ -1450,23 +1377,24 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod) } -void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long prod_begin, vec_begin, mat_begin, index, iVar, jVar, row_i; - /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ + /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ if ( (nVar != vec.GetNVar()) || (nVar != prod.GetNVar()) ) { - cerr << "CSysMatrix::MatrixVectorProduct(const CSysVector&, CSysVector): " + cerr << "CSysMatrix::MatrixVectorProduct(const CSysVector&, CSysVector): " << "nVar values incompatible." << endl; throw(-1); } if ( (nPoint != vec.GetNBlk()) || (nPoint != prod.GetNBlk()) ) { - cerr << "CSysMatrix::MatrixVectorProduct(const CSysVector&, CSysVector): " + cerr << "CSysMatrix::MatrixVectorProduct(const CSysVector&, CSysVector): " << "nPoint and nBlk values incompatible." << endl; throw(-1); } - prod = su2double(0.0); // set all entries of prod to zero + prod = ScalarType(0.0); // set all entries of prod to zero for (row_i = 0; row_i < nPointDomain; row_i++) { prod_begin = row_i*nVar; // offset to beginning of block row_i for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { @@ -1494,11 +1422,12 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, } -void CSysMatrix::MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long prod_begin, vec_begin, mat_begin, index, iVar, jVar , row_i; - /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ + /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ if ( (nVar != vec.GetNVar()) || (nVar != prod.GetNVar()) ) { SU2_MPI::Error("nVar values incompatible.", CURRENT_FUNCTION); } @@ -1506,7 +1435,7 @@ void CSysMatrix::MatrixVectorProductTransposed(const CSysVector & vec, CSysVecto SU2_MPI::Error("nPoint and nBlk values incompatible.", CURRENT_FUNCTION); } - prod = su2double(0.0); // set all entries of prod to zero + prod = ScalarType(0.0); // set all entries of prod to zero for (row_i = 0; row_i < nPointDomain; row_i++) { vec_begin = row_i*nVar; // offset to beginning of block col_ind[index] for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { @@ -1527,7 +1456,8 @@ void CSysMatrix::MatrixVectorProductTransposed(const CSysVector & vec, CSysVecto } -void CSysMatrix::GetMultBlockBlock(su2double *c, su2double *a, su2double *b) { +template +void CSysMatrix::GetMultBlockBlock(ScalarType *c, ScalarType *a, ScalarType *b) { unsigned long iVar, jVar, kVar; @@ -1540,7 +1470,8 @@ void CSysMatrix::GetMultBlockBlock(su2double *c, su2double *a, su2double *b) { } -void CSysMatrix::GetMultBlockVector(su2double *c, su2double *a, su2double *b) { +template +void CSysMatrix::GetMultBlockVector(ScalarType *c, ScalarType *a, ScalarType *b) { unsigned long iVar, jVar; @@ -1552,7 +1483,8 @@ void CSysMatrix::GetMultBlockVector(su2double *c, su2double *a, su2double *b) { } -void CSysMatrix::GetSubsBlock(su2double *c, su2double *a, su2double *b) { +template +void CSysMatrix::GetSubsBlock(ScalarType *c, ScalarType *a, ScalarType *b) { unsigned long iVar, jVar; @@ -1562,7 +1494,8 @@ void CSysMatrix::GetSubsBlock(su2double *c, su2double *a, su2double *b) { } -void CSysMatrix::GetSubsVector(su2double *c, su2double *a, su2double *b) { +template +void CSysMatrix::GetSubsVector(ScalarType *c, ScalarType *a, ScalarType *b) { unsigned long iVar; @@ -1571,7 +1504,8 @@ void CSysMatrix::GetSubsVector(su2double *c, su2double *a, su2double *b) { } -void CSysMatrix::InverseBlock(su2double *Block, su2double *invBlock) { +template +void CSysMatrix::InverseBlock(ScalarType *Block, ScalarType *invBlock) { unsigned long iVar, jVar; @@ -1589,7 +1523,8 @@ void CSysMatrix::InverseBlock(su2double *Block, su2double *invBlock) { } -void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, su2double *invBlock, bool transpose) { +template +void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose) { unsigned long iVar, jVar; @@ -1605,14 +1540,14 @@ void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, su2double *invBlock invBlock[jVar*nVar+iVar] = aux_vector[jVar]; } - // su2double Det, **Matrix, **CoFactor; - // su2double *Block = GetBlock(block_i, block_i); + // ScalarType Det, **Matrix, **CoFactor; + // ScalarType *Block = GetBlock(block_i, block_i); // - // Matrix = new su2double*[nVar]; - // CoFactor = new su2double*[nVar]; + // Matrix = new ScalarType*[nVar]; + // CoFactor = new ScalarType*[nVar]; // for (iVar=0;iVar +void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock) { unsigned long iVar, jVar; @@ -1656,14 +1591,14 @@ void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, su2double invBlock[jVar*nVar+iVar] = aux_vector[jVar]; } - // su2double Det, **Matrix, **CoFactor; - // su2double *Block = GetBlock_ILUMatrix(block_i, block_i); + // ScalarType Det, **Matrix, **CoFactor; + // ScalarType *Block = GetBlock_ILUMatrix(block_i, block_i); // - // Matrix = new su2double*[nVar]; - // CoFactor = new su2double*[nVar]; + // Matrix = new ScalarType*[nVar]; + // CoFactor = new ScalarType*[nVar]; // for (iVar=0;iVar +void CSysMatrix::BuildJacobiPreconditioner(bool transpose) { unsigned long iPoint, iVar, jVar; @@ -1708,8 +1644,8 @@ void CSysMatrix::BuildJacobiPreconditioner(bool transpose) { } - -void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long iPoint, iVar, jVar; @@ -1729,10 +1665,11 @@ void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector } -unsigned long CSysMatrix::Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config) { +template +unsigned long CSysMatrix::Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { unsigned long iPoint, iVar, jVar; - su2double norm_r = 0.0, norm0 = 0.0; + ScalarType norm_r = 0.0, norm0 = 0.0; int i = 0; /*--- Check the number of iterations requested ---*/ @@ -1747,8 +1684,8 @@ unsigned long CSysMatrix::Jacobi_Smoother(const CSysVector & b, CSysVector & x, of the Jacobian matrix with the current solution (x^k). These must be stored in order to perform multiple iterations of the smoother. ---*/ - CSysVector r(b); - CSysVector A_x(b); + CSysVector r(b); + CSysVector A_x(b); /*--- Calculate the initial residual, compute norm, and check if system is already solved. Recall, r holds b initially. ---*/ @@ -1834,10 +1771,11 @@ unsigned long CSysMatrix::Jacobi_Smoother(const CSysVector & b, CSysVector & x, } -void CSysMatrix::BuildILUPreconditioner(bool transposed) { +template +void CSysMatrix::BuildILUPreconditioner(bool transposed) { unsigned long index, index_, iVar; - su2double *Block_ij, *Block_jk; + ScalarType *Block_ij, *Block_jk; long iPoint, jPoint, kPoint; @@ -1917,10 +1855,11 @@ void CSysMatrix::BuildILUPreconditioner(bool transposed) { } -void CSysMatrix::ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long index; - su2double *Block_ij; + ScalarType *Block_ij; long iPoint, jPoint; unsigned short iVar; @@ -1982,13 +1921,14 @@ void CSysMatrix::ComputeILUPreconditioner(const CSysVector & vec, CSysVector & p } -unsigned long CSysMatrix::ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config) { +template +unsigned long CSysMatrix::ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { unsigned long index; - su2double *Block_ij, omega = 1.0; + ScalarType *Block_ij, omega = 1.0; long iPoint, jPoint; unsigned short iVar; - su2double norm_r = 0.0, norm0 = 0.0; + ScalarType norm_r = 0.0, norm0 = 0.0; int i = 0; /*--- Check the number of iterations requested ---*/ @@ -2003,8 +1943,8 @@ unsigned long CSysMatrix::ILU_Smoother(const CSysVector & b, CSysVector & x, CMa of the Jacobian matrix with the current solution (x^k). These must be stored in order to perform multiple iterations of the smoother. ---*/ - CSysVector r(b); - CSysVector A_x(b); + CSysVector r(b); + CSysVector A_x(b); /*--- Calculate the initial residual, compute norm, and check if system is already solved. Recall, r holds b initially. ---*/ @@ -2141,7 +2081,8 @@ unsigned long CSysMatrix::ILU_Smoother(const CSysVector & b, CSysVector & x, CMa } -void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { +template +void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long iPoint, iVar; /*--- First part of the symmetric iteration: (D+L).x* = b ---*/ @@ -2181,11 +2122,12 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector } -unsigned long CSysMatrix::LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config) { +template +unsigned long CSysMatrix::LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { unsigned long iPoint, iVar; - su2double omega = 1.0; - su2double norm_r = 0.0, norm0 = 0.0; + ScalarType omega = 1.0; + ScalarType norm_r = 0.0, norm0 = 0.0; int i = 0; /*--- Check the number of iterations requested ---*/ @@ -2200,9 +2142,9 @@ unsigned long CSysMatrix::LU_SGS_Smoother(const CSysVector & b, CSysVector & x, of the Jacobian matrix with the current solution (x^k). These must be stored in order to perform multiple iterations of the smoother. ---*/ - CSysVector r(b); - CSysVector A_x(b); - CSysVector xStar(x); + CSysVector r(b); + CSysVector A_x(b); + CSysVector xStar(x); /*--- Calculate the initial residual, compute norm, and check if system is already solved. Recall, r holds b initially. ---*/ @@ -2313,7 +2255,8 @@ unsigned long CSysMatrix::LU_SGS_Smoother(const CSysVector & b, CSysVector & x, } -unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geometry, CConfig *config) { +template +unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geometry, CConfig *config) { bool *check_Point, add_point; unsigned long iEdge, iPoint, jPoint, index_Point, iLinelet, iVertex, next_Point, counter, iElem; @@ -2474,28 +2417,28 @@ unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geometry, CConf SU2_MPI::Allreduce(&Local_nLineLets, &Global_nLineLets, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); #endif - MeanPoints = SU2_TYPE::Int(su2double(Global_nPoints)/su2double(Global_nLineLets)); + MeanPoints = SU2_TYPE::Int(ScalarType(Global_nPoints)/ScalarType(Global_nLineLets)); /*--- Memory allocation --*/ - UBlock = new su2double* [max_nElem]; - invUBlock = new su2double* [max_nElem]; - LBlock = new su2double* [max_nElem]; - yVector = new su2double* [max_nElem]; - zVector = new su2double* [max_nElem]; - rVector = new su2double* [max_nElem]; + UBlock = new ScalarType* [max_nElem]; + invUBlock = new ScalarType* [max_nElem]; + LBlock = new ScalarType* [max_nElem]; + yVector = new ScalarType* [max_nElem]; + zVector = new ScalarType* [max_nElem]; + rVector = new ScalarType* [max_nElem]; for (iElem = 0; iElem < max_nElem; iElem++) { - UBlock[iElem] = new su2double [nVar*nVar]; - invUBlock[iElem] = new su2double [nVar*nVar]; - LBlock[iElem] = new su2double [nVar*nVar]; - yVector[iElem] = new su2double [nVar]; - zVector[iElem] = new su2double [nVar]; - rVector[iElem] = new su2double [nVar]; + UBlock[iElem] = new ScalarType [nVar*nVar]; + invUBlock[iElem] = new ScalarType [nVar*nVar]; + LBlock[iElem] = new ScalarType [nVar*nVar]; + yVector[iElem] = new ScalarType [nVar]; + zVector[iElem] = new ScalarType [nVar]; + rVector[iElem] = new ScalarType [nVar]; } - LFBlock = new su2double [nVar*nVar]; - LyVector = new su2double [nVar]; - FzVector = new su2double [nVar]; + LFBlock = new ScalarType [nVar*nVar]; + LyVector = new ScalarType [nVar]; + FzVector = new ScalarType [nVar]; /*--- Memory deallocation --*/ @@ -2505,12 +2448,13 @@ unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geometry, CConf } -void CSysMatrix::ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, +template +void CSysMatrix::ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long iVar, jVar, nElem = 0, iLinelet, im1Point, iPoint, ip1Point, iElem; long iElemLoop; - su2double *block; + ScalarType *block; if (size == SINGLE_NODE) { @@ -2610,7 +2554,8 @@ void CSysMatrix::ComputeLineletPreconditioner(const CSysVector & vec, CSysVector } -void CSysMatrix::ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res) { +template +void CSysMatrix::ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res) { unsigned long iPoint, iVar; @@ -2622,3 +2567,31 @@ void CSysMatrix::ComputeResidual(const CSysVector & sol, const CSysVector & f, C } } + +/*--- Explicit instantiations ---*/ +template class CSysMatrix; +template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); +template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template class CSysMatrixVectorProduct; +template class CSysMatrixVectorProductTransposed; +template class CJacobiPreconditioner; +template class CILUPreconditioner; +template class CLU_SGSPreconditioner; +template class CLineletPreconditioner; + +#ifdef CODI_REVERSE_TYPE +template class CSysMatrix; +template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); +template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); +template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); +template class CSysMatrixVectorProduct; +template class CSysMatrixVectorProductTransposed; +template class CJacobiPreconditioner; +template class CILUPreconditioner; +template class CLU_SGSPreconditioner; +template class CLineletPreconditioner; +#endif diff --git a/Common/src/vector_structure.cpp b/Common/src/vector_structure.cpp index 289f0cf88dd0..77ac968741f6 100644 --- a/Common/src/vector_structure.cpp +++ b/Common/src/vector_structure.cpp @@ -37,13 +37,15 @@ #include "../include/vector_structure.hpp" -CSysVector::CSysVector(void) { +template +CSysVector::CSysVector(void) { vec_val = NULL; } -CSysVector::CSysVector(const unsigned long & size, const su2double & val) { +template +CSysVector::CSysVector(const unsigned long & size, const ScalarType & val) { nElm = size; nElmDomain = size; nBlk = nElm; nBlkDomain = nElmDomain; @@ -56,7 +58,7 @@ CSysVector::CSysVector(const unsigned long & size, const su2double & val) { SU2_MPI::Error(string(buf), CURRENT_FUNCTION); } - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned int i = 0; i < nElm; i++) vec_val[i] = val; @@ -67,8 +69,9 @@ CSysVector::CSysVector(const unsigned long & size, const su2double & val) { } -CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, - const su2double & val) { +template +CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, + const ScalarType & val) { nElm = numBlk*numVar; nElmDomain = numBlkDomain*numVar; nBlk = numBlk; nBlkDomain = numBlkDomain; @@ -81,7 +84,7 @@ CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBl SU2_MPI::Error(string(buf), CURRENT_FUNCTION); } - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned int i = 0; i < nElm; i++) vec_val[i] = val; @@ -92,14 +95,15 @@ CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBl } -CSysVector::CSysVector(const CSysVector & u) { +template +CSysVector::CSysVector(const CSysVector & u) { /*--- Copy size information, allocate memory, and initialize values ---*/ nElm = u.nElm; nElmDomain = u.nElmDomain; nBlk = u.nBlk; nBlkDomain = u.nBlkDomain; nVar = u.nVar; - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned long i = 0; i < nElm; i++) vec_val[i] = u.vec_val[i]; @@ -109,7 +113,8 @@ CSysVector::CSysVector(const CSysVector & u) { } -CSysVector::CSysVector(const unsigned long & size, const su2double* u_array) { +template +CSysVector::CSysVector(const unsigned long & size, const ScalarType* u_array) { nElm = size; nElmDomain = size; nBlk = nElm; nBlkDomain = nElmDomain; @@ -122,7 +127,7 @@ CSysVector::CSysVector(const unsigned long & size, const su2double* u_array) { SU2_MPI::Error(string(buf), CURRENT_FUNCTION); } - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned long i = 0; i < nElm; i++) vec_val[i] = u_array[i]; @@ -133,8 +138,9 @@ CSysVector::CSysVector(const unsigned long & size, const su2double* u_array) { } -CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, - const su2double* u_array) { +template +CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, + const ScalarType* u_array) { nElm = numBlk*numVar; nElmDomain = numBlkDomain*numVar; nBlk = numBlk; nBlkDomain = numBlkDomain; @@ -147,7 +153,7 @@ CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBl SU2_MPI::Error(string(buf), CURRENT_FUNCTION); } - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned long i = 0; i < nElm; i++) vec_val[i] = u_array[i]; @@ -158,7 +164,8 @@ CSysVector::CSysVector(const unsigned long & numBlk, const unsigned long & numBl } -CSysVector::~CSysVector() { +template +CSysVector::~CSysVector() { delete [] vec_val; nElm = 0; nElmDomain = 0; @@ -167,7 +174,8 @@ CSysVector::~CSysVector() { } -void CSysVector::Initialize(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const su2double & val) { +template +void CSysVector::Initialize(const unsigned long & numBlk, const unsigned long & numBlkDomain, const unsigned short & numVar, const ScalarType & val) { nElm = numBlk*numVar; nElmDomain = numBlkDomain*numVar; nBlk = numBlk; nBlkDomain = numBlkDomain; @@ -180,7 +188,7 @@ void CSysVector::Initialize(const unsigned long & numBlk, const unsigned long & SU2_MPI::Error(string(buf), CURRENT_FUNCTION); } - vec_val = new su2double[nElm]; + vec_val = new ScalarType[nElm]; for (unsigned long i = 0; i < nElm; i++) vec_val[i] = val; @@ -191,7 +199,8 @@ void CSysVector::Initialize(const unsigned long & numBlk, const unsigned long & } -void CSysVector::Equals_AX(const su2double & a, CSysVector & x) { +template +void CSysVector::Equals_AX(const ScalarType & a, CSysVector & x) { /*--- check that *this and x are compatible ---*/ if (nElm != x.nElm) { cerr << "CSysVector::Equals_AX(): " << "sizes do not match"; @@ -201,7 +210,8 @@ void CSysVector::Equals_AX(const su2double & a, CSysVector & x) { vec_val[i] = a * x.vec_val[i]; } -void CSysVector::Plus_AX(const su2double & a, CSysVector & x) { +template +void CSysVector::Plus_AX(const ScalarType & a, CSysVector & x) { /*--- check that *this and x are compatible ---*/ if (nElm != x.nElm) { SU2_MPI::Error("Sizes do not match", CURRENT_FUNCTION); @@ -210,7 +220,8 @@ void CSysVector::Plus_AX(const su2double & a, CSysVector & x) { vec_val[i] += a * x.vec_val[i]; } -void CSysVector::Equals_AX_Plus_BY(const su2double & a, CSysVector & x, const su2double & b, CSysVector & y) { +template +void CSysVector::Equals_AX_Plus_BY(const ScalarType & a, CSysVector & x, const ScalarType & b, CSysVector & y) { /*--- check that *this, x and y are compatible ---*/ if ((nElm != x.nElm) || (nElm != y.nElm)) { SU2_MPI::Error("Sizes do not match", CURRENT_FUNCTION); @@ -219,14 +230,15 @@ void CSysVector::Equals_AX_Plus_BY(const su2double & a, CSysVector & x, const su vec_val[i] = a * x.vec_val[i] + b * y.vec_val[i]; } -CSysVector & CSysVector::operator=(const CSysVector & u) { +template +CSysVector & CSysVector::operator=(const CSysVector & u) { /*--- check if self-assignment, otherwise perform deep copy ---*/ if (this == &u) return *this; /*--- determine if (re-)allocation is needed ---*/ if (nElm != u.nElm && vec_val != NULL) {delete [] vec_val; vec_val = NULL;} - if (vec_val == NULL) vec_val = new su2double[u.nElm]; + if (vec_val == NULL) vec_val = new ScalarType[u.nElm]; /*--- copy ---*/ nElm = u.nElm; @@ -245,21 +257,24 @@ CSysVector & CSysVector::operator=(const CSysVector & u) { return *this; } -CSysVector & CSysVector::operator=(const su2double & val) { +template +CSysVector & CSysVector::operator=(const ScalarType & val) { for (unsigned long i = 0; i < nElm; i++) vec_val[i] = val; return *this; } -CSysVector CSysVector::operator+(const CSysVector & u) const { +template +CSysVector CSysVector::operator+(const CSysVector & u) const { /*--- Use copy constructor and compound addition-assignment ---*/ - CSysVector sum(*this); + CSysVector sum(*this); sum += u; return sum; } -CSysVector & CSysVector::operator+=(const CSysVector & u) { +template +CSysVector & CSysVector::operator+=(const CSysVector & u) { /*--- Check for consistent sizes, then add elements ---*/ if (nElm != u.nElm) { @@ -270,15 +285,17 @@ CSysVector & CSysVector::operator+=(const CSysVector & u) { return *this; } -CSysVector CSysVector::operator-(const CSysVector & u) const { +template +CSysVector CSysVector::operator-(const CSysVector & u) const { /*--- Use copy constructor and compound subtraction-assignment ---*/ - CSysVector diff(*this); + CSysVector diff(*this); diff -= u; return diff; } -CSysVector & CSysVector::operator-=(const CSysVector & u) { +template +CSysVector & CSysVector::operator-=(const CSysVector & u) { /*--- Check for consistent sizes, then subtract elements ---*/ if (nElm != u.nElm) { @@ -289,32 +306,36 @@ CSysVector & CSysVector::operator-=(const CSysVector & u) { return *this; } -CSysVector CSysVector::operator*(const su2double & val) const { +template +CSysVector CSysVector::operator*(const ScalarType & val) const { /*--- use copy constructor and compound scalar multiplication-assignment ---*/ - CSysVector prod(*this); + CSysVector prod(*this); prod *= val; return prod; } -CSysVector operator*(const su2double & val, const CSysVector & u) { +template +CSysVector operator*(const ScalarType & val, const CSysVector & u) { /*--- use copy constructor and compound scalar multiplication-assignment ---*/ - CSysVector prod(u); + CSysVector prod(u); prod *= val; return prod; } -CSysVector & CSysVector::operator*=(const su2double & val) { +template +CSysVector & CSysVector::operator*=(const ScalarType & val) { for (unsigned long i = 0; i < nElm; i++) vec_val[i] *= val; return *this; } -CSysVector CSysVector::operator/(const su2double & val) const { +template +CSysVector CSysVector::operator/(const ScalarType & val) const { /*--- use copy constructor and compound scalar division-assignment ---*/ @@ -323,75 +344,87 @@ CSysVector CSysVector::operator/(const su2double & val) const { return quotient; } -CSysVector & CSysVector::operator/=(const su2double & val) { +template +CSysVector & CSysVector::operator/=(const ScalarType & val) { for (unsigned long i = 0; i < nElm; i++) vec_val[i] /= val; return *this; } -su2double CSysVector::norm() const { +template +ScalarType CSysVector::norm() const { /*--- just call dotProd on this*, then sqrt ---*/ - su2double val = dotProd(*this, *this); + ScalarType val = dotProd(*this, *this); if (val < 0.0) { SU2_MPI::Error("Inner product of CSysVector is negative", CURRENT_FUNCTION); } return sqrt(val); } -void CSysVector::CopyToArray(su2double* u_array) { +template +void CSysVector::CopyToArray(ScalarType* u_array) { for (unsigned long i = 0; i < nElm; i++) u_array[i] = vec_val[i]; } -void CSysVector::AddBlock(unsigned long val_ipoint, su2double *val_residual) { +template +void CSysVector::AddBlock(unsigned long val_ipoint, ScalarType *val_residual) { unsigned short iVar; for (iVar = 0; iVar < nVar; iVar++) vec_val[val_ipoint*nVar+iVar] += val_residual[iVar]; } -void CSysVector::SubtractBlock(unsigned long val_ipoint, su2double *val_residual) { +template +void CSysVector::SubtractBlock(unsigned long val_ipoint, ScalarType *val_residual) { unsigned short iVar; for (iVar = 0; iVar < nVar; iVar++) vec_val[val_ipoint*nVar+iVar] -= val_residual[iVar]; } -void CSysVector::SetBlock(unsigned long val_ipoint, su2double *val_residual) { +template +void CSysVector::SetBlock(unsigned long val_ipoint, ScalarType *val_residual) { unsigned short iVar; for (iVar = 0; iVar < nVar; iVar++) vec_val[val_ipoint*nVar+iVar] = val_residual[iVar]; } -void CSysVector::SetBlock(unsigned long val_ipoint, unsigned short val_var, su2double val_residual) { +template +void CSysVector::SetBlock(unsigned long val_ipoint, unsigned short val_var, ScalarType val_residual) { vec_val[val_ipoint*nVar+val_var] = val_residual; } -void CSysVector::SetBlock_Zero(unsigned long val_ipoint) { +template +void CSysVector::SetBlock_Zero(unsigned long val_ipoint) { unsigned short iVar; for (iVar = 0; iVar < nVar; iVar++) vec_val[val_ipoint*nVar+iVar] = 0.0; } -void CSysVector::SetBlock_Zero(unsigned long val_ipoint, unsigned short val_var) { +template +void CSysVector::SetBlock_Zero(unsigned long val_ipoint, unsigned short val_var) { vec_val[val_ipoint*nVar+val_var] = 0.0; } -su2double CSysVector::GetBlock(unsigned long val_ipoint, unsigned short val_var) { +template +ScalarType CSysVector::GetBlock(unsigned long val_ipoint, unsigned short val_var) { return vec_val[val_ipoint*nVar + val_var]; } -su2double *CSysVector::GetBlock(unsigned long val_ipoint) { +template +ScalarType *CSysVector::GetBlock(unsigned long val_ipoint) { return &vec_val[val_ipoint*nVar]; } -su2double dotProd(const CSysVector & u, const CSysVector & v) { +template +ScalarType dotProd(const CSysVector & u, const CSysVector & v) { /*--- check for consistent sizes ---*/ if (u.nElm != v.nElm) { @@ -400,16 +433,27 @@ su2double dotProd(const CSysVector & u, const CSysVector & v) { /*--- find local inner product and, if a parallel run, sum over all processors (we use nElemDomain instead of nElem) ---*/ - su2double loc_prod = 0.0; + ScalarType loc_prod = 0.0; for (unsigned long i = 0; i < u.nElmDomain; i++) loc_prod += u.vec_val[i]*v.vec_val[i]; - su2double prod = 0.0; + ScalarType prod = 0.0; #ifdef HAVE_MPI - SU2_MPI::Allreduce(&loc_prod, &prod, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + SelectMPIWrapper::W::Allreduce(&loc_prod, &prod, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #else prod = loc_prod; #endif return prod; } + +/*--- Explicit instantiations ---*/ +template class CSysVector; +template CSysVector operator*(const su2double&, const CSysVector&); +template su2double dotProd(const CSysVector & u, const CSysVector & v); + +#ifdef CODI_REVERSE_TYPE +template class CSysVector; +template CSysVector operator*(const passivedouble&, const CSysVector&); +template passivedouble dotProd(const CSysVector & u, const CSysVector & v); +#endif diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp index d248918f746e..7350c77538bc 100644 --- a/SU2_CFD/include/solver_structure.hpp +++ b/SU2_CFD/include/solver_structure.hpp @@ -146,15 +146,20 @@ class CSolver { public: - CSysVector LinSysSol; /*!< \brief vector to store iterative solution of implicit linear system. */ - CSysVector LinSysRes; /*!< \brief vector to store iterative residual of implicit linear system. */ - CSysVector LinSysAux; /*!< \brief vector to store iterative residual of implicit linear system. */ - CSysMatrix Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ - CSysSolve System; /*!< \brief Linear solver/smoother. */ + CSysVector LinSysSol; /*!< \brief vector to store iterative solution of implicit linear system. */ + CSysVector LinSysRes; /*!< \brief vector to store iterative residual of implicit linear system. */ + CSysVector LinSysAux; /*!< \brief vector to store iterative residual of implicit linear system. */ +#ifndef CODI_FORWARD_TYPE + CSysMatrix Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ + CSysSolve System; /*!< \brief Linear solver/smoother. */ +#else + CSysMatrix Jacobian; + CSysSolve System; +#endif - CSysMatrix StiffMatrix; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations, and grid movement. */ + CSysMatrix StiffMatrix; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations, and grid movement. */ - CSysVector OutputVariables; /*!< \brief vector to store the extra variables to be written. */ + CSysVector OutputVariables; /*!< \brief vector to store the extra variables to be written. */ string* OutputHeadingNames; /*< \brief vector of strings to store the headings for the exra variables */ CVariable** node; /*!< \brief Vector which the define the variables for each problem. */ @@ -10353,11 +10358,11 @@ class CTransLMSolver: public CTurbSolver { void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config); // Another set of matrix structures for the Lm equations - CSysMatrix JacobianItmc; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ + CSysMatrix JacobianItmc; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ su2double *LinSysSolItmc; /*!< \brief vector to store iterative solution of implicit linear system. */ su2double *LinSysResItmc; /*!< \brief vector to store iterative residual of implicit linear system. */ su2double *rhsItmc; /*!< \brief right hand side of implicit linear system. */ - CSysMatrix JacobianReth; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ + CSysMatrix JacobianReth; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ su2double *LinSysSolReth; /*!< \brief vector to store iterative solution of implicit linear system. */ su2double *LinSysResReth; /*!< \brief vector to store iterative residual of implicit linear system. */ su2double *rhsReth; /*!< \brief right hand side of implicit linear system. */ @@ -11354,8 +11359,8 @@ class CWaveSolver : public CSolver { su2double AllBound_CWave; /*!< \brief Total wave strength for all the boundaries. */ su2double Total_CWave; /*!< \brief Total wave strength for all the boundaries. */ - CSysMatrix StiffMatrixSpace; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations. */ - CSysMatrix StiffMatrixTime; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations. */ + CSysMatrix StiffMatrixSpace; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations. */ + CSysMatrix StiffMatrixTime; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations. */ su2double **StiffMatrix_Elem, /*!< \brief Auxiliary matrices for storing point to point Stiffness Matrices. */ **StiffMatrix_Node; /*!< \brief Auxiliary matrices for storing point to point Stiffness Matrices. */ @@ -11862,14 +11867,14 @@ class CFEASolver : public CSolver { public: - CSysVector TimeRes_Aux; /*!< \brief Auxiliary vector for adding mass and damping contributions to the residual. */ - CSysVector TimeRes; /*!< \brief Vector for adding mass and damping contributions to the residual */ - CSysVector LinSysReact; /*!< \brief Vector to store the residual before applying the BCs */ + CSysVector TimeRes_Aux; /*!< \brief Auxiliary vector for adding mass and damping contributions to the residual. */ + CSysVector TimeRes; /*!< \brief Vector for adding mass and damping contributions to the residual */ + CSysVector LinSysReact; /*!< \brief Vector to store the residual before applying the BCs */ - CSysVector LinSysSol_Adj; /*!< \brief Vector to store the solution of the adjoint problem */ - CSysVector LinSysRes_Adj; /*!< \brief Vector to store the residual of the adjoint problem */ + CSysVector LinSysSol_Adj; /*!< \brief Vector to store the solution of the adjoint problem */ + CSysVector LinSysRes_Adj; /*!< \brief Vector to store the residual of the adjoint problem */ - CSysMatrix MassMatrix; /*!< \brief Sparse structure for storing the mass matrix. */ + CSysMatrix MassMatrix; /*!< \brief Sparse structure for storing the mass matrix. */ CElement*** element_container; /*!< \brief Vector which the define the finite element structure for each problem. */ CElementProperty** element_properties; /*!< \brief Vector which stores the properties of each element */ diff --git a/SU2_CFD/src/solver_direct_elasticity.cpp b/SU2_CFD/src/solver_direct_elasticity.cpp index 18fa9b51fa48..a7802a9bb8ff 100644 --- a/SU2_CFD/src/solver_direct_elasticity.cpp +++ b/SU2_CFD/src/solver_direct_elasticity.cpp @@ -2749,8 +2749,19 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, else { /*--- If the problem is linear, the only check we do is the RMS of the displacements ---*/ /*--- Compute the residual Ax-f ---*/ - +#ifndef CODI_REVERSE_TYPE Jacobian.ComputeResidual(LinSysSol, LinSysRes, LinSysAux); +#else + /*--- We need temporaries to interface with the matrix ---*/ + { + CSysVector sol, res; + sol.PassiveCopy(LinSysSol); + res.PassiveCopy(LinSysRes); + CSysVector aux(res); + Jacobian.ComputeResidual(sol, res, aux); + LinSysAux.PassiveCopy(aux); + } +#endif /*--- Set maximum residual to zero ---*/ @@ -2858,8 +2869,19 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, /*--- If the problem is linear, the only check we do is the RMS of the displacements ---*/ /*--- Compute the residual Ax-f ---*/ - +#ifndef CODI_REVERSE_TYPE Jacobian.ComputeResidual(LinSysSol, LinSysRes, LinSysAux); +#else + /*--- We need temporaries to interface with the matrix ---*/ + { + CSysVector sol, res; + sol.PassiveCopy(LinSysSol); + res.PassiveCopy(LinSysRes); + CSysVector aux(res); + Jacobian.ComputeResidual(sol, res, aux); + LinSysAux.PassiveCopy(aux); + } +#endif /*--- Set maximum residual to zero ---*/