Add another RHS generation in the benchmark

- The benchmark script now has the option to change the initial guess - Add option to generate the RHS with: b = A * (s / |s|) with s(i) = sin(i)
ginkgo-project · thoasm · Feb 21, 2021 · Apr 27, 2020 · Apr 27, 2020 · Apr 29, 2020
commit 155c552411557df5f453cb7427bc93aabedfa2ae
diff --git a/BENCHMARKING.md b/BENCHMARKING.md
@@ -143,7 +143,7 @@ The benchmark suite can take a number of configuration parameters. Benchmarks
 can be run only for `sparse matrix vector products (spmv)`, for full solvers
 (with or without preconditioners), or for preconditioners only when supported.
 The benchmark suite also allows to target a sub-part of the SuiteSparse matrix
-collection. For details, see the [available benchmark options](### 5: Available
+collection. For details, see the [available benchmark options](### 6: Available
 benchmark options). Here are the most important options:
 * `BENCHMARK={spmv, solver, preconditioner}` - allows to select the type of
     benchmark to be ran.
@@ -303,9 +303,16 @@ The supported environment variables are described in the following list:
     the solver should stop. The default is `1e-6`.
 * `SOLVERS_MAX_ITERATION=<number>` - the maximum number of iterations with which
     a solver should be ran. The default is `10000`.
-* `SOLVERS_RHS={unit, random}` - whether to use a vector of all ones or random
-    values as the right-hand side in solver benchmarks. Default is `unit`.
+* `SOLVERS_RHS={1, random, sinus}` - whether to use a vector of all ones,
+    random values or $b = A * (s / |s|)$ with $s(i) = sin(i)$
+    as the right-hand side in solver benchmarks. Default is `1`.
+* `SOLVERS_INITIAL_GUESS`={rhs,0,random} - the initial guess generation of the
+    solvers. `rhs` uses the right-hand side, `0` uses a zero vector and `random`
+    generates a random vector as the initial guess.
 * `DETAILED={0,1}` - selects whether detailed benchmarks should be ran for the
     solver benchmarks, can be either `0` (off) or `1` (on). The default is `0`.
 * `GPU_TIMER={true, false}` - If set to `true`, use the gpu timer, which is
     valid for cuda/hip executor, to measure the timing. Default is `false`.
+* `SOLVERS_JACOBI_MAX_BS` - sets the maximum block size for the Jacobi
+    preconditioner (if used, otherwise, it does nothing) in the solvers
+    benchmark
diff --git a/benchmark/run_all_benchmarks.sh b/benchmark/run_all_benchmarks.sh
@@ -60,9 +60,15 @@ if [ ! "${DEVICE_ID}" ]; then
     DEVICE_ID="0"
 fi
 
+if [ ! "${SOLVERS_JACOBI_MAX_BS}" ]; then
+    SOLVERS_JACOBI_MAX_BS="32"
+    "SOLVERS_JACOBI_MAX_BS environment variable not set - assuming \"${SOLVERS_JACOBI_MAX_BS}\"" 1>&2
+fi
+
+
 if [ ! "${SOLVERS_RHS}" ]; then
-    echo "SOLVERS_RHS environment variable not set - assuming \"unit\"" 1>&2
-    SOLVERS_RHS="unit"
+    SOLVERS_RHS="1"
+    echo "SOLVERS_RHS environment variable not set - assuming \"${SOLVERS_RHS}\"" 1>&2
 fi
 
 if [ ! "${BENCHMARK_PRECISION}" ]; then
@@ -85,16 +91,39 @@ else
 fi
 
 if [ "${SOLVERS_RHS}" == "random" ]; then
-    SOLVERS_RHS_FLAG="--random_rhs=true"
+    SOLVERS_RHS_FLAG="--rhs_generation=random"
+elif [ "${SOLVERS_RHS}" == "1" ]; then
+    SOLVERS_RHS_FLAG="--rhs_generation=1"
+elif [ "${SOLVERS_RHS}" == "sinus" ]; then
+    SOLVERS_RHS_FLAG="--rhs_generation=sinus"
 else
-    SOLVERS_RHS_FLAG="--random_rhs=false"
+    echo "SOLVERS_RHS does not support the value \"${SOLVERS_RHS}\"." 1>&2
+    echo "The following values are supported: \"1\", \"random\" and \"sinus\"" 1>&2
+    exit 1
+fi
+
+if [ ! "${SOLVERS_INITIAL_GUESS}" ]; then
+    SOLVERS_INITIAL_GUESS="rhs"
+    echo "SOLVERS_RHS environment variable not set - assuming \"${SOLVERS_INITIAL_GUESS}\"" 1>&2
 fi
 
 if [ ! "${GPU_TIMER}" ]; then
     echo "GPU_TIMER    environment variable not set - assuming \"false\"" 1>&2
     GPU_TIMER="false"
 fi
 
+if [ "${SOLVERS_INITIAL_GUESS}" == "random" ]; then
+    SOLVERS_INITIAL_GUESS_FLAG="--initial_guess_generation=random"
+elif [ "${SOLVERS_INITIAL_GUESS}" == "0" ]; then
+    SOLVERS_INITIAL_GUESS_FLAG="--initial_guess_generation=0"
+elif [ "${SOLVERS_INITIAL_GUESS}" == "rhs" ]; then
+    SOLVERS_INITIAL_GUESS_FLAG="--initial_guess_generation=rhs"
+else
+    echo "SOLVERS_RHS does not support the value \"${SOLVERS_RHS}\"." 1>&2
+    echo "The following values are supported: \"1\", \"random\" and \"sinus\"" 1>&2
+    exit 1
+fi
+
 # Control whether to run detailed benchmarks or not.
 # Default setting is detailed=false. To activate, set DETAILED=1.
 if  [ ! "${DETAILED}" ] || [ "${DETAILED}" -eq 0 ]; then
@@ -202,7 +231,9 @@ run_solver_benchmarks() {
                     --executor="${EXECUTOR}" --solvers="${SOLVERS}" \
                     --preconditioners="${PRECONDS}" \
                     --max_iters=${SOLVERS_MAX_ITERATIONS} --rel_res_goal=${SOLVERS_PRECISION} \
-                    ${SOLVERS_RHS_FLAG} ${DETAILED_STR} --device_id="${DEVICE_ID}" --gpu_timer=${GPU_TIMER} \
+                    ${SOLVERS_RHS_FLAG} ${DETAILED_STR} ${SOLVERS_INITIAL_GUESS_FLAG} \
+                    --gpu_timer=${GPU_TIMER} \
+                    --jacobi_max_block_size=${SOLVERS_JACOBI_MAX_BS} --device_id="${DEVICE_ID}" \
                     <"$1.imd" 2>&1 >"$1"
     keep_latest "$1" "$1.bkp" "$1.bkp2" "$1.imd"
 }

diff --git a/benchmark/solver/solver.cpp b/benchmark/solver/solver.cpp
@@ -89,11 +89,19 @@ DEFINE_double(
     idr_kappa, 0.7,
     "the number to check whether Av_n and v_n are too close or not in IDR");
 
-DEFINE_bool(random_rhs, false,
-            "Use a random vector for the rhs (otherwise use all ones)");
-
-DEFINE_bool(random_initial_guess, false,
-            "Use a random vector for the initial guess (otherwise use rhs)");
+DEFINE_string(
+    rhs_generation, "1",
+    "Method used to generate the right hand side. Supported values are:"
+    "1, random, sinus. 1 sets all values of the right hand side to 1, "
+    "random assigns the values to a uniformly distributed random number "
+    "in [-1, 1), and sinus assigns b = A * (s / |s|) with A := system matrix, "
+    "s := vector with s(i) = sin(i).");
+
+DEFINE_string(
+    initial_guess_generation, "rhs",
+    "Method used to generate the initial guess. Supported values are: "
+    "random, rhs, 0. random uses a random vector, rhs uses the right "
+    "hand side, and 0 uses a zero vector as the initial guess.");
 
 // This allows to benchmark the overhead of a solver by using the following
 // data: A=[1.0], x=[0.0], b=[nan]. This data can be used to benchmark normal
@@ -118,6 +126,62 @@ DEFINE_bool(overhead, false,
 }
 
 
+template <typename Engine>
+std::unique_ptr<vec<etype>> generate_rhs(
+    std::shared_ptr<const gko::Executor> exec,
+    std::shared_ptr<const gko::LinOp> system_matrix, Engine engine)
+{
+    using rc_etype = gko::remove_complex<etype>;
+    gko::dim<2> vec_size{system_matrix->get_size()[0], FLAGS_nrhs};
+    if (FLAGS_rhs_generation == "1") {
+        return create_matrix<etype>(exec, vec_size, gko::one<etype>());
+    } else if (FLAGS_rhs_generation == "random") {
+        return create_matrix<etype>(exec, vec_size, engine);
+    } else if (FLAGS_rhs_generation == "sinus") {
+        auto rhs = vec<etype>::create(exec, vec_size);
+
+        auto tmp = create_matrix_sin<etype>(exec, vec_size);
+        auto scalar = gko::matrix::Dense<rc_etype>::create(
+            exec->get_master(), gko::dim<2>{1, vec_size[1]});
+        tmp->compute_norm2(scalar.get());
+        for (gko::size_type i = 0; i < vec_size[1]; ++i) {
+            scalar->at(0, i) = gko::one<rc_etype>() / scalar->at(0, i);
+        }
+        // normalize sin-vector
+        if (gko::is_complex_s<etype>::value) {
+            tmp->scale(scalar->make_complex().get());
+        } else {
+            tmp->scale(scalar.get());
+        }
+        system_matrix->apply(tmp.get(), rhs.get());
+        return rhs;
+    }
+    throw std::invalid_argument(std::string("\"rhs_generation\" = ") +
+                                FLAGS_rhs_generation + " is not supported!");
+}
+
+
+template <typename Engine>
+std::unique_ptr<vec<etype>> generate_initial_guess(
+    std::shared_ptr<const gko::Executor> exec,
+    std::shared_ptr<const gko::LinOp> system_matrix, const vec<etype> *rhs,
+    Engine engine)
+{
+    using rc_etype = gko::remove_complex<etype>;
+    gko::dim<2> vec_size{system_matrix->get_size()[1], FLAGS_nrhs};
+    if (FLAGS_initial_guess_generation == "0") {
+        return create_matrix<etype>(exec, vec_size, gko::zero<etype>());
+    } else if (FLAGS_initial_guess_generation == "random") {
+        return create_matrix<etype>(exec, vec_size, engine);
+    } else if (FLAGS_initial_guess_generation == "rhs") {
+        return rhs->clone();
+    }
+    throw std::invalid_argument(std::string("\"initial_guess_generation\" = ") +
+                                FLAGS_initial_guess_generation +
+                                " is not supported!");
+}
+
+
 void validate_option_object(const rapidjson::Value &value)
 {
     if (!value.IsObject() || !value.HasMember("optimal") ||
@@ -552,19 +616,10 @@ int main(int argc, char *argv[])
                     std::ifstream rhs_fd{test_case["rhs"].GetString()};
                     b = gko::read<Vec>(rhs_fd, exec);
                 } else {
-                    b = create_matrix<etype>(
-                        exec,
-                        gko::dim<2>{system_matrix->get_size()[0], FLAGS_nrhs},
-                        engine, FLAGS_random_rhs);
-                }
-                if (FLAGS_random_initial_guess) {
-                    x = create_matrix<etype>(
-                        exec,
-                        gko::dim<2>{system_matrix->get_size()[0], FLAGS_nrhs},
-                        engine);
-                } else {
-                    x = b->clone();
+                    b = generate_rhs(exec, system_matrix, engine);
                 }
+                x = generate_initial_guess(exec, system_matrix, b.get(),
+                                           engine);
             }
 
             std::clog << "Matrix is of size (" << system_matrix->get_size()[0]

diff --git a/benchmark/utils/general.hpp b/benchmark/utils/general.hpp
@@ -287,22 +287,51 @@ template <typename ValueType>
 using vec = gko::matrix::Dense<ValueType>;
 
 
-// creates a zero vector
+// Create a matrix with value indices s[i, j] = sin(i)
 template <typename ValueType>
-std::unique_ptr<vec<ValueType>> create_vector(
-    std::shared_ptr<const gko::Executor> exec, gko::size_type size)
+std::enable_if_t<!gko::is_complex_s<ValueType>::value,
+                 std::unique_ptr<vec<ValueType>>>
+create_matrix_sin(std::shared_ptr<const gko::Executor> exec, gko::dim<2> size)
 {
+    auto h_res = vec<ValueType>::create(exec->get_master(), size);
+    for (gko::size_type i = 0; i < size[0]; ++i) {
+        for (gko::size_type j = 0; j < size[1]; ++j) {
+            h_res->at(i, j) = std::sin(static_cast<ValueType>(i));
+        }
+    }
     auto res = vec<ValueType>::create(exec);
-    res->read(gko::matrix_data<ValueType>(gko::dim<2>{size, 1}));
+    h_res->move_to(res.get());
     return res;
 }
 
+// Note: complex values are assigned s[i, j] = {sin(2 * i, sin(2 * i + 1))}
+template <typename ValueType>
+std::enable_if_t<gko::is_complex_s<ValueType>::value,
+                 std::unique_ptr<vec<ValueType>>>
+create_matrix_sin(std::shared_ptr<const gko::Executor> exec, gko::dim<2> size)
+{
+    auto h_res = vec<ValueType>::create(exec->get_master(), size);
+    for (gko::size_type i = 0; i < size[0]; ++i) {
+        for (gko::size_type j = 0; j < size[1]; ++j) {
+            h_res->at(i, j) = ValueType{
+                std::sin(static_cast<gko::remove_complex<ValueType>>(2 * i)),
+                std::sin(
+                    static_cast<gko::remove_complex<ValueType>>(2 * i + 1))};
+        }
+    }
+    auto res = vec<ValueType>::create(exec);
+    h_res->move_to(res.get());
+    return res;
+}
+
+
 template <typename ValueType>
 std::unique_ptr<vec<ValueType>> create_matrix(
-    std::shared_ptr<const gko::Executor> exec, gko::dim<2> size)
+    std::shared_ptr<const gko::Executor> exec, gko::dim<2> size,
+    ValueType value)
 {
     auto res = vec<ValueType>::create(exec);
-    res->read(gko::matrix_data<ValueType>(size));
+    res->read(gko::matrix_data<ValueType>(size, value));
     return res;
 }
 
@@ -311,18 +340,25 @@ std::unique_ptr<vec<ValueType>> create_matrix(
 template <typename ValueType, typename RandomEngine>
 std::unique_ptr<vec<ValueType>> create_matrix(
     std::shared_ptr<const gko::Executor> exec, gko::dim<2> size,
-    RandomEngine &engine, bool random = true)
+    RandomEngine &engine)
 {
     auto res = vec<ValueType>::create(exec);
-    if (random) {
-        res->read(gko::matrix_data<ValueType>(
-            size,
-            std::uniform_real_distribution<gko::remove_complex<ValueType>>(-1.0,
-                                                                           1.0),
-            engine));
-    } else {
-        res->read(gko::matrix_data<ValueType>(size, gko::one<ValueType>()));
-    }
+    res->read(gko::matrix_data<ValueType>(
+        size,
+        std::uniform_real_distribution<gko::remove_complex<ValueType>>(-1.0,
+                                                                       1.0),
+        engine));
+    return res;
+}
+
+
+// creates a zero vector
+template <typename ValueType>
+std::unique_ptr<vec<ValueType>> create_vector(
+    std::shared_ptr<const gko::Executor> exec, gko::size_type size)
+{
+    auto res = vec<ValueType>::create(exec);
+    res->read(gko::matrix_data<ValueType>(gko::dim<2>{size, 1}));
     return res;
 }