Use FGCRODR instead of GMRES for "Newton-Krylov" adjoint (#2782)

* checkpoint

* optimize product

* nested parallel improvements

* add test, non flexible mode for FGCRODR

* fix the VkWk optimization

* simplify settings

* update

Author: pcarruscag

Common/include/linear_algebra/CSysSolve.hpp

@@ -41,6 +41,11 @@
 #include "CSysVector.hpp"
 #include "../option_structure.hpp"
 
+SU2_IGNORE_WARNING("-Wmaybe-uninitialized")
+#include "Eigen/Core"
+#include "Eigen/Dense"
+SU2_RESTORE_WARNING
+
 class CConfig;
 class CGeometry;
 template <class T>
@@ -110,6 +115,7 @@ class CSysSolve {
   mutable unsigned long k = 0;
   mutable std::vector<VectorType> Z, V; /*!< \brief Large matrices used by FGMRES, v^i+1 = A * z^i. */
   mutable std::vector<VectorType> W, T; /*!< \brief Large matrices used by FGCRODR for deflation vectors. */
+  mutable Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> VkWk;
 
   /*!< \brief Temporary used when it is necessary to interface between active and passive types. */
   VectorType LinSysSol_tmp;
@@ -292,8 +298,8 @@ class CSysSolve {
   template <class Dummy = int>
   unsigned long FGCRODR_LinSolverImpl(const VectorType& b, VectorType& x, const ProductType& mat_vec,
                                       const PrecondType& precond, ScalarType tol, unsigned long max_iter,
-                                      ScalarType& residual, bool monitoring, const CConfig* config,
-                                      FgcrodrMode mode) const;
+                                      ScalarType& residual, bool monitoring, const CConfig* config, FgcrodrMode mode,
+                                      unsigned long custom_m) const;
 
   /*!
    * \brief Creates the inner solver for nested preconditioning if the settings allow it.
@@ -316,7 +322,7 @@ class CSysSolve {
    * \param[in] tol - tolerance with which to solve the system
    * \param[in] m - maximum size of the search subspace
    * \param[out] residual - final normalized residual
-   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] monitoring - turn on priting residuals from solver to screen
    * \param[in] config - Definition of the particular problem.
    */
   unsigned long CG_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec, const PrecondType& precond,
@@ -332,7 +338,7 @@ class CSysSolve {
    * \param[in] tol - tolerance with which to solve the system
    * \param[in] m - maximum size of the search subspace
    * \param[out] residual - final normalized residual
-   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] monitoring - turn on priting residuals from solver to screen
    * \param[in] config - Definition of the particular problem.
    */
   unsigned long FGMRES_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
@@ -355,14 +361,15 @@ class CSysSolve {
    * \param[in] tol - tolerance with which to solve the system
    * \param[in] max_iter - maximum number of iterations
    * \param[out] residual - final normalized residual
-   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] monitoring - turn on priting residuals from solver to screen
    * \param[in] config - Definition of the particular problem.
    * \param[in] mode - See FgcrodrMode.
+   * \param[in] custom_m - alternative maximum size of the search subspace, overrides the config value if != 0.
    */
   unsigned long FGCRODR_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
                                   const PrecondType& precond, ScalarType tol, unsigned long max_iter,
                                   ScalarType& residual, bool monitoring, const CConfig* config,
-                                  FgcrodrMode mode = FgcrodrMode::NORMAL) const;
+                                  FgcrodrMode mode = FgcrodrMode::NORMAL, unsigned long custom_m = 0) const;
 
   /*!
    * \brief Biconjugate Gradient Stabilized Method (BCGSTAB)
@@ -373,7 +380,7 @@ class CSysSolve {
    * \param[in] tol - tolerance with which to solve the system
    * \param[in] m - maximum size of the search subspace
    * \param[out] residual - final normalized residual
-   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] monitoring - turn on priting residuals from solver to screen
    * \param[in] config - Definition of the particular problem.
    */
   unsigned long BCGSTAB_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
@@ -389,7 +396,7 @@ class CSysSolve {
    * \param[in] tol - tolerance with which to solve the system
    * \param[in] m - maximum number of iterations
    * \param[out] residual - final normalized residual
-   * \param[in] monitoring - turn on priting residuals from solver to screen.
+   * \param[in] monitoring - turn on priting residuals from solver to screen
    * \param[in] config - Definition of the particular problem.
    */
   unsigned long Smoother_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,

Common/include/linear_algebra/CSysVector.hpp

@@ -29,6 +29,7 @@
 #pragma once
 
 #include <memory>
+#include <vector>
 
 #include "../parallelization/mpi_structure.hpp"
 #include "../parallelization/omp_structure.hpp"
@@ -371,6 +372,18 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     return dot_scratch[0];
   }
 
+  /*!
+   * \brief Computes the product of V^T W efficiencly, where V and W are tall matrices stored as vectors of CSysVector.
+   * \param[in] V - Tall matrix.
+   * \param[in] i0 - First column of V to consider.
+   * \param[in] n - Number of columns to consider from V starting at i0.
+   * \param[in] W - Tall matrix.
+   * \param[in] m - Number of columns to consider from W.
+   * \return n by m matrix with the result of the product.
+   */
+  static const su2matrix<ScalarType>& multiDot(const std::vector<CSysVector>& V, size_t i0, size_t n,
+                                               const std::vector<CSysVector>& W, size_t m);
+
   /*!
    * \brief Squared L2 norm of the vector (via dot with self).
    * \return Squared L2 norm.