Update tutorials (mostly CoupCons3D one)

docs/Serena.rst

@@ -20,7 +20,7 @@ portrait is shown on the figure below.
 .. _examples/mm2bin: https://github.com/ddemidov/amgcl/blob/master/examples/mm2bin.cpp
 
 .. figure:: ../tutorial/2.Serena/Serena.png
-   :width: 80%
+   :width: 90%
    :name: serena_spy
 
    Serena matrix portrait
@@ -186,9 +186,11 @@ applications per iteration, while CG only does one.
 The problem description states that *the medium is strongly heterogeneous and
 characterized by a complex geometry consisting of alternating sequences of thin
 clay and sand layers*. This may result in high contrast between matrix
-coefficients in the neighboring rows, so lets see if scaling the matrix (so
-that it has the unit diagonal) helps with the convergence. The ``-s`` option
-tells the solver to do that::
+coefficients in the neighboring rows, which is confirmed by the plot of the
+matrix diagonal in :ref:`serena_spy`: the diagonal coefficients span more than
+10 orders of magnitude! Scaling the matrix (so that it has the unit diagonal)
+should help with the convergence. The ``-s`` option tells the solver to do
+that::
 
     $ solver -B -A Serena.bin -s \
           solver.type=cg solver.maxiter=200 \

docs/poisson3Db.rst

@@ -10,7 +10,7 @@ matrix has an interesting structure, presented on the figure below:
 .. _Matrix Market: https://math.nist.gov/MatrixMarket
 
 .. figure:: ../tutorial/1.poisson3Db/Poisson3D.png
-   :width: 80%
+   :width: 90%
 
    Poisson3Db matrix portrait
 

tutorial/1.poisson3Db/plot.py

@@ -3,7 +3,11 @@
 from scipy.io import mmread
 
 A = mmread('poisson3Db.mtx')
-figure(figsize=(8,8))
-spy(A, marker='.', markersize=0.25, alpha=0.2)
+
+fig, (ax1, ax2) = subplots(2, 1, sharex=True, figsize=(8,10), gridspec_kw=dict(height_ratios=[4,1]))
+ax1.spy(A, marker='.', markersize=0.25, alpha=0.2)
+ax2.semilogy(A.diagonal())
+ax2.set_ylabel('Diagonal')
 tight_layout()
+
 savefig('Poisson3D.png')

tutorial/2.Serena/plot.py

@@ -3,17 +3,21 @@
 from scipy.io import mmread
 
 A = mmread('Serena.mtx')
-figure(figsize=(8,8))
-spy(A, marker='.', markersize=0.25, alpha=0.2)
-ax = gca()
-axins = ax.inset_axes([0.55, 0.55, 0.4, 0.4])
+
+fig, (ax1, ax2) = subplots(2, 1, sharex=True, figsize=(8,10), gridspec_kw=dict(height_ratios=[4,1]))
+ax1.spy(A, marker='.', markersize=0.25, alpha=0.2)
+axins = ax1.inset_axes([0.55, 0.55, 0.4, 0.4])
 axins.spy(A, marker='o', markersize=3, alpha=0.5)
 n = (A.shape[0] // 3 // 3) * 3
 axins.set_xlim([n - 0.5, n + 44.5])
 axins.set_ylim([n - 0.5, n + 44.5])
 axins.invert_yaxis()
 axins.set_xticklabels('')
 axins.set_yticklabels('')
-ax.indicate_inset_zoom(axins)
+ax1.indicate_inset_zoom(axins)
+
+ax2.semilogy(A.diagonal())
+ax2.set_ylabel('Diagonal')
+
 tight_layout()
 savefig('Serena.png')

tutorial/3.CoupCons3D/CMakeLists.txt

@@ -1,7 +1,10 @@
 add_executable(coupcons3d coupcons3d.cpp)
 target_link_libraries(coupcons3d amgcl)
 
-#if (VexCL_FOUND)
-#    vexcl_add_executables(coupcons3d_vexcl coupcons3d_vexcl.cpp)
-#    target_link_libraries(coupcons3d_vexcl INTERFACE amgcl)
-#endif()
+add_executable(coupcons3d_spc coupcons3d_spc.cpp)
+target_link_libraries(coupcons3d_spc amgcl)
+
+if (VexCL_FOUND)
+    vexcl_add_executables(coupcons3d_vexcl coupcons3d_vexcl.cpp)
+    target_link_libraries(coupcons3d_vexcl INTERFACE amgcl)
+endif()

tutorial/3.CoupCons3D/coupcons3d.cpp

@@ -1,34 +1,28 @@
+#include <vector>
 #include <iostream>
-#include <string>
 
 #include <amgcl/backend/builtin.hpp>
 #include <amgcl/adapter/crs_tuple.hpp>
-#include <amgcl/value_type/static_matrix.hpp>
-#include <amgcl/adapter/block_matrix.hpp>
-#include <amgcl/preconditioner/schur_pressure_correction.hpp>
 #include <amgcl/make_solver.hpp>
-#include <amgcl/make_block_solver.hpp>
 #include <amgcl/amg.hpp>
-#include <amgcl/solver/bicgstab.hpp>
-#include <amgcl/solver/preonly.hpp>
-#include <amgcl/coarsening/aggregation.hpp>
+#include <amgcl/coarsening/smoothed_aggregation.hpp>
 #include <amgcl/relaxation/ilu0.hpp>
-#include <amgcl/relaxation/spai0.hpp>
-#include <amgcl/relaxation/as_preconditioner.hpp>
+#include <amgcl/solver/bicgstab.hpp>
+#include <amgcl/value_type/static_matrix.hpp>
+#include <amgcl/adapter/block_matrix.hpp>
 
 #include <amgcl/io/mm.hpp>
 #include <amgcl/profiler.hpp>
 
-//---------------------------------------------------------------------------
 int main(int argc, char *argv[]) {
     // The matrix and the RHS file names should be in the command line options:
-    if (argc < 3) {
-        std::cerr << "Usage: " << argv[0] << " <matrix.mtx> <nu>" << std::endl;
+    if (argc < 2) {
+        std::cerr << "Usage: " << argv[0] << " <matrix.mtx>" << std::endl;
         return 1;
     }
 
     // The profiler:
-    amgcl::profiler<> prof("CoupCons3D");
+    amgcl::profiler<> prof("Serena");
 
     // Read the system matrix:
     size_t rows, cols;
@@ -43,49 +37,54 @@ int main(int argc, char *argv[]) {
     // The RHS is filled with ones:
     std::vector<double> f(rows, 1.0);
 
-    // The number of unknowns in the U subsystem
-    size_t nu = std::stoi(argv[2]);
+    // Scale the matrix so that it has the unit diagonal.
+    // First, find the diagonal values:
+    std::vector<double> D(rows, 1.0);
+    for(size_t i = 0; i < rows; ++i) {
+        for(ptrdiff_t j = ptr[i], e = ptr[i+1]; j < e; ++j) {
+            if (col[j] == i) {
+                D[i] = 1 / sqrt(val[j]);
+                break;
+            }
+        }
+    }
+
+    // Then, apply the scaling in-place:
+    for(size_t i = 0; i < rows; ++i) {
+        for(ptrdiff_t j = ptr[i], e = ptr[i+1]; j < e; ++j) {
+            val[j] *= D[i] * D[col[j]];
+        }
+        f[i] *= D[i];
+    }
 
     // We use the tuple of CRS arrays to represent the system matrix.
     // Note that std::tie creates a tuple of references, so no data is actually
     // copied here:
     auto A = std::tie(rows, ptr, col, val);
 
     // Compose the solver type
-    typedef amgcl::backend::builtin<double> SBackend; // the outer iterative solver backend
-    typedef amgcl::backend::builtin<float> PBackend;  // the PSolver backend
-    typedef amgcl::backend::builtin<
-        amgcl::static_matrix<float,4,4>> UBackend;    // the USolver backend
+    typedef amgcl::static_matrix<double, 4, 4> dmat_type; // matrix value type in double precision
+    typedef amgcl::static_matrix<double, 4, 1> dvec_type; // the corresponding vector value type
+    typedef amgcl::static_matrix<float,  4, 4> smat_type; // matrix value type in single precision
+
+    typedef amgcl::backend::builtin<dmat_type> SBackend; // the solver backend
+    typedef amgcl::backend::builtin<smat_type> PBackend; // the preconditioner backend
 
     typedef amgcl::make_solver<
-        amgcl::preconditioner::schur_pressure_correction<
-            amgcl::make_block_solver<
-                amgcl::amg<
-                    UBackend,
-                    amgcl::coarsening::aggregation,
-                    amgcl::relaxation::ilu0
-                    >,
-                amgcl::solver::preonly<UBackend>
-                >,
-            amgcl::make_solver<
-                amgcl::relaxation::as_preconditioner<
-                    PBackend,
-                    amgcl::relaxation::spai0
-                    >,
-                amgcl::solver::preonly<PBackend>
-                >
+        amgcl::amg<
+            PBackend,
+            amgcl::coarsening::smoothed_aggregation,
+            amgcl::relaxation::ilu0
             >,
         amgcl::solver::bicgstab<SBackend>
         > Solver;
 
-    // Solver parameters
-    Solver::params prm;
-    prm.precond.pmask.resize(rows);
-    for(size_t i = 0; i < rows; ++i) prm.precond.pmask[i] = (i >= nu);
-
     // Initialize the solver with the system matrix.
+    // Use the block_matrix adapter to convert the matrix into
+    // the block format on the fly:
     prof.tic("setup");
-    Solver solve(A, prm);
+    auto Ab = amgcl::adapter::block_matrix<dmat_type>(A);
+    Solver solve(Ab);
     prof.toc("setup");
 
     // Show the mini-report on the constructed solver:
@@ -95,8 +94,15 @@ int main(int argc, char *argv[]) {
     int iters;
     double error;
     std::vector<double> x(rows, 0.0);
+
+    // Reinterpret both the RHS and the solution vectors as block-valued:
+    auto f_ptr = reinterpret_cast<dvec_type*>(f.data());
+    auto x_ptr = reinterpret_cast<dvec_type*>(x.data());
+    auto F = amgcl::make_iterator_range(f_ptr, f_ptr + rows / 4);
+    auto X = amgcl::make_iterator_range(x_ptr, x_ptr + rows / 4);
+
     prof.tic("solve");
-    std::tie(iters, error) = solve(A, f, x);
+    std::tie(iters, error) = solve(Ab, F, X);
     prof.toc("solve");
 
     // Output the number of iterations, the relative error,

tutorial/3.CoupCons3D/coupcons3d_spc.cpp

@@ -0,0 +1,107 @@
+#include <iostream>
+#include <string>
+
+#include <amgcl/backend/builtin.hpp>
+#include <amgcl/adapter/crs_tuple.hpp>
+#include <amgcl/value_type/static_matrix.hpp>
+#include <amgcl/adapter/block_matrix.hpp>
+#include <amgcl/preconditioner/schur_pressure_correction.hpp>
+#include <amgcl/make_solver.hpp>
+#include <amgcl/make_block_solver.hpp>
+#include <amgcl/amg.hpp>
+#include <amgcl/solver/bicgstab.hpp>
+#include <amgcl/solver/preonly.hpp>
+#include <amgcl/coarsening/aggregation.hpp>
+#include <amgcl/relaxation/ilu0.hpp>
+#include <amgcl/relaxation/spai0.hpp>
+#include <amgcl/relaxation/as_preconditioner.hpp>
+
+#include <amgcl/io/mm.hpp>
+#include <amgcl/profiler.hpp>
+
+//---------------------------------------------------------------------------
+int main(int argc, char *argv[]) {
+    // The matrix and the RHS file names should be in the command line options:
+    if (argc < 3) {
+        std::cerr << "Usage: " << argv[0] << " <matrix.mtx> <nu>" << std::endl;
+        return 1;
+    }
+
+    // The profiler:
+    amgcl::profiler<> prof("CoupCons3D");
+
+    // Read the system matrix:
+    size_t rows, cols;
+    std::vector<ptrdiff_t> ptr, col;
+    std::vector<double> val;
+
+    prof.tic("read");
+    std::tie(rows, cols) = amgcl::io::mm_reader(argv[1])(ptr, col, val);
+    std::cout << "Matrix " << argv[1] << ": " << rows << "x" << cols << std::endl;
+    prof.toc("read");
+
+    // The RHS is filled with ones:
+    std::vector<double> f(rows, 1.0);
+
+    // The number of unknowns in the U subsystem
+    size_t nu = std::stoi(argv[2]);
+
+    // We use the tuple of CRS arrays to represent the system matrix.
+    // Note that std::tie creates a tuple of references, so no data is actually
+    // copied here:
+    auto A = std::tie(rows, ptr, col, val);
+
+    // Compose the solver type
+    typedef amgcl::backend::builtin<double> SBackend; // the outer iterative solver backend
+    typedef amgcl::backend::builtin<float> PBackend;  // the PSolver backend
+    typedef amgcl::backend::builtin<
+        amgcl::static_matrix<float,4,4>> UBackend;    // the USolver backend
+
+    typedef amgcl::make_solver<
+        amgcl::preconditioner::schur_pressure_correction<
+            amgcl::make_block_solver<
+                amgcl::amg<
+                    UBackend,
+                    amgcl::coarsening::aggregation,
+                    amgcl::relaxation::ilu0
+                    >,
+                amgcl::solver::preonly<UBackend>
+                >,
+            amgcl::make_solver<
+                amgcl::relaxation::as_preconditioner<
+                    PBackend,
+                    amgcl::relaxation::spai0
+                    >,
+                amgcl::solver::preonly<PBackend>
+                >
+            >,
+        amgcl::solver::bicgstab<SBackend>
+        > Solver;
+
+    // Solver parameters
+    Solver::params prm;
+    prm.precond.pmask.resize(rows);
+    for(size_t i = 0; i < rows; ++i) prm.precond.pmask[i] = (i >= nu);
+
+    // Initialize the solver with the system matrix.
+    prof.tic("setup");
+    Solver solve(A, prm);
+    prof.toc("setup");
+
+    // Show the mini-report on the constructed solver:
+    std::cout << solve << std::endl;
+
+    // Solve the system with the zero initial approximation:
+    int iters;
+    double error;
+    std::vector<double> x(rows, 0.0);
+    prof.tic("solve");
+    std::tie(iters, error) = solve(A, f, x);
+    prof.toc("solve");
+
+    // Output the number of iterations, the relative error,
+    // and the profiling data:
+    std::cout << "Iters: " << iters << std::endl
+              << "Error: " << error << std::endl
+              << prof << std::endl;
+}