diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp
index a0cecaa63d7..e6829461f47 100644
--- a/Common/include/linear_algebra/CMatrixVectorProduct.hpp
+++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp
@@ -72,6 +72,7 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
   const CSysMatrix<ScalarType>& matrix; /*!< \brief pointer to matrix that defines the product. */
   CGeometry* geometry;                  /*!< \brief geometry associated with the matrix. */
   const CConfig* config;                /*!< \brief config of the problem. */
+  VecExpr::CDeviceExpressionContext device_context;
   mutable bool matrix_uploaded = false; /*!< \brief Upload the matrix lazily on the first actual GPU matvec. */
 
  public:
@@ -83,7 +84,7 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct<ScalarType> {
    */
   inline CSysMatrixVectorProduct(const CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref,
                                  const CConfig* config_ref)
-      : matrix(matrix_ref), geometry(geometry_ref), config(config_ref) {}
+      : matrix(matrix_ref), geometry(geometry_ref), config(config_ref), device_context(config_ref->GetCUDA()) {}
 
   /*!
    * \note This class cannot be default constructed as that would leave us with invalid pointers.
diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp
index e4fc7cf159f..4eff5f947fb 100644
--- a/Common/include/linear_algebra/CPreconditioner.hpp
+++ b/Common/include/linear_algebra/CPreconditioner.hpp
@@ -77,6 +77,18 @@ class CPreconditioner {
 template <class ScalarType>
 CPreconditioner<ScalarType>::~CPreconditioner() {}
 
+/*!
+ * \class CIdentityPreconditioner
+ * \brief No-op preconditioner used when Krylov solvers run without preconditioning.
+ */
+template <class ScalarType>
+class CIdentityPreconditioner final : public CPreconditioner<ScalarType> {
+ public:
+  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { v = u; }
+
+  inline bool IsIdentity() const override { return true; }
+};
+
 /*!
  * \class CJacobiPreconditioner
  * \brief Specialization of preconditioner that uses CSysMatrix class.
@@ -332,6 +344,9 @@ CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOL
   CPreconditioner<ScalarType>* prec = nullptr;
 
   switch (kind) {
+    case NO_PRECONDITIONER:
+      prec = new CIdentityPreconditioner<ScalarType>();
+      break;
     case JACOBI:
       prec = new CJacobiPreconditioner<ScalarType>(jacobian, geometry, config);
       break;
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
index ecb26959a06..6368670587e 100644
--- a/Common/include/linear_algebra/CSysMatrix.hpp
+++ b/Common/include/linear_algebra/CSysMatrix.hpp
@@ -108,6 +108,16 @@ struct CSysMatrixComms {
                        MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX);
 };
 
+/*!
+ * \brief Opaque storage for CUDA/cuSPARSE resources used by the GPU SpMV path.
+ */
+struct CudaSpMVResources;
+
+/*!
+ * \brief Release cached CUDA/cuSPARSE resources used by the GPU SpMV path.
+ */
+void ReleaseCudaSpMVResources(CudaSpMVResources*& resources);
+
 /*!
  * \class CSysMatrix
  * \ingroup SpLinSys
@@ -148,8 +158,10 @@ class CSysMatrix {
   ScalarType* d_matrix;           /*!< \brief Device Pointer to store the matrix values on the GPU. */
   const unsigned long* d_row_ptr; /*!< \brief Device Pointers to the first element in each row. */
   const unsigned long* d_col_ind; /*!< \brief Device Column index for each of the elements in val(). */
+  ScalarType* d_invM = nullptr;   /*!< \brief Device inverse diagonal blocks for the Jacobi preconditioner. */
   bool useCuda = false;           /*!< \brief Boolean that indicates whether user has enabled CUDA or not.
                                      Mainly used to conditionally free GPU memory in the class destructor. */
+  mutable CudaSpMVResources* spmv_resources = nullptr; /*!< \brief Cached cuSPARSE resources for GPU SpMV. */
 
   ScalarType* ILU_matrix;           /*!< \brief Entries of the ILU sparse matrix. */
   unsigned long nnz_ilu;            /*!< \brief Number of possible nonzero entries in the matrix (ILU). */
@@ -924,6 +936,13 @@ class CSysMatrix {
    */
   void BuildJacobiPreconditioner();
 
+  /*!
+   * \brief Build the Jacobi preconditioner on the GPU/device side.
+   * \note This helper is intended as the implementation hook for GPU-resident Krylov solvers.
+   *       The actual implementation belongs in CSysMatrixGPU.cu.
+   */
+  void BuildJacobiPreconditionerGPU();
+
   /*!
    * \brief Multiply CSysVector by the preconditioner
    * \param[in] vec - CSysVector to be multiplied by the preconditioner.
@@ -934,6 +953,14 @@ class CSysMatrix {
   void ComputeJacobiPreconditioner(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod, CGeometry* geometry,
                                    const CConfig* config) const;
 
+  /*!
+   * \brief Apply the Jacobi preconditioner on the GPU/device side.
+   * \note This helper is intended as the implementation hook for GPU-resident Krylov solvers.
+   *       The actual implementation belongs in CSysMatrixGPU.cu.
+   */
+  void ComputeJacobiPreconditionerGPU(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod,
+                                      CGeometry* geometry, const CConfig* config) const;
+
   /*!
    * \brief Build the ILU preconditioner.
    */
diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp
index 1498b549bbb..aabb1aa69a4 100644
--- a/Common/include/linear_algebra/CSysVector.hpp
+++ b/Common/include/linear_algebra/CSysVector.hpp
@@ -37,6 +37,16 @@
 #include "vector_expressions.hpp"
 #include "../../include/CConfig.hpp"
 
+#ifdef __CUDACC__
+#include "GPUComms.cuh"
+#define SU2_CUDA_HOST_DEVICE __host__ __device__
+#else
+#define SU2_CUDA_HOST_DEVICE
+#endif
+
+template <class ScalarType>
+class CSysVector;
+
 /*!
  * \brief OpenMP worksharing construct used in CSysVector for loops.
  * \note The loop will only run in parallel if methods are called from a
@@ -59,6 +69,132 @@
 #define END_CSYSVEC_PARFOR
 #endif
 
+namespace VecExpr {
+
+enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide };
+
+#ifdef HAVE_CUDA
+bool DeviceExpressionsEnabled();
+void BeginDeviceExpressionContext(bool enabled);
+void EndDeviceExpressionContext();
+void FlushDeviceExpressionContext();
+unsigned long CurrentDeviceExpressionContextId();
+void RegisterDeviceModifiedVector(const void* vector, void (*sync)(const void*));
+void UnregisterDeviceModifiedVector(const void* vector);
+#else
+inline bool DeviceExpressionsEnabled() { return false; }
+inline void BeginDeviceExpressionContext(bool) {}
+inline void EndDeviceExpressionContext() {}
+inline void FlushDeviceExpressionContext() {}
+inline unsigned long CurrentDeviceExpressionContextId() { return 0; }
+inline void RegisterDeviceModifiedVector(const void*, void (*)(const void*)) {}
+inline void UnregisterDeviceModifiedVector(const void*) {}
+#endif
+
+class CDeviceExpressionContext {
+ private:
+  bool enabled = false;
+
+ public:
+  explicit CDeviceExpressionContext(bool enabled_) : enabled(enabled_) { BeginDeviceExpressionContext(enabled); }
+  ~CDeviceExpressionContext() { EndDeviceExpressionContext(); }
+
+  CDeviceExpressionContext(const CDeviceExpressionContext&) = delete;
+  CDeviceExpressionContext& operator=(const CDeviceExpressionContext&) = delete;
+};
+
+template <class Scalar>
+class CVectorView : public CVecExpr<CVectorView<Scalar>, Scalar> {
+ private:
+  const Scalar* data = nullptr;
+
+ public:
+  static constexpr bool StoreAsRef = false;
+
+  CVectorView(const Scalar* data_) : data(data_) {}
+  CVectorView(const CSysVector<Scalar>& vector);
+
+  SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t i) const { return data[i]; }
+};
+
+template <class Scalar>
+struct store_type<CSysVector<Scalar>> {
+  using type = CVectorView<Scalar>;
+};
+
+template <class Scalar>
+struct store_type<const CSysVector<Scalar>> {
+  using type = CVectorView<Scalar>;
+};
+
+template <class Scalar>
+struct is_device_assignable<CSysVector<Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<const CSysVector<Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<CVectorView<Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<Bcast<Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<minus_<CSysVector<Scalar>, Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<
+    add_<mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>>
+    : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<
+    add_<add_<mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>,
+         mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>> : std::true_type {};
+
+template <class Scalar>
+struct is_device_assignable<add_<
+    add_<add_<mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>,
+         mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>,
+    mul_<CSysVector<Scalar>, Bcast<Scalar>, Scalar>, Scalar>> : std::true_type {};
+
+template <DeviceAssignOp Op, class Scalar, class T>
+void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr<T, Scalar>& expr);
+
+#ifdef __CUDACC__
+template <DeviceAssignOp Op, class Scalar, class T>
+__global__ void DeviceAssignKernel(Scalar* data, unsigned long size, T expr) {
+  const unsigned long i = static_cast<unsigned long>(blockIdx.x) * blockDim.x + threadIdx.x;
+  if (i >= size) return;
+
+  if constexpr (Op == DeviceAssignOp::Assign) {
+    data[i] = expr[i];
+  } else if constexpr (Op == DeviceAssignOp::Add) {
+    data[i] += expr[i];
+  } else if constexpr (Op == DeviceAssignOp::Subtract) {
+    data[i] -= expr[i];
+  } else if constexpr (Op == DeviceAssignOp::Multiply) {
+    data[i] *= expr[i];
+  } else {
+    data[i] /= expr[i];
+  }
+}
+
+template <DeviceAssignOp Op, class Scalar, class T>
+inline void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr<T, Scalar>& expr) {
+  if (size == 0) return;
+  constexpr unsigned block_size = 256;
+  const auto grid_size = static_cast<unsigned>((size + block_size - 1) / block_size);
+  DeviceAssignKernel<Op><<<grid_size, block_size>>>(data, size, expr.derived());
+  gpuErrChk(cudaPeekAtLastError());
+}
+#endif
+
+}  // namespace VecExpr
+
 /*!
  * \class CSysVector
  * \ingroup SpLinSys
@@ -77,6 +213,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   unsigned long nVar = 1;       /*!< \brief Number of elements in a block. */
 
   ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */
+  mutable bool device_data_valid = false;
+  mutable bool host_data_valid = true;
+  mutable unsigned long device_context_id = 0;
 
 #ifdef HAVE_OMP
   mutable std::unique_ptr<ScalarType[]>
@@ -110,6 +249,34 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     }
   }
 
+  template <VecExpr::DeviceAssignOp Op, class T>
+  CSysVector& AssignDevice(const VecExpr::CVecExpr<T, ScalarType>& expr) {
+#ifdef HAVE_CUDA
+    if constexpr (Op != VecExpr::DeviceAssignOp::Assign) EnsureDeviceData();
+    VecExpr::store_t<const T> stored_expr(expr.derived());
+    VecExpr::AssignDeviceExpression<Op>(d_vec_val, nElm, stored_expr);
+    MarkDeviceDataModified();
+#endif
+    return *this;
+  }
+
+  template <VecExpr::DeviceAssignOp Op>
+  CSysVector& AssignDevice(ScalarType val) {
+#ifdef HAVE_CUDA
+    if constexpr (Op != VecExpr::DeviceAssignOp::Assign) EnsureDeviceData();
+    VecExpr::AssignDeviceExpression<Op>(d_vec_val, nElm, VecExpr::Bcast<ScalarType>(val));
+    MarkDeviceDataModified();
+#endif
+    return *this;
+  }
+
+  void MarkHostDataModified() const {
+#ifdef HAVE_CUDA
+    host_data_valid = true;
+    device_data_valid = false;
+#endif
+  }
+
  public:
   static constexpr bool StoreAsRef = true; /*! \brief Required by CVecExpr. */
 
@@ -173,6 +340,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     std::swap(omp_chunk_size, other.omp_chunk_size);
     std::swap(vec_val, other.vec_val);
     std::swap(d_vec_val, other.d_vec_val);
+    std::swap(device_data_valid, other.device_data_valid);
+    std::swap(host_data_valid, other.host_data_valid);
+    std::swap(device_context_id, other.device_context_id);
     std::swap(nElm, other.nElm);
     std::swap(nElmDomain, other.nElmDomain);
     std::swap(nVar, other.nVar);
@@ -220,6 +390,8 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]);
     END_CSYSVEC_PARFOR
+
+    MarkHostDataModified();
   }
 
   /*!
@@ -234,17 +406,66 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   void DtHTransfer(bool trigger = true) const;
 
+  void EnsureDeviceData() const {
+#ifdef HAVE_CUDA
+    const auto context_id = VecExpr::CurrentDeviceExpressionContextId();
+    if (VecExpr::DeviceExpressionsEnabled() && (!device_data_valid || device_context_id != context_id)) {
+      HtDTransfer();
+      device_data_valid = true;
+      host_data_valid = true;
+      device_context_id = context_id;
+    }
+#endif
+  }
+
+  void MarkDeviceDataModified() const {
+#ifdef HAVE_CUDA
+    if (VecExpr::DeviceExpressionsEnabled()) {
+      device_data_valid = true;
+      host_data_valid = false;
+      device_context_id = VecExpr::CurrentDeviceExpressionContextId();
+      VecExpr::RegisterDeviceModifiedVector(
+          this, [](const void* vector) { static_cast<const CSysVector*>(vector)->SyncHostFromDevice(); });
+    }
+#endif
+  }
+
+  void SyncHostFromDevice() const {
+#ifdef HAVE_CUDA
+    if (device_data_valid && !host_data_valid) {
+      DtHTransfer();
+      host_data_valid = true;
+    }
+#endif
+  }
+
   /*!
-   * \brief Sets all the elements of the GPU vector to a certain value
-   * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
+   * \brief Dot product between this vector and another vector on the device.
+   * \note Explicit GPU helper for solver-side reductions.
+   * \param[in] other - Input vector.
+   * \return Dot product result.
+   */
+  ScalarType GPUDot(const CSysVector& other) const;
+
+  /*!
+   * \brief L2 norm of this vector on the device.
+   * \note Explicit GPU helper for solver-side reductions.
+   * \return L2 norm result.
    */
-  void GPUSetVal(ScalarType val, bool trigger = true) const;
+  ScalarType GPUNorm() const;
 
   /*!
    * \brief return device pointer that points to the CSysVector values in GPU memory
    */
   inline ScalarType* GetDevicePointer() const { return d_vec_val; }
 
+  /*!
+   * \brief Return the pointer used by expression-template views.
+   */
+  inline const ScalarType* GetExpressionPointer() const {
+    return VecExpr::DeviceExpressionsEnabled() ? d_vec_val : vec_val;
+  }
+
   /*!
    * \brief return the number of local elements in the CSysVector
    */
@@ -301,9 +522,19 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] other - Another vector.
    */
   CSysVector& operator=(const CSysVector& other) {
+#ifdef HAVE_CUDA
+    if (VecExpr::DeviceExpressionsEnabled()) {
+      other.EnsureDeviceData();
+      VecExpr::CVectorView<ScalarType> view(other);
+      VecExpr::AssignDeviceExpression<VecExpr::DeviceAssignOp::Assign>(d_vec_val, nElm, view);
+      MarkDeviceDataModified();
+      return *this;
+    }
+#endif
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i];
     END_CSYSVEC_PARFOR
+    MarkHostDataModified();
     return *this;
   }
 
@@ -311,25 +542,35 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \brief Compound assignement operations with scalars and expressions.
    * \param[in] val/expr - Scalar value or expression.
    */
-#define MAKE_COMPOUND(OP)                                                 \
-  CSysVector& operator OP(ScalarType val) {                               \
-    CSYSVEC_PARFOR                                                        \
-    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val;                  \
-    END_CSYSVEC_PARFOR                                                    \
-    return *this;                                                         \
-  }                                                                       \
-  template <class T>                                                      \
-  CSysVector& operator OP(const VecExpr::CVecExpr<T, ScalarType>& expr) { \
-    CSYSVEC_PARFOR                                                        \
-    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i];    \
-    END_CSYSVEC_PARFOR                                                    \
-    return *this;                                                         \
+#define MAKE_COMPOUND(OP, ASSIGN_OP)                                              \
+  CSysVector& operator OP(ScalarType val) {                                       \
+    if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice<ASSIGN_OP>(val); \
+    CSYSVEC_PARFOR                                                                \
+    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val;                          \
+    END_CSYSVEC_PARFOR                                                            \
+    MarkHostDataModified();                                                       \
+    return *this;                                                                 \
+  }                                                                               \
+  template <class T>                                                              \
+  CSysVector& operator OP(const VecExpr::CVecExpr<T, ScalarType>& expr) {         \
+    if (VecExpr::DeviceExpressionsEnabled()) {                                    \
+      using DeviceExpr = std::remove_cv_t<VecExpr::remove_reference_t<T>>;        \
+      if constexpr (VecExpr::is_device_assignable<DeviceExpr>::value) {           \
+        return AssignDevice<ASSIGN_OP>(expr);                                     \
+      }                                                                           \
+      VecExpr::FlushDeviceExpressionContext();                                    \
+    }                                                                             \
+    CSYSVEC_PARFOR                                                                \
+    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i];            \
+    END_CSYSVEC_PARFOR                                                            \
+    MarkHostDataModified();                                                       \
+    return *this;                                                                 \
   }
-  MAKE_COMPOUND(=)
-  MAKE_COMPOUND(+=)
-  MAKE_COMPOUND(-=)
-  MAKE_COMPOUND(*=)
-  MAKE_COMPOUND(/=)
+  MAKE_COMPOUND(=, VecExpr::DeviceAssignOp::Assign)
+  MAKE_COMPOUND(+=, VecExpr::DeviceAssignOp::Add)
+  MAKE_COMPOUND(-=, VecExpr::DeviceAssignOp::Subtract)
+  MAKE_COMPOUND(*=, VecExpr::DeviceAssignOp::Multiply)
+  MAKE_COMPOUND(/=, VecExpr::DeviceAssignOp::Divide)
 #undef MAKE_COMPOUND
 
   /*!
@@ -344,6 +585,18 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   template <class T>
   ScalarType dot(const VecExpr::CVecExpr<T, ScalarType>& expr) const {
+#ifdef HAVE_CUDA
+    if (VecExpr::DeviceExpressionsEnabled()) {
+      using DeviceExpr = std::remove_cv_t<VecExpr::remove_reference_t<T>>;
+      if constexpr (std::is_same_v<DeviceExpr, CSysVector>) {
+        EnsureDeviceData();
+        expr.derived().EnsureDeviceData();
+        return GPUDot(expr.derived());
+      }
+      VecExpr::FlushDeviceExpressionContext();
+    }
+#endif
+
     /*--- All threads get the same "view" of the vectors. ---*/
     SU2_OMP_BARRIER
 
@@ -502,5 +755,16 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   }
 };
 
+namespace VecExpr {
+
+template <class Scalar>
+CVectorView<Scalar>::CVectorView(const CSysVector<Scalar>& vector) {
+  vector.EnsureDeviceData();
+  data = vector.GetExpressionPointer();
+}
+
+}  // namespace VecExpr
+
 #undef CSYSVEC_PARFOR
 #undef END_CSYSVEC_PARFOR
+#undef SU2_CUDA_HOST_DEVICE
diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh
index 13854381877..eb727477f21 100644
--- a/Common/include/linear_algebra/GPUComms.cuh
+++ b/Common/include/linear_algebra/GPUComms.cuh
@@ -25,6 +25,11 @@
 * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
 */
 
+#pragma once
+
+#ifndef SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH
+#define SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH
+
 #include<cuda_runtime.h>
 #include<iostream>
 
@@ -51,3 +56,5 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t
 }
 
 #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
+
+#endif  // SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH
diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp
index a0d0ce28901..f1f2639e22a 100644
--- a/Common/include/linear_algebra/vector_expressions.hpp
+++ b/Common/include/linear_algebra/vector_expressions.hpp
@@ -39,6 +39,12 @@ namespace VecExpr {
 /// \addtogroup VecExpr
 /// @{
 
+#ifdef __CUDACC__
+#define SU2_CUDA_HOST_DEVICE __host__ __device__
+#else
+#define SU2_CUDA_HOST_DEVICE
+#endif
+
 /*!
  * \brief Base vector expression class.
  * \ingroup BLAS
@@ -59,7 +65,7 @@ class CVecExpr {
   /*!
    * \brief Cast the expression to Derived, usually to allow evaluation via operator[].
    */
-  FORCEINLINE const Derived& derived() const { return static_cast<const Derived&>(*this); }
+  SU2_CUDA_HOST_DEVICE FORCEINLINE const Derived& derived() const { return static_cast<const Derived&>(*this); }
 
   // Allowed from C++14, allows nested expression propagation without
   // manually calling derived() on the expression being evaluated.
@@ -76,8 +82,8 @@ class Bcast : public CVecExpr<Bcast<Scalar>, Scalar> {
 
  public:
   static constexpr bool StoreAsRef = false;
-  FORCEINLINE Bcast(const Scalar& x_) : x(x_) {}
-  FORCEINLINE const Scalar& operator[](size_t) const { return x; }
+  SU2_CUDA_HOST_DEVICE FORCEINLINE Bcast(const Scalar& x_) : x(x_) {}
+  SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t) const { return x; }
 };
 
 /*!
@@ -101,7 +107,14 @@ struct add_lref_if<T, true> {
   using type = remove_reference_t<T>&;
 };
 template <class T>
-using store_t = typename add_lref_if<T, T::StoreAsRef>::type;
+struct store_type {
+  using type = typename add_lref_if<T, T::StoreAsRef>::type;
+};
+template <class T>
+using store_t = typename store_type<T>::type;
+
+template <class T>
+struct is_device_assignable : std::false_type {};
 
 /*--- Namespace from which the math function implementations come. ---*/
 
@@ -120,19 +133,19 @@ namespace math = ::std;
 /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary
  * functions, based on their coefficient-wise implementation (IMPL). ---*/
 
-#define MAKE_UNARY_FUN(FUN, EXPR, IMPL)                             \
-  /*!--- Expression class. ---*/                                    \
-  template <class U, class Scalar>                                  \
-  class EXPR : public CVecExpr<EXPR<U, Scalar>, Scalar> {           \
-    store_t<const U> u;                                             \
-                                                                    \
-   public:                                                          \
-    static constexpr bool StoreAsRef = false;                       \
-    FORCEINLINE EXPR(const U& u_) : u(u_) {}                        \
-    FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \
-  };                                                                \
-  /*!--- Function overload, returns an expression object. ---*/     \
-  template <class U, class S>                                       \
+#define MAKE_UNARY_FUN(FUN, EXPR, IMPL)                                                  \
+  /*!--- Expression class. ---*/                                                         \
+  template <class U, class Scalar>                                                       \
+  class EXPR : public CVecExpr<EXPR<U, Scalar>, Scalar> {                                \
+    store_t<const U> u;                                                                  \
+                                                                                         \
+   public:                                                                               \
+    static constexpr bool StoreAsRef = false;                                            \
+    FORCEINLINE EXPR(const U& u_) : u(u_) {}                                             \
+    SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \
+  };                                                                                     \
+  /*!--- Function overload, returns an expression object. ---*/                          \
+  template <class U, class S>                                                            \
   FORCEINLINE auto FUN(const CVecExpr<U, S>& u) RETURNS(EXPR<U, S>(u.derived()))
 
 #define sign_impl(x) Scalar(1 - 2 * (x < 0))
@@ -158,7 +171,7 @@ MAKE_UNARY_FUN(sign, sign_, sign_impl)
    public:                                                                                                            \
     static constexpr bool StoreAsRef = false;                                                                         \
     FORCEINLINE EXPR(const U& u_, const V& v_) : u(u_), v(v_) {}                                                      \
-    FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i]))                                             \
+    SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i]))                       \
   };                                                                                                                  \
   /*!--- Vector with vector function overload. ---*/                                                                  \
   template <class U, class V, class S>                                                                                \
@@ -241,4 +254,5 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl)
 #undef MAKE_BINARY_FUN
 
 /// @}
+#undef SU2_CUDA_HOST_DEVICE
 }  // namespace VecExpr
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 4c521c29239..65ba6dfa1ab 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -2517,6 +2517,7 @@ enum ENUM_LINEAR_SOLVER_PREC {
   PASTIX_ILU=10,  /*!< \brief PaStiX ILU(k) preconditioner. */
   PASTIX_LU_P,    /*!< \brief PaStiX LU as preconditioner. */
   PASTIX_LDLT_P,  /*!< \brief PaStiX LDLT as preconditioner. */
+  NO_PRECONDITIONER=20, /*!< \brief No preconditioner. */
 };
 static const MapType<std::string, ENUM_LINEAR_SOLVER_PREC> Linear_Solver_Prec_Map = {
   MakePair("JACOBI", JACOBI)
@@ -2526,6 +2527,7 @@ static const MapType<std::string, ENUM_LINEAR_SOLVER_PREC> Linear_Solver_Prec_Ma
   MakePair("PASTIX_ILU", PASTIX_ILU)
   MakePair("PASTIX_LU", PASTIX_LU_P)
   MakePair("PASTIX_LDLT", PASTIX_LDLT_P)
+  MakePair("NONE", NO_PRECONDITIONER)
 };
 
 /*!
diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
index 00ef51c2023..fafa093734e 100644
--- a/Common/src/linear_algebra/CSysMatrix.cpp
+++ b/Common/src/linear_algebra/CSysMatrix.cpp
@@ -54,6 +54,7 @@ CSysMatrix<ScalarType>::CSysMatrix() : rank(SU2_MPI::GetRank()), size(SU2_MPI::G
   col_ind_ilu = nullptr;
 
   invM = nullptr;
+  d_invM = nullptr;
 
 #ifdef USE_MKL
   MatrixMatrixProductJitter = nullptr;
@@ -67,6 +68,8 @@ template <class ScalarType>
 CSysMatrix<ScalarType>::~CSysMatrix() {
   SU2_ZONE_SCOPED
 
+  ReleaseCudaSpMVResources(spmv_resources);
+
   delete[] omp_partitions;
   MemoryAllocation::aligned_free(ILU_matrix);
   MemoryAllocation::aligned_free(matrix);
@@ -76,6 +79,7 @@ CSysMatrix<ScalarType>::~CSysMatrix() {
     GPUMemoryAllocation::gpu_free(d_matrix);
     GPUMemoryAllocation::gpu_free(d_row_ptr);
     GPUMemoryAllocation::gpu_free(d_col_ind);
+    GPUMemoryAllocation::gpu_free(d_invM);
   }
 
 #ifdef USE_MKL
@@ -86,6 +90,10 @@ CSysMatrix<ScalarType>::~CSysMatrix() {
 #endif
 }
 
+#ifndef HAVE_CUDA
+void ReleaseCudaSpMVResources(CudaSpMVResources*& resources) { resources = nullptr; }
+#endif
+
 template <class ScalarType>
 void CSysMatrix<ScalarType>::Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar,
                                         unsigned short neqn, bool EdgeConnect, CGeometry* geometry,
@@ -196,6 +204,10 @@ void CSysMatrix<ScalarType>::Initialize(unsigned long npoint, unsigned long npoi
 
   if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn);
 
+  if (useCuda && diag_needed) {
+    d_invM = GPUMemoryAllocation::gpu_alloc<ScalarType, true>(nPointDomain * nVar * nEqn * sizeof(ScalarType));
+  }
+
   /*--- Thread parallel initialization. ---*/
 
   int num_threads = omp_get_max_threads();
@@ -672,6 +684,19 @@ void CSysMatrix<ScalarType>::MatrixVectorProduct(const CSysVector<ScalarType>& v
 template <class ScalarType>
 void CSysMatrix<ScalarType>::BuildJacobiPreconditioner() {
   SU2_ZONE_SCOPED
+
+  if (useCuda) {
+#ifdef HAVE_CUDA
+    BuildJacobiPreconditionerGPU();
+    return;
+#else
+    SU2_MPI::Error(
+        "\nError in building Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options "
+        "enabled in Meson to access GPU Functions",
+        CURRENT_FUNCTION);
+#endif
+  }
+
   /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/
   SU2_OMP_FOR_DYN(omp_heavy_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++)
@@ -684,6 +709,19 @@ void CSysMatrix<ScalarType>::ComputeJacobiPreconditioner(const CSysVector<Scalar
                                                          CSysVector<ScalarType>& prod, CGeometry* geometry,
                                                          const CConfig* config) const {
   SU2_ZONE_SCOPED
+
+  if (config->GetCUDA()) {
+#ifdef HAVE_CUDA
+    ComputeJacobiPreconditionerGPU(vec, prod, geometry, config);
+    return;
+#else
+    SU2_MPI::Error(
+        "\nError in applying Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options "
+        "enabled in Meson to access GPU Functions",
+        CURRENT_FUNCTION);
+#endif
+  }
+
   /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known. ---*/
   SU2_OMP_BARRIER
   SU2_OMP_FOR_DYN(omp_heavy_size)
diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu
index a3f1a77ca40..35e3c644496 100644
--- a/Common/src/linear_algebra/CSysMatrixGPU.cu
+++ b/Common/src/linear_algebra/CSysMatrixGPU.cu
@@ -55,6 +55,38 @@ inline cusparseIndexType_t GetCusparseIndexType() {
   }
 }
 
+struct CudaSpMVResources {
+  cusparseHandle_t handle = nullptr;
+  cusparseConstSpMatDescr_t mat = nullptr;
+  size_t buffer_size = 0;
+  void* buffer = nullptr;
+};
+
+void ReleaseCudaSpMVResources(CudaSpMVResources*& resources) {
+  if (resources == nullptr) {
+    return;
+  }
+
+  if (resources->buffer != nullptr) {
+    gpuErrChk(cudaFree(resources->buffer));
+    resources->buffer = nullptr;
+    resources->buffer_size = 0;
+  }
+
+  if (resources->mat != nullptr) {
+    cusparseErrChk(cusparseDestroySpMat(resources->mat));
+    resources->mat = nullptr;
+  }
+
+  if (resources->handle != nullptr) {
+    cusparseErrChk(cusparseDestroy(resources->handle));
+    resources->handle = nullptr;
+  }
+
+  delete resources;
+  resources = nullptr;
+}
+
 template <class ScalarType>
 constexpr cudaDataType GetCudaDataType() {
   if constexpr (std::is_same<ScalarType, float>::value) {
@@ -76,6 +108,8 @@ void CSysMatrix<ScalarType>::HtDTransfer(bool trigger) const
 template <class ScalarType>
 void CSysMatrix<ScalarType>::GPUMatrixVectorProduct(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod,
                                                     CGeometry* geometry, const CConfig* config) const {
+  vec.EnsureDeviceData();
+
   if (nVar != nEqn) {
     SU2_MPI::Error("CUDA CSysMatrix matvec with cuSPARSE BSR requires square blocks.", CURRENT_FUNCTION);
   }
@@ -83,8 +117,6 @@ void CSysMatrix<ScalarType>::GPUMatrixVectorProduct(const CSysVector<ScalarType>
   ScalarType* d_vec = vec.GetDevicePointer();
   ScalarType* d_prod = prod.GetDevicePointer();
 
-  vec.HtDTransfer();
-
   const auto indexType = GetCusparseIndexType();
   const auto valueType = GetCudaDataType<ScalarType>();
 
@@ -100,42 +132,58 @@ void CSysMatrix<ScalarType>::GPUMatrixVectorProduct(const CSysVector<ScalarType>
   const ScalarType alpha = 1.0;
   const ScalarType beta = 0.0;
 
-  cusparseHandle_t handle = nullptr;
-  cusparseConstSpMatDescr_t matA = nullptr;
-  cusparseDnVecDescr_t vecX = nullptr;
-  cusparseDnVecDescr_t vecY = nullptr;
+  if (spmv_resources == nullptr) {
+    spmv_resources = new CudaSpMVResources;
 
-  cusparseErrChk(cusparseCreate(&handle));
+    cusparseErrChk(cusparseCreate(&spmv_resources->handle));
 
-  cusparseErrChk(cusparseCreateConstBsr(&matA, brows, bcols, bnnz, blockSize, blockSize, d_row_ptr, d_col_ind, d_matrix,
-                                        indexType, indexType, CUSPARSE_INDEX_BASE_ZERO, valueType, CUSPARSE_ORDER_ROW));
+    cusparseErrChk(cusparseCreateConstBsr(&spmv_resources->mat, brows, bcols, bnnz, blockSize, blockSize, d_row_ptr,
+                                          d_col_ind, d_matrix, indexType, indexType, CUSPARSE_INDEX_BASE_ZERO,
+                                          valueType, CUSPARSE_ORDER_ROW));
+  }
+
+  cusparseDnVecDescr_t vecX = nullptr;
+  cusparseDnVecDescr_t vecY = nullptr;
 
   cusparseErrChk(cusparseCreateDnVec(&vecX, xSize, d_vec, valueType));
   cusparseErrChk(cusparseCreateDnVec(&vecY, ySize, d_prod, valueType));
 
-  size_t bufferSize = 0;
-  void* dBuffer = nullptr;
+  size_t required_buffer_size = 0;
 
-  cusparseErrChk(cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY,
-                                         valueType, CUSPARSE_SPMV_BSR_ALG1, &bufferSize));
+  cusparseErrChk(cusparseSpMV_bufferSize(spmv_resources->handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
+                                         spmv_resources->mat, vecX, &beta, vecY, valueType, CUSPARSE_SPMV_BSR_ALG1,
+                                         &required_buffer_size));
 
-  if (bufferSize > 0) {
-    gpuErrChk(cudaMalloc(&dBuffer, bufferSize));
-  }
+  if (required_buffer_size > spmv_resources->buffer_size) {
+    if (spmv_resources->buffer != nullptr) {
+      gpuErrChk(cudaFree(spmv_resources->buffer));
+    }
 
-  cusparseErrChk(cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, valueType,
-                              CUSPARSE_SPMV_BSR_ALG1, dBuffer));
+    if (required_buffer_size > 0) {
+      gpuErrChk(cudaMalloc(&spmv_resources->buffer, required_buffer_size));
+    } else {
+      spmv_resources->buffer = nullptr;
+    }
 
-  if (dBuffer != nullptr) {
-    gpuErrChk(cudaFree(dBuffer));
+    spmv_resources->buffer_size = required_buffer_size;
   }
 
+  cusparseErrChk(cusparseSpMV(spmv_resources->handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmv_resources->mat,
+                              vecX, &beta, vecY, valueType, CUSPARSE_SPMV_BSR_ALG1, spmv_resources->buffer));
+
   cusparseErrChk(cusparseDestroyDnVec(vecY));
   cusparseErrChk(cusparseDestroyDnVec(vecX));
-  cusparseErrChk(cusparseDestroySpMat(matA));
-  cusparseErrChk(cusparseDestroy(handle));
 
-  prod.DtHTransfer();
+  prod.MarkDeviceDataModified();
 }
-
-template class CSysMatrix<su2mixedfloat>; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.
+template void CSysMatrix<su2mixedfloat>::HtDTransfer(bool trigger) const;
+template void CSysMatrix<su2mixedfloat>::GPUMatrixVectorProduct(const CSysVector<su2mixedfloat>& vec,
+                                                                CSysVector<su2mixedfloat>& prod, CGeometry* geometry,
+                                                                const CConfig* config) const;
+
+#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION)
+template void CSysMatrix<passivedouble>::HtDTransfer(bool trigger) const;
+template void CSysMatrix<passivedouble>::GPUMatrixVectorProduct(const CSysVector<passivedouble>& vec,
+                                                                CSysVector<passivedouble>& prod, CGeometry* geometry,
+                                                                const CConfig* config) const;
+#endif
diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu
new file mode 100644
index 00000000000..4241a0fd538
--- /dev/null
+++ b/Common/src/linear_algebra/CSysPreconditionerGPU.cu
@@ -0,0 +1,113 @@
+/*!
+ * \file CSysPreconditionerGPU.cu
+ * \brief CUDA/GPU skeleton implementations for matrix-based preconditioners.
+ * \author Jesse Li
+ * \version 8.5.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "../../include/linear_algebra/CSysMatrix.inl"
+#include "../../include/linear_algebra/GPUComms.cuh"
+
+namespace {
+
+template <class ScalarType>
+__global__ void ApplyJacobiPreconditionerKernel(const ScalarType* invM, const ScalarType* vec, ScalarType* prod,
+                                                unsigned long nPointDomain, unsigned long nVar) {
+  const auto iPoint = static_cast<unsigned long>(blockIdx.x) * blockDim.x + threadIdx.x;
+  if (iPoint >= nPointDomain) return;
+
+  const auto block = &invM[iPoint * nVar * nVar];
+  const auto rhs = &vec[iPoint * nVar];
+  auto out = &prod[iPoint * nVar];
+
+  for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+    ScalarType sum = ScalarType(0);
+    for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+      sum += block[iVar * nVar + jVar] * rhs[jVar];
+    }
+    out[iVar] = sum;
+  }
+}
+
+}  // namespace
+
+template <class ScalarType>
+void CSysMatrix<ScalarType>::BuildJacobiPreconditionerGPU() {
+  SU2_ZONE_SCOPED
+
+  if (nVar != nEqn) {
+    SU2_MPI::Error("CUDA Jacobi preconditioner requires square blocks.", CURRENT_FUNCTION);
+  }
+
+  if (invM == nullptr) {
+    SU2_MPI::Error("CUDA Jacobi preconditioner was requested without host inverse block storage.", CURRENT_FUNCTION);
+  }
+
+  if (d_invM == nullptr) {
+    d_invM = GPUMemoryAllocation::gpu_alloc<ScalarType, true>(nPointDomain * nVar * nEqn * sizeof(ScalarType));
+  }
+
+  for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
+    InverseDiagonalBlock(iPoint, &(invM[iPoint * nVar * nVar]));
+  }
+
+  gpuErrChk(cudaMemcpy(d_invM, invM, nPointDomain * nVar * nVar * sizeof(ScalarType), cudaMemcpyHostToDevice));
+}
+
+template <class ScalarType>
+void CSysMatrix<ScalarType>::ComputeJacobiPreconditionerGPU(const CSysVector<ScalarType>& vec,
+                                                            CSysVector<ScalarType>& prod, CGeometry* geometry,
+                                                            const CConfig* config) const {
+  (void)geometry;
+  (void)config;
+
+  SU2_ZONE_SCOPED
+
+  vec.EnsureDeviceData();
+
+  if (d_invM == nullptr) {
+    SU2_MPI::Error("CUDA Jacobi preconditioner used before BuildJacobiPreconditionerGPU.", CURRENT_FUNCTION);
+  }
+
+  constexpr unsigned threadsPerBlock = 128;
+  const auto blocks = static_cast<unsigned>((nPointDomain + threadsPerBlock - 1) / threadsPerBlock);
+  ApplyJacobiPreconditionerKernel<<<blocks, threadsPerBlock>>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(),
+                                                               nPointDomain, nVar);
+  gpuErrChk(cudaPeekAtLastError());
+
+  prod.MarkDeviceDataModified();
+}
+
+template void CSysMatrix<su2mixedfloat>::BuildJacobiPreconditionerGPU();
+template void CSysMatrix<su2mixedfloat>::ComputeJacobiPreconditionerGPU(const CSysVector<su2mixedfloat>& vec,
+                                                                        CSysVector<su2mixedfloat>& prod,
+                                                                        CGeometry* geometry,
+                                                                        const CConfig* config) const;
+
+#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION)
+template void CSysMatrix<passivedouble>::BuildJacobiPreconditionerGPU();
+template void CSysMatrix<passivedouble>::ComputeJacobiPreconditionerGPU(const CSysVector<passivedouble>& vec,
+                                                                        CSysVector<passivedouble>& prod,
+                                                                        CGeometry* geometry,
+                                                                        const CConfig* config) const;
+#endif
diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
index cd100f311ab..f6deffea4df 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -41,7 +41,6 @@ SU2_RESTORE_WARNING
 #include <complex>
 #include <iostream>
 #include <limits>
-#include <type_traits>
 #include <memory>
 
 namespace {
@@ -423,11 +422,9 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
 
   const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   const bool flexible = !precond.IsIdentity();
-  /*--- If we call the solver outside of a parallel region, but the number of threads allows,
-   * we still want to parallelize some of the expensive operations. ---*/
   const bool nestedParallel = !omp_in_parallel() && omp_get_max_threads() > 1;
 
-  /*---  Check the subspace size ---*/
+  /*--- Check the subspace size. ---*/
 
   if (m < 1) {
     SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION);
@@ -437,7 +434,7 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
     SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION);
   }
 
-  /*--- Allocate if not allocated yet ---*/
+  /*--- Allocate if not allocated yet. ---*/
 
   if (V.size() <= m || (flexible && Z.size() <= m)) {
     BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
@@ -460,6 +457,7 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   sn = ScalarType(0);
   cs = ScalarType(0);
   y = ScalarType(0);
+
   su2matrix<ScalarType> H(m + 1, m);
   H = ScalarType(0);
 
@@ -485,12 +483,12 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   ResetDeflation();
   END_SU2_OMP_MASTER
 
-  /*--- Set the norm to the initial initial residual value ---*/
+  /*--- Set the norm to the initial initial residual value. ---*/
 
   if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta;
 
   if (beta < tol * norm0 || beta < eps) {
-    /*--- System is already solved ---*/
+    /*--- System is already solved. ---*/
 
     if (masterRank) {
       SU2_OMP_MASTER
@@ -506,11 +504,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
 
   V[0] /= -beta;
 
-  /*--- Initialize the RHS of the reduced system ---*/
+  /*--- Initialize the RHS of the reduced system. ---*/
 
   g[0] = beta;
 
-  /*--- Output header information including initial residual ---*/
+  /*--- Output header information including initial residual. ---*/
 
   unsigned long i = 0;
   if (monitoring && masterRank) {
@@ -521,30 +519,30 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
     END_SU2_OMP_MASTER
   }
 
-  /*---  Loop over all search directions ---*/
+  /*--- Loop over all search directions. ---*/
 
   for (i = 0; i < m; i++) {
-    /*---  Check if solution has converged ---*/
+    /*--- Check if solution has converged. ---*/
 
     if (beta < tol * norm0) break;
 
     if (flexible) {
-      /*---  Precondition the CSysVector v[i] and store result in z[i] ---*/
+      /*--- Precondition the CSysVector v[i] and store result in z[i]. ---*/
 
       precond(V[i], Z[i]);
 
-      /*---  Add to Krylov subspace ---*/
+      /*--- Add to Krylov subspace. ---*/
 
       mat_vec(Z[i], V[i + 1]);
     } else {
       mat_vec(V[i], V[i + 1]);
     }
 
-    /*---  Modified Gram-Schmidt orthogonalization ---*/
+    /*--- Modified Gram-Schmidt orthogonalization. ---*/
 
     bool orthog_ok;
     if (nestedParallel) {
-      /*--- "omp parallel if" does not work well here ---*/
+      /*--- "omp parallel if" does not work well here. ---*/
       SU2_OMP_PARALLEL
       orthog_ok = ModGramSchmidt(true, i, H, V);
       END_SU2_OMP_PARALLEL
@@ -561,18 +559,18 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
       break;
     }
 
-    /*---  Apply old Givens rotations to new column of the Hessenberg matrix then generate the
-     new Givens rotation matrix and apply it to the last two elements of H[:][i] and g ---*/
+    /*--- Apply old Givens rotations to new column of the Hessenberg matrix then generate the
+     new Givens rotation matrix and apply it to the last two elements of H[:][i] and g. ---*/
 
     for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i));
     GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]);
     ApplyGivens(sn[i], cs[i], g[i], g[i + 1]);
 
-    /*---  Set L2 norm of residual and check if solution has converged ---*/
+    /*--- Set L2 norm of residual and check if solution has converged. ---*/
 
     beta = fabs(g[i + 1]);
 
-    /*---  Output the relative residual if necessary ---*/
+    /*--- Output the relative residual if necessary. ---*/
 
     if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) {
       SU2_OMP_MASTER
@@ -581,15 +579,14 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
     }
   }
 
-  /*---  Solve the least-squares system and update solution ---*/
+  /*--- Solve the least-squares system and update solution. ---*/
 
   SolveReduced(i, H, g, y);
 
   const auto& basis = flexible ? Z : V;
-
   LinearCombination(nestedParallel, i, basis, y, x, true);
 
-  /*---  Recalculate final (neg.) residual (this should be optional) ---*/
+  /*--- Recalculate final (neg.) residual (this should be optional). ---*/
 
   if (monitoring && config->GetComm_Level() == COMM_FULL) {
     if (masterRank) {
@@ -1464,84 +1461,86 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
 
     HandleTemporariesIn(LinSysRes, LinSysSol);
 
-    auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
-
-    /*--- Build preconditioner. ---*/
-
-    const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
-    auto* normal_prec = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
-    normal_prec->Build();
-
-    CPreconditioner<ScalarType>* nested_prec = nullptr;
-    if (nested) {
-      auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
-        /*--- Initialize to 0 to be safe. ---*/
-        v = ScalarType{};
-        ScalarType res{};
-        /*--- Handle other types here if desired but do not call Solve because
-         * that will create issues with the AD external function. ---*/
-        if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
-          inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
-        } else {
-          const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
-          inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
-        }
-      };
-      nested_prec = new CAbstractPreconditioner<ScalarType>(f);
-    }
-    const auto* precond = nested ? nested_prec : normal_prec;
+    {
+      auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
+
+      /*--- Build preconditioner. ---*/
+
+      const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+      auto* normal_prec = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
+      normal_prec->Build();
+
+      CPreconditioner<ScalarType>* nested_prec = nullptr;
+      if (nested) {
+        auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
+          /*--- Initialize to 0 to be safe. ---*/
+          v = ScalarType{};
+          ScalarType res{};
+          /*--- Handle other types here if desired but do not call Solve because
+           * that will create issues with the AD external function. ---*/
+          if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
+            inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
+          } else {
+            const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
+            inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
+          }
+        };
+        nested_prec = new CAbstractPreconditioner<ScalarType>(f);
+      }
+      const auto* precond = nested ? nested_prec : normal_prec;
 
-    /*--- Solve system. ---*/
+      /*--- Solve system. ---*/
 
-    ScalarType residual = 0.0;
+      ScalarType residual = 0.0;
 
-    switch (KindSolver) {
-      case BCGSTAB:
-        IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
-                                       ScreenOutput, config);
-        break;
-      case FGMRES:
-        IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
-                                      ScreenOutput, config);
-        break;
-      case FGCRODR:
-        IterLinSol = FGCRODR_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
-                                       ScreenOutput, config);
-        break;
-      case RESTARTED_FGMRES:
-        IterLinSol = RFGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
-                                       ScreenOutput, config);
-        break;
-      case CONJUGATE_GRADIENT:
-        IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
-                                  ScreenOutput, config);
-        break;
-      case SMOOTHER:
-        IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
+      switch (KindSolver) {
+        case BCGSTAB:
+          IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
+                                         residual, ScreenOutput, config);
+          break;
+        case FGMRES:
+          IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
                                         ScreenOutput, config);
-        break;
-      case PASTIX_LDLT:
-      case PASTIX_LU:
-        Jacobian.BuildPastixPreconditioner(geometry, config, KindSolver);
-        Jacobian.ComputePastixPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config);
-        IterLinSol = 1;
-        residual = 1e-20;
-        break;
-      default:
-        SU2_MPI::Error("Unknown type of linear solver.", CURRENT_FUNCTION);
-    }
+          break;
+        case FGCRODR:
+          IterLinSol = FGCRODR_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
+                                         residual, ScreenOutput, config);
+          break;
+        case RESTARTED_FGMRES:
+          IterLinSol = RFGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
+                                         residual, ScreenOutput, config);
+          break;
+        case CONJUGATE_GRADIENT:
+          IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
+                                    ScreenOutput, config);
+          break;
+        case SMOOTHER:
+          IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
+                                          residual, ScreenOutput, config);
+          break;
+        case PASTIX_LDLT:
+        case PASTIX_LU:
+          Jacobian.BuildPastixPreconditioner(geometry, config, KindSolver);
+          Jacobian.ComputePastixPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config);
+          IterLinSol = 1;
+          residual = 1e-20;
+          break;
+        default:
+          SU2_MPI::Error("Unknown type of linear solver.", CURRENT_FUNCTION);
+      }
 
-    SU2_OMP_MASTER {
-      Residual = residual;
-      Iterations = IterLinSol;
+      SU2_OMP_MASTER {
+        Residual = residual;
+        Iterations = IterLinSol;
+      }
+      END_SU2_OMP_MASTER
+
+      delete normal_prec;
+      delete nested_prec;
     }
-    END_SU2_OMP_MASTER
 
     HandleTemporariesOut(LinSysSol);
 
-    delete normal_prec;
-    delete nested_prec;
-
     if (TapeActive) {
       /*--- To keep the behavior of SU2_DOT, but not strictly required since jacobian is symmetric(?). ---*/
       const bool RequiresTranspose =
diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp
index f7df34633a0..c5c4dc9a3cb 100644
--- a/Common/src/linear_algebra/CSysVector.cpp
+++ b/Common/src/linear_algebra/CSysVector.cpp
@@ -65,6 +65,10 @@ void CSysVector<ScalarType>::Initialize(unsigned long numBlk, unsigned long numB
       for (auto i = 0ul; i < nElm; i++) vec_val[i] = val[i];
     }
   }
+
+  host_data_valid = true;
+  device_data_valid = false;
+  device_context_id = 0;
 }
 
 template <class ScalarType>
@@ -78,6 +82,23 @@ const su2matrix<ScalarType>& CSysVector<ScalarType>::multiDot(const std::vector<
 
   if (n == 0 || m == 0) return shared;
 
+#ifdef HAVE_CUDA
+  if (VecExpr::DeviceExpressionsEnabled()) {
+    SU2_OMP_MASTER {
+      shared.resize(n, m);
+      for (size_t i = 0; i < n; ++i) {
+        for (size_t j = 0; j < m; ++j) {
+          shared(i, j) = V[i0 + i].GPUDot(W[j]);
+        }
+      }
+    }
+    END_SU2_OMP_MASTER
+
+    SU2_OMP_BARRIER
+    return shared;
+  }
+#endif
+
   SU2_OMP_BARRIER
   const size_t size = V[0].nElmDomain;
 
@@ -134,6 +155,7 @@ const su2matrix<ScalarType>& CSysVector<ScalarType>::multiDot(const std::vector<
 
 template <class ScalarType>
 CSysVector<ScalarType>::~CSysVector() {
+  VecExpr::UnregisterDeviceModifiedVector(this);
   if constexpr (!std::is_trivial_v<ScalarType>) {
     for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType();
   }
diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu
index 94ec17bb88f..e6429fcb351 100644
--- a/Common/src/linear_algebra/CSysVectorGPU.cu
+++ b/Common/src/linear_algebra/CSysVectorGPU.cu
@@ -27,23 +27,257 @@
 
 #include "../../include/linear_algebra/CSysVector.hpp"
 #include "../../include/linear_algebra/GPUComms.cuh"
+#include <cublas_v2.h>
+#include <algorithm>
+#include <type_traits>
+#include <vector>
 
-template<class ScalarType>
-void CSysVector<ScalarType>::HtDTransfer(bool trigger) const
-{
-   if(trigger) gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType)*nElm), cudaMemcpyHostToDevice));
+namespace {
+
+thread_local cublasHandle_t active_solver_blas_handle = nullptr;
+thread_local unsigned active_solver_blas_depth = 0;
+
+cublasHandle_t GetBlasHandle(const char* creation_error, bool& owns_handle) {
+  if (active_solver_blas_handle != nullptr) {
+    owns_handle = false;
+    return active_solver_blas_handle;
+  }
+
+  cublasHandle_t handle = nullptr;
+  owns_handle = true;
+  if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS) {
+    SU2_MPI::Error(creation_error, CURRENT_FUNCTION);
+  }
+  return handle;
+}
+
+void ReleaseBlasHandle(cublasHandle_t handle, bool owns_handle, const char* destruction_error) {
+  if (owns_handle && handle != nullptr && cublasDestroy(handle) != CUBLAS_STATUS_SUCCESS) {
+    SU2_MPI::Error(destruction_error, CURRENT_FUNCTION);
+  }
+}
+
+void BeginDeviceBlasContext() {
+  if (active_solver_blas_depth == 0) {
+    cublasHandle_t handle = nullptr;
+    if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS) {
+      SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear algebra context.", CURRENT_FUNCTION);
+    }
+    active_solver_blas_handle = handle;
+  }
+  ++active_solver_blas_depth;
 }
 
-template<class ScalarType>
-void CSysVector<ScalarType>::DtHTransfer(bool trigger) const
-{
-   if(trigger) gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType)*nElm), cudaMemcpyDeviceToHost));
+void EndDeviceBlasContext() {
+  if (active_solver_blas_depth == 0) {
+    SU2_MPI::Error("GPU linear algebra BLAS context ended without a matching begin.", CURRENT_FUNCTION);
+  }
+
+  --active_solver_blas_depth;
+  if (active_solver_blas_depth == 0) {
+    auto status = cublasDestroy(active_solver_blas_handle);
+    active_solver_blas_handle = nullptr;
+    if (status != CUBLAS_STATUS_SUCCESS) {
+      SU2_MPI::Error("cuBLAS handle destruction failed for the GPU linear algebra context.", CURRENT_FUNCTION);
+    }
+  }
+}
+
+}  // namespace
+
+namespace VecExpr {
+
+namespace {
+thread_local bool device_expressions_enabled = false;
+struct DeviceExpressionContextState {
+  bool enabled = false;
+  bool previous = false;
+  unsigned long previous_context_id = 0;
+};
+struct ModifiedVector {
+  const void* vector = nullptr;
+  void (*sync)(const void*) = nullptr;
+};
+thread_local std::vector<DeviceExpressionContextState> device_context_stack;
+thread_local std::vector<ModifiedVector> modified_vectors;
+thread_local unsigned long next_device_context_id = 0;
+thread_local unsigned long current_device_context_id = 0;
+}  // namespace
+
+bool DeviceExpressionsEnabled() { return device_expressions_enabled; }
+
+void BeginDeviceExpressionContext(bool enabled) {
+  device_context_stack.push_back({enabled, device_expressions_enabled, current_device_context_id});
+  if (!enabled) return;
+
+  BeginDeviceBlasContext();
+  device_expressions_enabled = true;
+  current_device_context_id = ++next_device_context_id;
+}
+
+void FlushDeviceExpressionContext() {
+  for (const auto& entry : modified_vectors) entry.sync(entry.vector);
+  modified_vectors.clear();
+}
+
+void EndDeviceExpressionContext() {
+  if (device_context_stack.empty()) {
+    SU2_MPI::Error("Device expression context ended without a matching begin.", CURRENT_FUNCTION);
+  }
+
+  const auto state = device_context_stack.back();
+  device_context_stack.pop_back();
+
+  if (state.enabled) {
+    FlushDeviceExpressionContext();
+    device_expressions_enabled = state.previous;
+    current_device_context_id = state.previous_context_id;
+    EndDeviceBlasContext();
+  }
+}
+
+unsigned long CurrentDeviceExpressionContextId() { return current_device_context_id; }
+
+void RegisterDeviceModifiedVector(const void* vector, void (*sync)(const void*)) {
+  if (!DeviceExpressionsEnabled()) return;
+  const auto exists = std::any_of(modified_vectors.begin(), modified_vectors.end(),
+                                  [vector](const auto& entry) { return entry.vector == vector; });
+  if (!exists) modified_vectors.push_back({vector, sync});
+}
+
+void UnregisterDeviceModifiedVector(const void* vector) {
+  modified_vectors.erase(std::remove_if(modified_vectors.begin(), modified_vectors.end(),
+                                        [vector](const auto& entry) { return entry.vector == vector; }),
+                         modified_vectors.end());
+}
+
+}  // namespace VecExpr
+
+template <class ScalarType>
+void CSysVector<ScalarType>::HtDTransfer(bool trigger) const {
+  if (trigger)
+    gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType) * nElm), cudaMemcpyHostToDevice));
 }
 
-template<class ScalarType>
-void CSysVector<ScalarType>::GPUSetVal(ScalarType val, bool trigger) const
-{
-   if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm)));
+template <class ScalarType>
+void CSysVector<ScalarType>::DtHTransfer(bool trigger) const {
+  if (trigger)
+    gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType) * nElm), cudaMemcpyDeviceToHost));
 }
 
-template class CSysVector<su2double>; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.
+template <class ScalarType>
+ScalarType CSysVector<ScalarType>::GPUDot(const CSysVector& other) const {
+  EnsureDeviceData();
+  other.EnsureDeviceData();
+
+  bool owns_handle = false;
+  cublasHandle_t handle = GetBlasHandle("cuBLAS handle creation failed in CSysVector::GPUDot.", owns_handle);
+  cublasStatus_t status = CUBLAS_STATUS_SUCCESS;
+
+  ScalarType local_dot = ScalarType(0);
+
+  if constexpr (std::is_same_v<ScalarType, float>) {
+    status = cublasSdot(handle, static_cast<int>(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1,
+                        &local_dot);
+  } else if constexpr (std::is_same_v<ScalarType, double>) {
+    status = cublasDdot(handle, static_cast<int>(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1,
+                        &local_dot);
+  } else {
+    ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot.");
+    SU2_MPI::Error("Unsupported ScalarType in CSysVector::GPUDot.", CURRENT_FUNCTION);
+    return ScalarType(0);
+  }
+
+  if (status != CUBLAS_STATUS_SUCCESS) {
+    ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot.");
+    SU2_MPI::Error("cuBLAS dot failed in CSysVector::GPUDot.", CURRENT_FUNCTION);
+    return ScalarType(0);
+  }
+
+  ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot.");
+
+  ScalarType global_dot = ScalarType(0);
+  const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE;
+  SelectMPIWrapper<ScalarType>::W::Allreduce(&local_dot, &global_dot, 1, mpi_type, MPI_SUM, SU2_MPI::GetComm());
+
+  return global_dot;
+}
+
+template <class ScalarType>
+ScalarType CSysVector<ScalarType>::GPUNorm() const {
+  return sqrt(GPUDot(*this));
+}
+
+namespace {
+
+template <class Scalar>
+using DeviceView = VecExpr::CVectorView<Scalar>;
+
+template <class Scalar>
+using DeviceBcast = VecExpr::Bcast<Scalar>;
+
+template <class Scalar>
+using DeviceNeg = VecExpr::minus_<CSysVector<Scalar>, Scalar>;
+
+template <class Scalar>
+using DeviceScale = VecExpr::mul_<CSysVector<Scalar>, DeviceBcast<Scalar>, Scalar>;
+
+template <class Scalar>
+using DeviceScale2 = VecExpr::add_<DeviceScale<Scalar>, DeviceScale<Scalar>, Scalar>;
+
+template <class Scalar>
+using DeviceScale3 = VecExpr::add_<DeviceScale2<Scalar>, DeviceScale<Scalar>, Scalar>;
+
+template <class Scalar>
+using DeviceScale4 = VecExpr::add_<DeviceScale3<Scalar>, DeviceScale<Scalar>, Scalar>;
+
+}  // namespace
+
+#define INSTANTIATE_DEVICE_ASSIGN(SCALAR, OP, EXPR)                                                           \
+  template void VecExpr::AssignDeviceExpression<VecExpr::DeviceAssignOp::OP, SCALAR, EXPR<SCALAR>>(          \
+      SCALAR*, unsigned long, const VecExpr::CVecExpr<EXPR<SCALAR>, SCALAR>&)
+
+#define INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, EXPR)      \
+  INSTANTIATE_DEVICE_ASSIGN(SCALAR, Assign, EXPR);        \
+  INSTANTIATE_DEVICE_ASSIGN(SCALAR, Add, EXPR);           \
+  INSTANTIATE_DEVICE_ASSIGN(SCALAR, Subtract, EXPR);      \
+  INSTANTIATE_DEVICE_ASSIGN(SCALAR, Multiply, EXPR);      \
+  INSTANTIATE_DEVICE_ASSIGN(SCALAR, Divide, EXPR)
+
+#define INSTANTIATE_DEVICE_ASSIGN_SCALAR(SCALAR)          \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceBcast);    \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceView);     \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceNeg);      \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale);    \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale2);   \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale3);   \
+  INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale4)
+
+INSTANTIATE_DEVICE_ASSIGN_SCALAR(su2mixedfloat);
+
+#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION)
+INSTANTIATE_DEVICE_ASSIGN_SCALAR(passivedouble);
+#endif
+
+#undef INSTANTIATE_DEVICE_ASSIGN_SCALAR
+#undef INSTANTIATE_DEVICE_ASSIGN_EXPR
+#undef INSTANTIATE_DEVICE_ASSIGN
+
+template void CSysVector<su2mixedfloat>::HtDTransfer(bool trigger) const;
+template void CSysVector<su2mixedfloat>::DtHTransfer(bool trigger) const;
+template su2mixedfloat CSysVector<su2mixedfloat>::GPUDot(const CSysVector<su2mixedfloat>& other) const;
+template su2mixedfloat CSysVector<su2mixedfloat>::GPUNorm() const;
+
+#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION)
+template void CSysVector<passivedouble>::HtDTransfer(bool trigger) const;
+template void CSysVector<passivedouble>::DtHTransfer(bool trigger) const;
+template passivedouble CSysVector<passivedouble>::GPUDot(const CSysVector<passivedouble>& other) const;
+template passivedouble CSysVector<passivedouble>::GPUNorm() const;
+#endif
+
+#ifdef CODI_REVERSE_TYPE
+template void CSysVector<su2double>::HtDTransfer(bool trigger) const;
+template void CSysVector<su2double>::DtHTransfer(bool trigger) const;
+template su2double CSysVector<su2double>::GPUDot(const CSysVector<su2double>& other) const;
+template su2double CSysVector<su2double>::GPUNorm() const;
+#endif
diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build
index 7b880b29c1e..0142e24ec4d 100644
--- a/Common/src/linear_algebra/meson.build
+++ b/Common/src/linear_algebra/meson.build
@@ -6,5 +6,5 @@ common_src += files(['CSysSolve_b.cpp',
                      'blas_structure.cpp'])
 
   if get_option('enable-cuda')
-    common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu',])
+    common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu', 'CSysPreconditionerGPU.cu',])
 endif
diff --git a/Common/src/meson.build b/Common/src/meson.build
index f385c4a32ed..710e764303f 100644
--- a/Common/src/meson.build
+++ b/Common/src/meson.build
@@ -6,6 +6,27 @@ common_src =files(['graph_coloring_structure.cpp',
            '../include/parallelization/mpi_structure.cpp',
            '../include/parallelization/omp_structure.cpp'])
 
+common_cuda_cpp_args = []
+foreach arg : su2_cpp_args
+  if arg.startswith('-D') or arg.startswith('-U')
+    common_cuda_cpp_args += arg
+  endif
+endforeach
+
+common_cuda_rev_args = common_cuda_cpp_args
+foreach arg : codi_rev_args
+  if arg.startswith('-D') or arg.startswith('-U')
+    common_cuda_rev_args += arg
+  endif
+endforeach
+
+common_cuda_for_args = common_cuda_cpp_args
+foreach arg : codi_for_args
+  if arg.startswith('-D') or arg.startswith('-U')
+    common_cuda_for_args += arg
+  endif
+endforeach
+
 subdir('linear_algebra')
 subdir('toolboxes')
 subdir('geometry')
@@ -24,7 +45,8 @@ if get_option('enable-normal')
                           common_src,
                           install : false,
                           dependencies : su2_deps,
-                          cpp_args: [default_warning_flags, su2_cpp_args])
+                          cpp_args: [default_warning_flags, su2_cpp_args],
+                          cuda_args: common_cuda_cpp_args)
 
   common_dep = declare_dependency(link_with: common,
                                   include_directories : common_include)
@@ -36,7 +58,8 @@ if get_option('enable-autodiff')
                             common_src,
                             install : false,
                             dependencies : [su2_deps, codi_dep],
-                            cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args])
+                            cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args],
+                            cuda_args: common_cuda_rev_args)
 
   commonAD_dep = declare_dependency(link_with: commonAD,
                                     include_directories : common_include)
@@ -49,7 +72,8 @@ if get_option('enable-directdiff')
                             common_src,
                             install : false,
                             dependencies : [su2_deps, codi_dep],
-                            cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args])
+                            cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args],
+                            cuda_args: common_cuda_for_args)
 
   commonDD_dep = declare_dependency(link_with: commonDD,
                                     include_directories : common_include)
diff --git a/meson.build b/meson.build
index d43e2b031d8..43786c100e9 100644
--- a/meson.build
+++ b/meson.build
@@ -20,13 +20,18 @@ python = pymod.find_installation()
 if get_option('enable-cuda')
   add_languages('cuda')
   add_global_arguments('-arch=sm_86', language : 'cuda')
-  cuda_deps = [meson.get_compiler('cuda').find_library('cusparse', required : true)]
+  cuda_deps = [
+    meson.get_compiler('cuda').find_library('cusparse', required : true),
+    meson.get_compiler('cuda').find_library('cublas', required : true),
+  ]
 else
   cuda_deps = []
 endif
 
 su2_cpp_args = []
 su2_deps     = [declare_dependency(include_directories: 'externals/CLI11')] + cuda_deps
+codi_rev_args = []
+codi_for_args = []
 
 default_warning_flags = []
 if build_machine.system() != 'windows'