Move GMRES from example to sparse experimental (#1620)

* Add GMRES boilerplate * Remove old par_ilut code * Remove gmres_symbolic * Use correct API for GMRES * Take CSR, not views * Better handling of GmresStats * Add KokkosSparse_Preconditioner.hpp * Examples working * Fix for layoutleft ETI issues * Fix intel warning * Add preconditioner test and refactor prec handling in gmres * Put correct author as contact * Move nrows out of the GMRES handle, get it from the CRS arg * Documentation for handle * Move tol to handle constructor arg * Get rid of unnecessary items on GMRES handle * Improve/update documentation * Update licenses for files I made * Rename gmres files without numeric
kokkos · Dec 12, 2022 · 61ff6a3 · 61ff6a3
1 parent 9347873
commit 61ff6a3
Show file tree

Hide file tree

Showing 42 changed files with 1,733 additions and 1,036 deletions.
diff --git a/example/gmres/CMakeLists.txt b/example/gmres/CMakeLists.txt
@@ -6,19 +6,6 @@ KOKKOSKERNELS_ADD_EXECUTABLE(
   SOURCES ex_real_A.cpp
   )
 
-# FIXME_SYCL CUDA_ERROR_INVALID_ADDRESS_SPACE
-IF(NOT KOKKOS_ENABLE_SYCL)
-KOKKOSKERNELS_ADD_EXECUTABLE_AND_TEST(
-  gmres_test_real_A
-  SOURCES test_real_A.cpp
-  )
-ENDIF()
-
-KOKKOSKERNELS_ADD_EXECUTABLE(
-  gmres_test_cmplx_A
-  SOURCES test_cmplx_A.cpp
-  )
-
 # FIXME_SYCL CUDA_ERROR_INVALID_ADDRESS_SPACE
 IF(NOT KOKKOS_ENABLE_SYCL)
 KOKKOSKERNELS_ADD_EXECUTABLE_AND_TEST(

diff --git a/example/gmres/README.md b/example/gmres/README.md
@@ -5,62 +5,52 @@ GMRES (the Generalized Minimum RESidual method) is an iterative solver for spars
 ### Command-Line parameters for ex\_real\_A:
 Current solver parameters for the real-valued example are as follows:
 
-"**--filename**   :  The name of a matrix market (.mtx) file for matrix A (Default bcsstk09.mtx)."    
-"**--max-subsp**   :  The maximum size of the Kyrlov subspace before restarting (Default 50)."   
-"**--max-restarts:**  Maximum number of GMRES restarts (Default 50)."   
-"**--tol        :**  Convergence tolerance.  (Default 1e-10)."   
-"**--ortho       :**  Type of orthogonalization. Use 'CGS2' or 'MGS'. (Default 'CGS2')"    
+"**--filename**   :  The name of a matrix market (.mtx) file for matrix A (Default bcsstk09.mtx)."
+"**--max-subsp**   :  The maximum size of the Kyrlov subspace before restarting (Default 50)."
+"**--max-restarts:**  Maximum number of GMRES restarts (Default 50)."
+"**--tol        :**  Convergence tolerance.  (Default 1e-10)."
+"**--ortho       :**  Type of orthogonalization. Use 'CGS2' or 'MGS'. (Default 'CGS2')"
 "**--rand\_rhs**    :  Generate a random right-hand side b.  (Without this option, the solver default generates b = vector of ones.)"
 
-
-### Solver template paramters:
-The GMRES solver has the following template parameters:   
-**ScalarType:** Type of scalars in the matrix A. (double, float, half, Kokkos::complex<double>, etc.)   
-**Layout:** Layout of vectors X and B, Kokkos::LayoutLeft, or other Kokkos layouts.   
-**EXSP:** A Kokkos execution space.   
-**OrdinalType:** The ordinal type from the Kokkos::CrsMatrix A.
-
-
 ### Solver input parameters:
-The gmres function takes the following input paramters:   
-**A:** A Kokkos::CrsMatrix for the linar system Ax=b.   
-**B:** A Kokkos::View that is the system right-hand side. Must have B.extent(1)=1. (Currently only one right-hand side is supported.)   
-**X:** A Kokkos::View that is used as both the initial vector for the GMRES iteration and the output for the solution vector.  (Must have X.extent(1)=1.)   
-**opts:** A 'GmresOpts' struct as described below.   
-**M:** A pointer to a KokkosSparse::Experimental::Preconditioner. Only right preconditioning is supported at this time.  
-
-The solver has a 'GmresOpts' struct to pass in solver options.  Available options are:   
-**tol:** The convergence tolerance for GMRES.  Based upon the relative residual. The solver will terminate when norm(b-Ax)/norm(b) <= tol. (Default: 1e-8)   
-**m:** The restart length (maximum subspace size) for GMRES.  (Default: 50)   
-**maxRestart:** The maximum number of restarts (or 'cycles') that GMRES is to perform. (Default: 50)   
-**ortho:** The orthogonalization type.  Can be "CGS2" (Default) or "MGS".  (Two iterations of Classical Gram-Schmidt, or one iteration of Modified Gram-Schmidt.)   
+The gmres function takes the following input paramters:
+**A:** A Kokkos::CrsMatrix for the linar system Ax=b.
+**B:** A Kokkos::View that is the system right-hand side. Must have B.extent(1)=1. (Currently only one right-hand side is supported.)
+**X:** A Kokkos::View that is used as both the initial vector for the GMRES iteration and the output for the solution vector.  (Must have X.extent(1)=1.)
+**M:** A pointer to a KokkosSparse::Experimental::Preconditioner. Only right preconditioning is supported at this time.
+
+### Handle input parameters:
+The solver has a GMRESHandle struct to pass in solver options.  Available options are:
+**tol:** The convergence tolerance for GMRES.  Based upon the relative residual. The solver will terminate when norm(b-Ax)/norm(b) <= tol. (Default: 1e-8)
+**m:** The restart length (maximum subspace size) for GMRES.  (Default: 50)
+**maxRestart:** The maximum number of restarts (or 'cycles') that GMRES is to perform. (Default: 50)
+**ortho:** The orthogonalization type.  Can be "CGS2" (Default) or "MGS".  (Two iterations of Classical Gram-Schmidt, or one iteration of Modified Gram-Schmidt.)
+**verbose:** Tells solve to print more information
 
 ### Solver Output:
-The solver outputs a 'GmresStats' struct with solver statistics.  These include:
-
-**numIters**: The number of iterations the GMRES solver took before terminating.   
-**endRelRes**: The ending relative residual norm attained in the GMRES solve.   
-**convFlagVal**: An enum FLAG value which will return "Conv" if the solver converged, "NoConv" if the solver did not converge, or "LOA" if the solver converged with respect to the internally computed residual norm but lost accuracy from the true residual norm.    
-**convFlag**: A string giving the convergence status of the solver upon exit.  
+The GMRESHandle struct is also used to pass back solver statistics. These include:
+**numIters**: The number of iterations the GMRES solver took before terminating.
+**endRelRes**: The ending relative residual norm attained in the GMRES solve.
+**convFlagVal**: An enum FLAG value which will return "Conv" if the solver converged, "NoConv" if the solver did not converge, or "LOA" if the solver converged with respect to the internally computed residual norm but lost accuracy from the true residual norm.
 
 ### Test files:
 
-You can download the matrices for the real example and the complex test from the [SuiteSparse](https://sparse.tamu.edu/) matrix collection.  The two matrices needed are:  
+You can download the matrices for the real example and the complex test from the [SuiteSparse](https://sparse.tamu.edu/) matrix collection.  The two matrices needed are:
 * **young1c** (for complex-valued test)
 * **bcsstk09** (for real-valued example)
 
 The real-valued test uses a matrix generated directly by Kokkos Kernels.
 
 ### Test Measurements:
-These measurements were taken on 7/23/21, running on an NVIDIA V100 GPU on Weaver7.  
+These measurements were taken on 7/23/21, running on an NVIDIA V100 GPU on Weaver7.
 (Timings based upon the GMRES::TotalTime profiling region.)
 
 **ex\_real\_A:** Converges in 2271 iterations and 0.9629 seconds.
 
-(The two following timings total the time for the CGS2 and MGS tests.)   
+(The two following timings total the time for the CGS2 and MGS tests.)
 **test\_real\_A:** Converges in 30 iterations (with a restart size of 15) and 0.2536 seconds.
 
-**test\_cmplx\_A:** Converges in 652 iterations (to a tolerance of 1e-5) in 2.822 seconds.  
+**test\_cmplx\_A:** Converges in 652 iterations (to a tolerance of 1e-5) in 2.822 seconds.
 
 ### Concerns, enhancements, or bug reporting:
 If you wish to suggest an enhancement or make a bug report for this solver code, please post an issue at https://github.com/kokkos/kokkos-kernels/issues or email jloe@sandia.gov.

diff --git a/example/gmres/READMEPreconditioners.md b/example/gmres/READMEPreconditioners.md
@@ -1,44 +1,40 @@
 ## KokkosSparse Preconditioner Interface:
 
-The `KokkosSparse_Preconditioner` class provides an abstract base class interface to use Kokkos-based preconditioners with iterative linear solvers. In particular, this class is designed to work with the Kokkos-based GMRES implementation in `examples/gmres`. It may also be useful for integrating Kokkos-based preconditioners into other solver codes and packages. (For Trilinos users: This class is loosely based upon the IfPack2::Preconditioner class.)  The Preconditioner class and derived classes sit in the `KokkosSparse::Experimental` namespace. 
+The `KokkosSparse_Preconditioner` class provides an abstract base class interface to use Kokkos-based preconditioners with iterative linear solvers. In particular, this class is designed to work with the Kokkos-based GMRES implementation in `examples/gmres`. It may also be useful for integrating Kokkos-based preconditioners into other solver codes and packages. (For Trilinos users: This class is loosely based upon the IfPack2::Preconditioner class.)  The Preconditioner class and derived classes sit in the `KokkosSparse::Experimental` namespace.
 
 ### Input parameters:
 
 ##### Preconditioner template paramters:
-The KokkosSparse Preconditioner has the following template parameters:   
-**ScalarType:** Type of scalars in the preconditioner apply. (double, float, half, Kokkos::complex<double>, etc.)   
-**Layout:** Layout of vectors X to which the preconditioner will be applied. (Kokkos::LayoutLeft, or other Kokkos layouts.)   
-**EXSP:** A Kokkos execution space.   
-**OrdinalType:** The ordinal type from the Kokkos::CrsMatrix A.   
-**Note:** If using this preconditioner with the Kokkos-based GMRES example, these template parameters must match those used to initialize the GMRES solver.  
+The KokkosSparse Preconditioner has the following template parameters:
+**CRS:** The type of the compressed row sparse matrix. All key types should be derivable from this.
 
 ### Preconditioner Base Class Functions (See code for full details.):
 
-**constructor**: Empty in the base class.   
+**constructor**: Empty in the base class.
 **apply**`( Kokkos::View<ScalarType*> &X, Kokkos::View<ScalarType*> &Y, const char transM[] = "N", ScalarType alpha = 1.0, ScalarType beta =  0):`
-Returns `y = beta * y + alpha * M * x` where `M` is the preconditioning operator.  (May apply `M` transposed if `transM = 'T' or 'C'` for 'transpose' and 'conjugate transpose' respectively.)   
-**setParameters():** Used to set preconditioner parameters.     
-**initialize():** Set up the graph structure of this preconditioner. If the graph structure of the constructor's input matrix has changed, or if you have not yet called `initialize()`, you must call `initialize()` before you may call `compute()` or `apply()`. Thus, `initialize()` corresponds to the "symbolic factorization" step of a sparse factorization, whether or not the specific preconditioner actually does a sparse factorization.   
-**compute():** Set up the numerical values in this preconditioner. If the values of the constructor's input matrix have changed, or if you have not yet called `compute()`, you must call `compute()` before you may call `apply()`. Thus, `compute()` corresponds to the "numeric factorization" step of a sparse factorization, whether or not the specific preconditioner actually does a sparse factorization.    
-**isInitialized()** and **isComputed()** return whether the preconditioner has been initialized or computed, respectively.  
+Returns `y = beta * y + alpha * M * x` where `M` is the preconditioning operator.  (May apply `M` transposed if `transM = 'T' or 'C'` for 'transpose' and 'conjugate transpose' respectively.)
+**setParameters():** Used to set preconditioner parameters.
+**initialize():** Set up the graph structure of this preconditioner. If the graph structure of the constructor's input matrix has changed, or if you have not yet called `initialize()`, you must call `initialize()` before you may call `compute()` or `apply()`. Thus, `initialize()` corresponds to the "symbolic factorization" step of a sparse factorization, whether or not the specific preconditioner actually does a sparse factorization.
+**compute():** Set up the numerical values in this preconditioner. If the values of the constructor's input matrix have changed, or if you have not yet called `compute()`, you must call `compute()` before you may call `apply()`. Thus, `compute()` corresponds to the "numeric factorization" step of a sparse factorization, whether or not the specific preconditioner actually does a sparse factorization.
+**isInitialized()** and **isComputed()** return whether the preconditioner has been initialized or computed, respectively.
 **hasTransposeApply()** Returns true if the transposed (or conjugate transposed) operator apply is available for this preconditioner.  Base class function returns `false`.
 
 ### Derived Preconditioner Classes:
 
 #### Matrix Prec:
-A simple class that takes a `KokkosSparse::CrsMatrix` to apply as the preconditioner `M`.  This matrix is given to the preconditioner constructor function.  There are no parameters to set.  The functions `initialize` and `compute` do not perform any operations.   
+A simple class that takes a `KokkosSparse::CrsMatrix` to apply as the preconditioner `M`.  This matrix is given to the preconditioner constructor function.  There are no parameters to set.  The functions `initialize` and `compute` do not perform any operations.
 
 ### Testing
 
 ##### Command-Line parameters for test\_prec:
 Current solver parameters for the GMRES preconditioning test are as follows:
 
-"**--mat-size**   :   The size of the n x n test matrix. (Default: n=1000.)   
-"**--max-subsp**   :  The maximum size of the Kyrlov subspace before restarting (Default 50)."   
-"**--max-restarts:**  Maximum number of GMRES restarts (Default 50)."   
-"**--tol        :**  Convergence tolerance.  (Default 1e-10)."   
-"**--ortho       :**  Type of orthogonalization. Use 'CGS2' or 'MGS'. (Default 'CGS2')"   
-"**--rand\_rhs**    :  Generate a random right-hand side b.  (Without this option, the solver default generates b = vector of ones.)"   
+"**--mat-size**   :   The size of the n x n test matrix. (Default: n=1000.)
+"**--max-subsp**   :  The maximum size of the Kyrlov subspace before restarting (Default 50)."
+"**--max-restarts:**  Maximum number of GMRES restarts (Default 50)."
+"**--tol        :**  Convergence tolerance.  (Default 1e-10)."
+"**--ortho       :**  Type of orthogonalization. Use 'CGS2' or 'MGS'. (Default 'CGS2')"
+"**--rand\_rhs**    :  Generate a random right-hand side b.  (Without this option, the solver default generates b = vector of ones.)"
 "**--help**         : Display a help message."
 
 #### Test Measurements:
@@ -48,4 +44,3 @@ Current solver parameters for the GMRES preconditioning test are as follows:
 If you wish to suggest an enhancement or make a bug report for this preconditioner code, please post an issue at https://github.com/kokkos/kokkos-kernels/issues or email jloe@sandia.gov.
 
 SAND2021-9744 O
-
diff --git a/example/gmres/ex_real_A.cpp b/example/gmres/ex_real_A.cpp
@@ -49,15 +49,19 @@
 #include <KokkosBlas.hpp>
 #include <KokkosBlas3_trsm.hpp>
 #include <KokkosSparse_spmv.hpp>
-
-#include "gmres.hpp"
+#include <KokkosSparse_gmres.hpp>
 
 int main(int argc, char* argv[]) {
-  typedef double ST;
-  typedef int OT;
-  typedef Kokkos::DefaultExecutionSpace EXSP;
+  using ST   = double;
+  using OT   = int;
+  using EXSP = Kokkos::DefaultExecutionSpace;
+  using MESP = typename EXSP::memory_space;
+  using CRS  = KokkosSparse::CrsMatrix<ST, OT, EXSP, void, OT>;
 
   using ViewVectorType = Kokkos::View<ST*, Kokkos::LayoutLeft, EXSP>;
+  using KernelHandle =
+      KokkosKernels::Experimental::KokkosKernelsHandle<OT, OT, ST, EXSP, MESP,
+                                                       MESP>;
 
   std::string filename("bcsstk09.mtx");  // example matrix
   std::string ortho("CGS2");             // orthog type
@@ -104,28 +108,29 @@ int main(int argc, char* argv[]) {
   std::cout << "File to process is: " << filename << std::endl;
   std::cout << "Convergence tolerance is: " << convTol << std::endl;
 
-  // Set GMRES options:
-  GmresOpts<ST> solverOpts;
-  solverOpts.tol        = convTol;
-  solverOpts.m          = m;
-  solverOpts.maxRestart = cycLim;
-  solverOpts.ortho      = ortho;
-  solverOpts.verbose    = false;  // No verbosity needed for most testing
-
   // Initialize Kokkos AFTER parsing parameters:
   Kokkos::initialize();
   {
     // Read in a matrix Market file and use it to test the Kokkos Operator.
-    KokkosSparse::CrsMatrix<ST, OT, EXSP> A =
-        KokkosSparse::Impl::read_kokkos_crst_matrix<
-            KokkosSparse::CrsMatrix<ST, OT, EXSP>>(filename.c_str());
+    CRS A = KokkosSparse::Impl::read_kokkos_crst_matrix<CRS>(filename.c_str());
 
     int n = A.numRows();
     ViewVectorType X("X", n);    // Solution and initial guess
     ViewVectorType Wj("Wj", n);  // For checking residuals at end.
     ViewVectorType B(Kokkos::view_alloc(Kokkos::WithoutInitializing, "B"),
                      n);  // right-hand side vec
 
+    // Make kernel handles and set options
+    KernelHandle kh;
+    kh.create_gmres_handle(m, convTol, cycLim);
+    auto gmres_handle = kh.get_gmres_handle();
+    // Get full gmres handle type using decltype. Deferencing a pointer gives a
+    // reference, so we need to strip that too.
+    using GMRESHandle =
+        typename std::remove_reference<decltype(*gmres_handle)>::type;
+    gmres_handle->set_ortho(ortho == "CGS2" ? GMRESHandle::Ortho::CGS2
+                                            : GMRESHandle::Ortho::MGS);
+
     if (rand_rhs) {
       // Make rhs random.
       int rand_seed = 123;
@@ -137,8 +142,11 @@ int main(int argc, char* argv[]) {
     }
 
     // Run GMRS solve:
-    GmresStats solveStats =
-        gmres<ST, Kokkos::LayoutLeft, EXSP>(A, B, X, solverOpts);
+    KokkosSparse::Experimental::gmres(&kh, A, B, X);
+
+    const auto numIters  = gmres_handle->get_num_iters();
+    const auto convFlag  = gmres_handle->get_conv_flag_val();
+    const auto endRelRes = gmres_handle->get_end_rel_res();
 
     // Double check residuals at end of solve:
     ST nrmB = KokkosBlas::nrm2(B);
@@ -147,11 +155,10 @@ int main(int argc, char* argv[]) {
     ST endRes = KokkosBlas::nrm2(B) / nrmB;
     std::cout << "=========================================" << std::endl;
     std::cout << "Verify from main: Ending residual is " << endRes << std::endl;
-    std::cout << "Number of iterations is: " << solveStats.numIters
-              << std::endl;
+    std::cout << "Number of iterations is: " << numIters << std::endl;
     std::cout << "Diff of residual from main - residual from solver: "
-              << solveStats.endRelRes - endRes << std::endl;
-    std::cout << "Convergence flag is : " << solveStats.convFlag() << std::endl;
+              << endRelRes - endRes << std::endl;
+    std::cout << "Convergence flag is : " << convFlag << std::endl;
   }
   Kokkos::finalize();
 }