Poisson in a cube

cube.cpp

#include <iostream>
#include <iomanip>
#include <vector>
#include <vexcl/vexcl.hpp>

#include <amgcl/amgcl.hpp>
#include <amgcl/interp_smoothed_aggr.hpp>
#include <amgcl/aggr_plain.hpp>
#include <amgcl/level_vexcl.hpp>
#include <amgcl/gmres.hpp>

const double bot = 1;
const double top = 2;

int main() {
    const int n = 100;
    const int nx = n;
    const int ny = n;
    const int nz = n;
    const int n3 = nx * ny * nz;

    std::vector<int>    row;
    std::vector<int>    col;
    std::vector<double> val;
    std::vector<double> rhs;

    row.reserve(n3 + 1);
    col.reserve(n3 * 7);
    val.reserve(n3 * 7);
    rhs.reserve(n3);

    row.push_back(0);

    auto append = [&](int c, double v) {
        col.push_back(c);
        val.push_back(v);
    };

    for(size_t k = 0, idx = 0; k < nz; ++k) {
        for(size_t j = 0; j < ny; ++j) {
            for(size_t i = 0; i < nx; ++i, ++idx) {
                if (k == 0) {
                    rhs.push_back(bot);
                    append(idx, 1);
                } else if (k + 1 == nz) {
                    rhs.push_back(top);
                    append(idx, 1);
                } else {
                    rhs.push_back(0);
                    append(idx, 6);

                    if (i == 0) {
                        append(idx + 1, -2);
                    } else if (i + 1 == nx) {
                        append(idx - 1, -2);
                    } else {
                        append(idx - 1, -1);
                        append(idx + 1, -1);
                    }

                    if (j == 0) {
                        append(idx + nx, -2);
                    } else if (j + 1 == ny) {
                        append(idx - nx, -2);
                    } else {
                        append(idx - nx, -1);
                        append(idx + nx, -1);
                    }

                    append(idx - nx * ny, -1);
                    append(idx + nx * ny, -1);
                }

                row.push_back(col.size());
            }
        }
    }

    vex::Context ctx(vex::Filter::Env);
    std::cout << ctx << std::endl;

    typedef amgcl::solver<
        double /* matrix value type */, int /* matrix index type */,
        amgcl::interp::smoothed_aggregation<amgcl::aggr::plain>,
        amgcl::level::vexcl<amgcl::relax::spai0>
    > AMG;

    AMG::params prm;
    prm.level.ctx = &ctx;

    AMG amg(amgcl::sparse::map(n3, n3, row.data(), col.data(), val.data()), prm);

    std::cout << amg << std::endl;

    vex::vector<double> x(ctx, n3);
    x = 1;

    std::cout << "err = " <<
        std::get<1>(
                amgcl::solve(
                    amg.top_matrix(), vex::vector<double>(ctx, rhs), amg,
                    x, amgcl::gmres_tag()
                    )
                )
        << std::endl;

    std::vector<double> X(n3);
    vex::copy(x, X);

    std::cout << "x[:][0][0] = {";
    for(size_t k = 0, idx = 0; k < nz; ++k, idx += nx * ny) {
        if (k % 5 == 0) std::cout << "\n" << std::setw(4) << k << ":";
        std::cout << std::scientific << std::setw(14) << X[idx];
    }
    std::cout << "\n}" << std::endl;
}

Makefile

cube: cube.cpp
	g++ -o $@ $^ -std=c++11 -O3 \
	    -I$(HOME)/work/vexcl -I$(HOME)/work/amgcl \
	    -lOpenCL -lboost_system

output.txt

$ ./cube 
1. Tesla K20c (NVIDIA CUDA)

Number of levels:    5
Operator complexity: 1.61
Grid complexity:     1.13

level     unknowns       nonzeros
---------------------------------
    0      1000000        6840800 (62.11%)
    1       118916        3644604 (33.09%)
    2         7389         506507 ( 4.60%)
    3          366          21476 ( 0.19%)
    4           21            423 ( 0.00%)

err = 6.78417e-09
x[:][0][0] = {
   0:  1.000000e+00  1.010101e+00  1.020202e+00  1.030303e+00  1.040404e+00
   5:  1.050505e+00  1.060606e+00  1.070707e+00  1.080808e+00  1.090909e+00
  10:  1.101010e+00  1.111111e+00  1.121212e+00  1.131313e+00  1.141414e+00
  15:  1.151515e+00  1.161616e+00  1.171717e+00  1.181818e+00  1.191919e+00
  20:  1.202020e+00  1.212121e+00  1.222222e+00  1.232323e+00  1.242424e+00
  25:  1.252525e+00  1.262626e+00  1.272727e+00  1.282828e+00  1.292929e+00
  30:  1.303030e+00  1.313131e+00  1.323232e+00  1.333333e+00  1.343434e+00
  35:  1.353535e+00  1.363636e+00  1.373737e+00  1.383838e+00  1.393939e+00
  40:  1.404040e+00  1.414141e+00  1.424242e+00  1.434343e+00  1.444444e+00
  45:  1.454545e+00  1.464646e+00  1.474747e+00  1.484848e+00  1.494949e+00
  50:  1.505051e+00  1.515152e+00  1.525253e+00  1.535354e+00  1.545455e+00
  55:  1.555556e+00  1.565657e+00  1.575758e+00  1.585859e+00  1.595960e+00
  60:  1.606061e+00  1.616162e+00  1.626263e+00  1.636364e+00  1.646465e+00
  65:  1.656566e+00  1.666667e+00  1.676768e+00  1.686869e+00  1.696970e+00
  70:  1.707071e+00  1.717172e+00  1.727273e+00  1.737374e+00  1.747475e+00
  75:  1.757576e+00  1.767677e+00  1.777778e+00  1.787879e+00  1.797980e+00
  80:  1.808081e+00  1.818182e+00  1.828283e+00  1.838384e+00  1.848485e+00
  85:  1.858586e+00  1.868687e+00  1.878788e+00  1.888889e+00  1.898990e+00
  90:  1.909091e+00  1.919192e+00  1.929293e+00  1.939394e+00  1.949495e+00
  95:  1.959596e+00  1.969697e+00  1.979798e+00  1.989899e+00  2.000000e+00
}