add example

examples/CMakeLists.txt

@@ -30,6 +30,7 @@ set(EXAMPLES_LIST
     par-ilu-convergence
     performance-debugging
     preconditioner-export
+    reordered-preconditioned-solver
     simple-solver-logging)
 
 if(GINKGO_BUILD_CUDA AND GINKGO_BUILD_OMP)

examples/reordered-preconditioned-solver/CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 3.16)
+project(reordered-preconditioned-solver)
+
+# We only need to find Ginkgo if we build this example stand-alone
+if (NOT GINKGO_BUILD_EXAMPLES)
+    find_package(Ginkgo 1.8.0 REQUIRED)
+endif()
+add_executable(reordered-preconditioned-solver reordered-preconditioned-solver.cpp)
+target_link_libraries(reordered-preconditioned-solver Ginkgo::ginkgo)
+
+# Copy the data files to the execution directory
+configure_file(data/A.mtx data/A.mtx COPYONLY)
+configure_file(data/b.mtx data/b.mtx COPYONLY)
+configure_file(data/x0.mtx data/x0.mtx COPYONLY)

examples/reordered-preconditioned-solver/build.sh

@@ -0,0 +1,17 @@
+#!/bin/bash
+
+# set up script
+if [ $# -ne 1 ]; then
+    echo -e "Usage: $0 GINKGO_BUILD_DIRECTORY"
+    exit 1
+fi
+BUILD_DIR=$1
+THIS_DIR=$( cd "$( dirname "${BASH_SOURCE[0]}" )" &>/dev/null && pwd )
+
+source ${THIS_DIR}/../build-setup.sh
+
+# build
+${CXX} -std=c++14 -o ${THIS_DIR}/preconditioned-solver \
+    ${THIS_DIR}/reordered-preconditioned-solver.cpp \
+    -I${THIS_DIR}/../../include -I${BUILD_DIR}/include \
+    -L${THIS_DIR} ${LINK_FLAGS}

examples/reordered-preconditioned-solver/data/A.mtx

@@ -0,0 +1,114 @@
+%%MatrixMarket matrix coordinate integer symmetric
+%-------------------------------------------------------------------------------
+% UF Sparse Matrix Collection, Tim Davis
+% http://www.cise.ufl.edu/research/sparse/matrices/JGD_Trefethen/Trefethen_20b
+% name: JGD_Trefethen/Trefethen_20b
+% [Diagonal matrices with primes, Nick Trefethen, Oxford Univ.]
+% id: 2203
+% date: 2008
+% author: N. Trefethen
+% ed: J.-G. Dumas
+% fields: name title A id date author ed kind notes
+% kind: combinatorial problem
+%-------------------------------------------------------------------------------
+% notes:
+% Diagonal matrices with primes, Nick Trefethen, Oxford Univ.          
+% From Jean-Guillaume Dumas' Sparse Integer Matrix Collection,         
+% http://ljk.imag.fr/membres/Jean-Guillaume.Dumas/simc.html            
+%                                                                      
+% Problem 7 of the Hundred-dollar, Hundred-digit Challenge Problems,   
+% SIAM News, vol 35, no. 1.                                            
+%                                                                      
+% 7. Let A be the 20,000 x 20,000 matrix whose entries are zero        
+% everywhere except for the primes 2, 3, 5, 7, . . . , 224737 along the
+% main diagonal and the number 1 in all the positions A(i,j) with      
+% |i-j| = 1,2,4,8, . . . ,16384.  What is the (1,1) entry of inv(A)?   
+%                                                                      
+% http://www.siam.org/news/news.php?id=388                             
+%                                                                      
+% Filename in JGD collection: Trefethen/trefethen_20__19_minor.sms     
+%-------------------------------------------------------------------------------
+19 19 83
+1 1 3
+2 1 1
+3 1 1
+5 1 1
+9 1 1
+17 1 1
+2 2 5
+3 2 1
+4 2 1
+6 2 1
+10 2 1
+18 2 1
+3 3 7
+4 3 1
+5 3 1
+7 3 1
+11 3 1
+19 3 1
+4 4 11
+5 4 1
+6 4 1
+8 4 1
+12 4 1
+5 5 13
+6 5 1
+7 5 1
+9 5 1
+13 5 1
+6 6 17
+7 6 1
+8 6 1
+10 6 1
+14 6 1
+7 7 19
+8 7 1
+9 7 1
+11 7 1
+15 7 1
+8 8 23
+9 8 1
+10 8 1
+12 8 1
+16 8 1
+9 9 29
+10 9 1
+11 9 1
+13 9 1
+17 9 1
+10 10 31
+11 10 1
+12 10 1
+14 10 1
+18 10 1
+11 11 37
+12 11 1
+13 11 1
+15 11 1
+19 11 1
+12 12 41
+13 12 1
+14 12 1
+16 12 1
+13 13 43
+14 13 1
+15 13 1
+17 13 1
+14 14 47
+15 14 1
+16 14 1
+18 14 1
+15 15 53
+16 15 1
+17 15 1
+19 15 1
+16 16 59
+17 16 1
+18 16 1
+17 17 61
+18 17 1
+19 17 1
+18 18 67
+19 18 1
+19 19 71

examples/reordered-preconditioned-solver/data/b.mtx

@@ -0,0 +1,21 @@
+%%MatrixMarket matrix array real general
+19 1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1

examples/reordered-preconditioned-solver/data/x0.mtx

@@ -0,0 +1,21 @@
+%%MatrixMarket matrix array real general
+19 1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0

examples/reordered-preconditioned-solver/doc/builds-on

@@ -0,0 +1 @@
+preconditioned-solver

examples/reordered-preconditioned-solver/doc/intro.dox

@@ -0,0 +1,5 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+<h3> About the example </h3>
+

examples/reordered-preconditioned-solver/doc/kind

@@ -0,0 +1 @@
+preconditioners

examples/reordered-preconditioned-solver/doc/results.dox

@@ -0,0 +1,35 @@
+<h1>Results</h1>
+This is the expected output:
+
+@code{.cpp}
+
+Solution (x):
+%%MatrixMarket matrix array real general
+19 1
+0.252218
+0.108645
+0.0662811
+0.0630433
+0.0384088
+0.0396536
+0.0402648
+0.0338935
+0.0193098
+0.0234653
+0.0211499
+0.0196413
+0.0199151
+0.0181674
+0.0162722
+0.0150714
+0.0107016
+0.0121141
+0.0123025
+Residual norm sqrt(r^T r):
+%%MatrixMarket matrix array real general
+1 1
+4.82005e-08
+
+@endcode
+
+<h3> Comments about programming and debugging </h3>

examples/reordered-preconditioned-solver/doc/short-intro

@@ -0,0 +1 @@
+The reordered preconditioned solver example.

examples/reordered-preconditioned-solver/doc/tooltip

@@ -0,0 +1 @@
+Use a reordering in Ginkgo. Solve a linear system with a preconditioner.

examples/reordered-preconditioned-solver/reordered-preconditioned-solver.cpp

@@ -0,0 +1,111 @@
+// SPDX-FileCopyrightText: 2017-2023 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+
+#include <ginkgo/ginkgo.hpp>
+
+
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <string>
+#include "ginkgo/core/matrix/permutation.hpp"
+
+
+int main(int argc, char* argv[])
+{
+    // Some shortcuts
+    using ValueType = double;
+    using RealValueType = gko::remove_complex<ValueType>;
+    using IndexType = int;
+
+    using vec = gko::matrix::Dense<ValueType>;
+    using real_vec = gko::matrix::Dense<RealValueType>;
+    using mtx = gko::matrix::Csr<ValueType, IndexType>;
+    using cg = gko::solver::Cg<ValueType>;
+    using bj = gko::preconditioner::Jacobi<ValueType, IndexType>;
+
+    // Print version information
+    std::cout << gko::version_info::get() << std::endl;
+
+    // Figure out where to run the code
+    if (argc == 2 && (std::string(argv[1]) == "--help")) {
+        std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
+        std::exit(-1);
+    }
+
+    // Figure out where to run the code
+    const auto executor_string = argc >= 2 ? argv[1] : "reference";
+    std::map<std::string, std::function<std::shared_ptr<gko::Executor>()>>
+        exec_map{
+            {"omp", [] { return gko::OmpExecutor::create(); }},
+            {"cuda",
+             [] {
+                 return gko::CudaExecutor::create(0,
+                                                  gko::OmpExecutor::create());
+             }},
+            {"hip",
+             [] {
+                 return gko::HipExecutor::create(0, gko::OmpExecutor::create());
+             }},
+            {"dpcpp",
+             [] {
+                 return gko::DpcppExecutor::create(0,
+                                                   gko::OmpExecutor::create());
+             }},
+            {"reference", [] { return gko::ReferenceExecutor::create(); }}};
+
+    // executor where Ginkgo will perform the computation
+    const auto exec = exec_map.at(executor_string)();  // throws if not valid
+
+    // Read data
+    auto A = share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));
+    auto b = gko::read<vec>(std::ifstream("data/b.mtx"), exec);
+    auto x = gko::read<vec>(std::ifstream("data/x0.mtx"), exec);
+
+    auto reordering =
+        gko::experimental::reorder::Rcm<IndexType>::build().on(exec)->generate(
+            A);
+    // Permute matrix and vectors
+    auto A_reordered = share(A->permute(reordering));
+    auto b_reordered = b->permute(reordering, gko::matrix::permute_mode::rows);
+    // this reordering is not necessary, but it maps the initial guess to the
+    // unreordered case
+    auto x_reordered = x->permute(reordering, gko::matrix::permute_mode::rows);
+
+    const RealValueType reduction_factor{1e-7};
+    // Create solver factory
+    auto solver_gen =
+        cg::build()
+            .with_criteria(gko::stop::Iteration::build().with_max_iters(20u),
+                           gko::stop::ResidualNorm<ValueType>::build()
+                               .with_reduction_factor(reduction_factor))
+            // Add preconditioner, these 2 lines are the only
+            // difference from the simple solver example
+            .with_preconditioner(bj::build().with_max_block_size(8u))
+            .on(exec);
+    // Create solver
+    auto solver = solver_gen->generate(A_reordered);
+
+    // Solve system
+    solver->apply(b_reordered, x_reordered);
+
+    // Revert permutation to get the unpermuted solution
+    x_reordered->permute(reordering, x,
+                         gko::matrix::permute_mode::inverse_rows);
+
+    // Print solution
+    std::cout << "Solution (x):\n";
+    write(std::cout, x);
+
+    // Calculate residual
+    auto one = gko::initialize<vec>({1.0}, exec);
+    auto neg_one = gko::initialize<vec>({-1.0}, exec);
+    auto res = gko::initialize<real_vec>({0.0}, exec);
+    A->apply(one, x, neg_one, b);
+    b->compute_norm2(res);
+
+    std::cout << "Residual norm sqrt(r^T r):\n";
+    write(std::cout, res);
+}