version 7.4.2 (electronic-structure#809)

Source code:
* introduce r3:: namespace for 3D vectors and matrices
* introduce fft:: namespace for the FFT related functions
* introduce mpi:: namaspace for the MPI related functions
* introduce la:: namespace for linear algebra related functions

Containers:
* all docker files are stored in ./dockerfile folder
* added base container for OpenMPI + gcc-12

Fixes:
* use conversion operator semantics to wrap MPI types
* mdarray wrapper for pybind11


Co-authored-by: Simon Pintarelli <1237199+simonpintarelli@users.noreply.github.com>

.github/workflows/build.yml

@@ -36,6 +36,26 @@ jobs:
           cd ${GITHUB_WORKSPACE}/build
           spack --color always -e gcc-build-env build-env $SPEC_GCC -- ctest --output-on-failure -L cpu_serial
 
+  build_openmpi:
+    runs-on: ubuntu-latest
+    container: electronicstructure/sirius_openmpi
+    env:
+      SPEC_GCC: "sirius@develop %gcc@12: +python +tests +apps +vdwxc +scalapack +fortran build_type=RelWithDebInfo ^openblas ^openmpi"
+    steps:
+      - uses: actions/checkout@v2
+      - name: Show the spec
+        run: spack --color always -e gcc-build-env spec -I $SPEC_GCC
+      - name: Configure SIRIUS
+        run: |
+          cd ${GITHUB_WORKSPACE}
+          mkdir build
+          cd build
+          spack --color always -e gcc-build-env build-env $SPEC_GCC -- cmake .. -DUSE_SCALAPACK=1 -DUSE_VDWXC=1 -DBUILD_TESTING=1 -DCREATE_FORTRAN_BINDINGS=1 -DCREATE_PYTHON_MODULE=1 -DCMAKE_BUILD_TYPE=RelWithDebInfo
+      - name: Build SIRIUS
+        run: |
+          cd ${GITHUB_WORKSPACE}/build
+          spack --color always -e gcc-build-env build-env $SPEC_GCC -- make
+
   build_elpa:
     runs-on: ubuntu-latest
     container: electronicstructure/sirius

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.18)
 
-project(SIRIUS VERSION 7.4.1)
+project(SIRIUS VERSION 7.4.2)
 
 # user variables
 set(CREATE_PYTHON_MODULE    OFF CACHE BOOL "create sirius Python module")

apps/atoms/atom.cpp

@@ -677,7 +677,7 @@ void generate_atom_file(Free_atom&         a,
     sirius::Atom_symmetry_class atom_class(0, a);
     atom_class.set_spherical_potential(veff);
     atom_class.generate_radial_functions(relativity_t::none);
-    pstdout pout(Communicator::self());
+    mpi::pstdout pout(mpi::Communicator::self());
     atom_class.write_enu(pout);
     std::cout << pout.flush(0);
 

apps/dft_loop/sirius.scf.cpp

@@ -25,7 +25,7 @@ struct task_t
 void json_output_common(json& dict__)
 {
     dict__["git_hash"] = sirius::git_hash();
-    dict__["comm_world_size"] = sddk::Communicator::world().size();
+    dict__["comm_world_size"] = mpi::Communicator::world().size();
     dict__["threads_per_rank"] = omp_get_max_threads();
 }
 
@@ -72,7 +72,7 @@ create_sim_ctx(std::string fname__, cmd_args const& args__)
 {
     auto json = preprocess_json_input(fname__);
 
-    auto ctx_ptr = std::make_unique<Simulation_context>(json.dump(), Communicator::world());
+    auto ctx_ptr = std::make_unique<Simulation_context>(json.dump(), mpi::Communicator::world());
     Simulation_context& ctx = *ctx_ptr;
 
     auto& inp = ctx.cfg().parameters();
@@ -360,7 +360,7 @@ void run_tasks(cmd_args const& args)
     }
 
     if (!fs::exists(fpath)) {
-        if (Communicator::world().rank() == 0) {
+        if (mpi::Communicator::world().rank() == 0) {
             std::printf("input file does not exist\n");
         }
         return;
@@ -508,7 +508,7 @@ void run_tasks(cmd_args const& args)
         band.solve<double, double>(ks, H0, ctx->cfg().iterative_solver().energy_tolerance());
 
         ks.sync_band<double, sync_band_t::energy>();
-        if (Communicator::world().rank() == 0) {
+        if (mpi::Communicator::world().rank() == 0) {
             json dict;
             dict["header"] = {};
             dict["header"]["x_axis"] = x_axis;
@@ -587,7 +587,7 @@ int main(int argn, char** argv)
 
     run_tasks(args);
 
-    int my_rank = Communicator::world().rank();
+    int my_rank = mpi::Communicator::world().rank();
 
     sirius::finalize(1);
 

apps/nlcg/sirius.nlcg.cpp

@@ -18,14 +18,14 @@ void json_output_common(json& dict__)
 {
     dict__["git_hash"] = sirius::git_hash();
     //dict__["build_date"] = build_date;
-    dict__["comm_world_size"] = sddk::Communicator::world().size();
+    dict__["comm_world_size"] = mpi::Communicator::world().size();
     dict__["threads_per_rank"] = omp_get_max_threads();
 }
 
 std::unique_ptr<Simulation_context> create_sim_ctx(std::string     fname__,
                                                    cmd_args const& args__)
 {
-    auto ctx_ptr = std::unique_ptr<Simulation_context>(new Simulation_context(fname__, sddk::Communicator::world()));
+    auto ctx_ptr = std::unique_ptr<Simulation_context>(new Simulation_context(fname__, mpi::Communicator::world()));
     Simulation_context& ctx = *ctx_ptr;
 
     auto& inp = ctx.cfg().parameters();
@@ -247,7 +247,7 @@ void run_tasks(cmd_args const& args)
     /* get the input file name */
     std::string fname = args.value<std::string>("input", "sirius.json");
     if (!utils::file_exists(fname)) {
-        if (sddk::Communicator::world().rank() == 0) {
+        if (mpi::Communicator::world().rank() == 0) {
             std::printf("input file does not exist\n");
         }
         return;
@@ -297,7 +297,7 @@ int main(int argn, char** argv)
 
     run_tasks(args);
 
-    int my_rank = sddk::Communicator::world().rank();
+    int my_rank = mpi::Communicator::world().rank();
 
     sirius::finalize(1);
 

apps/tests/test_allgather.cpp

@@ -7,37 +7,37 @@ void test_allgather()
     int N = 11;
     std::vector<double> vec(N, 0.0);
 
-    sddk::splindex<sddk::splindex_t::block> spl(N, sddk::Communicator::world().size(), sddk::Communicator::world().rank());
+    sddk::splindex<sddk::splindex_t::block> spl(N, mpi::Communicator::world().size(), mpi::Communicator::world().rank());
 
     for (int i = 0; i < spl.local_size(); i++) {
-        vec[spl[i]] = sddk::Communicator::world().rank() + 1.0;
+        vec[spl[i]] = mpi::Communicator::world().rank() + 1.0;
     }
 
     {
-        sddk::pstdout pout(sddk::Communicator::world());
-        if (sddk::Communicator::world().rank() == 0) {
+        mpi::pstdout pout(mpi::Communicator::world());
+        if (mpi::Communicator::world().rank() == 0) {
             pout << "before" << std::endl;
         }
-        pout << "rank : " << sddk::Communicator::world().rank() << " array : ";
+        pout << "rank : " << mpi::Communicator::world().rank() << " array : ";
         for (int i = 0; i < N; i++) {
             pout << vec[i];
         }
         pout << std::endl;
         std::cout << pout.flush(0);
 
-        sddk::Communicator::world().allgather(&vec[0], spl.local_size(), spl.global_offset());
+        mpi::Communicator::world().allgather(&vec[0], spl.local_size(), spl.global_offset());
 
-        if (sddk::Communicator::world().rank() == 0) {
+        if (mpi::Communicator::world().rank() == 0) {
             pout << "after" << std::endl;
         }
-        pout << "rank : " << sddk::Communicator::world().rank() << " array : ";
+        pout << "rank : " << mpi::Communicator::world().rank() << " array : ";
         for (int i = 0; i < N; i++) {
             pout << vec[i];
         }
         pout << std::endl;
         std::cout << pout.flush(0);
     }
-    sddk::Communicator::world().barrier();
+    mpi::Communicator::world().barrier();
 }
 
 int main(int argn, char** argv)

apps/tests/test_bcast_v2.cpp

@@ -9,11 +9,11 @@ int test1(int size)
     sddk::mdarray<char, 1> buf(size);
     buf.zero();
 
-    for (int r = 0; r < Communicator::world().size(); r++) {
-        Communicator::world().bcast(buf.at(memory_t::host), size, r);
+    for (int r = 0; r < mpi::Communicator::world().size(); r++) {
+        mpi::Communicator::world().bcast(buf.at(memory_t::host), size, r);
     }
     t += utils::wtime();
-    if (Communicator::world().rank() == 0) {
+    if (mpi::Communicator::world().rank() == 0) {
         printf("time : %f sec.\n", t);
     }
     return 0;

apps/tests/test_blacs.cpp

@@ -5,11 +5,11 @@ int main(int argn, char** argv)
     sirius::initialize(true);
 
 #if defined(SIRIUS_SCALAPACK)
-    std::cout << sddk::Communicator::self().size() << " " << sddk::Communicator::self().rank() << std::endl;
-    std::cout << sddk::Communicator::world().size() << " " << sddk::Communicator::world().rank() << std::endl;
+    std::cout << mpi::Communicator::self().size() << " " << mpi::Communicator::self().rank() << std::endl;
+    std::cout << mpi::Communicator::world().size() << " " << mpi::Communicator::world().rank() << std::endl;
 
-    auto blacs_handler = sddk::linalg_base::create_blacs_handler(sddk::Communicator::self().mpi_comm());
-    blacs_handler = sddk::linalg_base::create_blacs_handler(sddk::Communicator::world().mpi_comm());
+    auto blacs_handler = la::linalg_base::create_blacs_handler(mpi::Communicator::self().native());
+    blacs_handler = la::linalg_base::create_blacs_handler(mpi::Communicator::world().native());
     std::cout << blacs_handler << std::endl;
 
     sirius::finalize(true);

apps/tests/test_comm_split.cpp

@@ -5,19 +5,19 @@ using namespace sddk;
 
 void test_comm_split(int comm_size)
 {
-    if (Communicator::world().rank() == 0) {
+    if (mpi::Communicator::world().rank() == 0) {
         printf("sub comm size: %i\n", comm_size);
     }
 
-    auto c1 = Communicator::world().split(Communicator::world().rank() / comm_size);
-    auto c2 = Communicator::world().split(Communicator::world().rank() % comm_size);
+    auto c1 = mpi::Communicator::world().split(mpi::Communicator::world().rank() / comm_size);
+    auto c2 = mpi::Communicator::world().split(mpi::Communicator::world().rank() % comm_size);
 
     auto c3 = c1;
-    Communicator c4(c2);
+    mpi::Communicator c4(c2);
 
-    pstdout pout(Communicator::world());
+    mpi::pstdout pout(mpi::Communicator::world());
 
-    pout << "global rank: " << Communicator::world().rank() << ", c1.rank: " << c1.rank()
+    pout << "global rank: " << mpi::Communicator::world().rank() << ", c1.rank: " << c1.rank()
          << ", c2.rank: " << c2.rank() << std::endl;
     std::cout << pout.flush(0);
 }

apps/tests/test_davidson.cpp

@@ -70,7 +70,7 @@ diagonalize(Simulation_context& ctx__, std::array<double, 3> vk__, Potential& po
             wf::num_mag_dims(ctx__.num_mag_dims()), kp.spinor_wave_functions(), [&](int i, int ispn){return eval_tol__;}, res_tol__,
             60, locking, subspace_size__, estimate_eval__, extra_ortho__, std::cout, 2);
 
-    if (sddk::Communicator::world().rank() == 0 && only_kin__) {
+    if (mpi::Communicator::world().rank() == 0 && only_kin__) {
         std::vector<double> ekin(kp.num_gkvec());
         for (int i = 0; i < kp.num_gkvec(); i++) {
             ekin[i] = 0.5 * kp.gkvec().template gkvec_cart<sddk::index_domain_t::global>(i).length2();
@@ -85,7 +85,7 @@ diagonalize(Simulation_context& ctx__, std::array<double, 3> vk__, Potential& po
         printf("maximum eigen-value difference: %20.16e\n", max_diff);
     }
 
-    if (sddk::Communicator::world().rank() == 0 && !only_kin__) {
+    if (mpi::Communicator::world().rank() == 0 && !only_kin__) {
         std::cout << "Converged eigen-values" << std::endl;
         for (int i = 0; i < ctx__.num_bands(); i++) {
             printf("e[%i] = %20.16f\n", i, result.eval(i, 0));
@@ -212,7 +212,7 @@ int main(int argn, char** argv)
 
     sirius::initialize(1);
     test_davidson(args);
-    int rank = sddk::Communicator::world().rank();
+    int rank = mpi::Communicator::world().rank();
     sirius::finalize();
     if (rank == 0)  {
         const auto timing_result = ::utils::global_rtgraph_timer.process();

apps/tests/test_eigen.cpp

@@ -5,20 +5,20 @@ using namespace sirius;
 
 template <typename T>
 double
-test_diag(sddk::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int bs__, bool test_gen__, std::string name__,
-          Eigensolver& solver)
+test_diag(la::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int bs__, bool test_gen__, std::string name__,
+          la::Eigensolver& solver)
 {
     auto A_ref = random_symmetric<T>(N__, bs__, blacs_grid__);
-    sddk::dmatrix<T> A(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
+    la::dmatrix<T> A(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
     A_ref >> A;
 
-    sddk::dmatrix<T> Z(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
+    la::dmatrix<T> Z(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
 
-    sddk::dmatrix<T> B;
-    sddk::dmatrix<T> B_ref;
+    la::dmatrix<T> B;
+    la::dmatrix<T> B_ref;
     if (test_gen__) {
         B_ref = random_positive_definite<T>(N__, bs__, blacs_grid__);
-        B = sddk::dmatrix<T>(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
+        B = la::dmatrix<T>(N__, N__, blacs_grid__, bs__, bs__, solver.host_memory_t());
         B_ref >> B;
     }
 
@@ -92,22 +92,22 @@ test_diag(sddk::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int
     if (test_gen__) {
         /* lambda * B * Z */
 #if defined(SIRIUS_SCALAPACK)
-        sddk::linalg(sddk::linalg_t::scalapack).gemm('N', 'N', n__, nev__, n__, &sddk::linalg_const<T>::one(),
-                B_ref, 0, 0, A, 0, 0, &sddk::linalg_const<T>::zero(), B, 0, 0);
+        la::wrap(la::lib_t::scalapack).gemm('N', 'N', n__, nev__, n__, &la::constant<T>::one(),
+                B_ref, 0, 0, A, 0, 0, &la::constant<T>::zero(), B, 0, 0);
 #else
-        sddk::linalg(sddk::linalg_t::blas).gemm('N', 'N', n__, nev__, n__, &sddk::linalg_const<T>::one(),
-                &B_ref(0, 0), B_ref.ld(), &A(0, 0), A.ld(), &sddk::linalg_const<T>::zero(), &B(0, 0), B.ld());
+        la::wrap(la::lib_t::blas).gemm('N', 'N', n__, nev__, n__, &la::constant<T>::one(),
+                &B_ref(0, 0), B_ref.ld(), &A(0, 0), A.ld(), &la::constant<T>::zero(), &B(0, 0), B.ld());
 #endif
         B >> A;
     }
 
     /* A * Z - lambda * B * Z */
 #if defined(SIRIUS_SCALAPACK)
-    sddk::linalg(sddk::linalg_t::scalapack).gemm('N', 'N', n__, nev__, n__, &sddk::linalg_const<T>::one(), A_ref, 0, 0, Z, 0, 0,
-        &sddk::linalg_const<T>::m_one(), A, 0, 0);
+    la::wrap(la::lib_t::scalapack).gemm('N', 'N', n__, nev__, n__, &la::constant<T>::one(), A_ref, 0, 0, Z, 0, 0,
+        &la::constant<T>::m_one(), A, 0, 0);
 #else
-    sddk::linalg(sddk::linalg_t::blas).gemm('N', 'N', n__, nev__, n__, &sddk::linalg_const<T>::one(), &A_ref(0, 0), A_ref.ld(),
-            &Z(0, 0), Z.ld(), &sddk::linalg_const<T>::m_one(), &A(0, 0), A.ld());
+    la::wrap(la::lib_t::blas).gemm('N', 'N', n__, nev__, n__, &la::constant<T>::one(), &A_ref(0, 0), A_ref.ld(),
+            &Z(0, 0), Z.ld(), &la::constant<T>::m_one(), &A(0, 0), A.ld());
 #endif
     double diff{0};
     for (int j = 0; j < A.num_cols_local(); j++) {
@@ -117,7 +117,7 @@ test_diag(sddk::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int
             }
         }
     }
-    blacs_grid__.comm().template allreduce<double, sddk::mpi_op_t::max>(&diff, 1);
+    blacs_grid__.comm().template allreduce<double, mpi::op_t::max>(&diff, 1);
     if (blacs_grid__.comm().rank() == 0) {
         printf("maximum difference: %22.18f\n", diff);
     }
@@ -131,12 +131,12 @@ test_diag(sddk::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int
     return t;
 }
 
-void test_diag2(sddk::BLACS_grid const& blacs_grid__,
+void test_diag2(la::BLACS_grid const& blacs_grid__,
                 int bs__,
                 std::string name__,
                 std::string fname__)
 {
-    auto solver = Eigensolver_factory(name__);
+    auto solver = la::Eigensolver_factory(name__);
 
     sddk::matrix<std::complex<double>> full_mtrx;
     int n;
@@ -162,8 +162,8 @@ void test_diag2(sddk::BLACS_grid const& blacs_grid__,
     }
 
     std::vector<double> eval(n);
-    sddk::dmatrix<std::complex<double>> A(n, n, blacs_grid__, bs__, bs__);
-    sddk::dmatrix<std::complex<double>> Z(n, n, blacs_grid__, bs__, bs__);
+    la::dmatrix<std::complex<double>> A(n, n, blacs_grid__, bs__, bs__);
+    la::dmatrix<std::complex<double>> Z(n, n, blacs_grid__, bs__, bs__);
 
     for (int j = 0; j < A.num_cols_local(); j++) {
         for (int i = 0; i < A.num_rows_local(); i++) {
@@ -190,8 +190,8 @@ void call_test(std::vector<int> mpi_grid__,
                int repeat__,
                int type__)
 {
-    auto solver = Eigensolver_factory(name__);
-    sddk::BLACS_grid blacs_grid(sddk::Communicator::world(), mpi_grid__[0], mpi_grid__[1]);
+    auto solver = la::Eigensolver_factory(name__);
+    la::BLACS_grid blacs_grid(mpi::Communicator::world(), mpi_grid__[0], mpi_grid__[1]);
     if (fname__.length() == 0) {
         Measurement m;
         for (int i = 0; i < repeat__; i++) {