A skeleton of a C API for the new interface

Makefile.am

@@ -425,6 +425,10 @@ noinst_PROGRAMS = \
   test_real_double_2stage \
   test_complex_double_1stage \
   test_complex_double_2stage \
+  test_c_real_double_1stage \
+  test_c_real_double_2stage \
+  test_c_complex_double_1stage \
+  test_c_complex_double_2stage \
   double_instance@SUFFIX@ \
   real_2stage_banded@SUFFIX@ \
   complex_2stage_banded@SUFFIX@ \
@@ -437,6 +441,8 @@ if WANT_SINGLE_PRECISION_REAL
 noinst_PROGRAMS += \
   test_single_real_1stage \
   test_single_real_2stage \
+  test_c_single_real_1stage \
+  test_c_single_real_2stage \
   single_real_2stage_banded@SUFFIX@ \
   single_real_2stage@SUFFIX@ \
   single_real_1stage@SUFFIX@
@@ -446,6 +452,8 @@ if WANT_SINGLE_PRECISION_COMPLEX
 noinst_PROGRAMS += \
   test_single_complex_1stage \
   test_single_complex_2stage \
+  test_c_single_complex_1stage \
+  test_c_single_complex_2stage \
   single_complex_2stage_banded@SUFFIX@ \
   single_complex_2stage@SUFFIX@ \
   single_complex_1stage@SUFFIX@
@@ -647,6 +655,74 @@ test_single_complex_2stage_FCFLAGS = $(AM_FCFLAGS) $(FC_MODOUT)private_modules $
         -DTEST_GPU=0
 endif
 
+test_c_real_double_1stage_SOURCES = test/C/test.c
+test_c_real_double_1stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_real_double_1stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_REAL \
+        -DTEST_DOUBLE \
+        -DTEST_SOLVER_1STAGE \
+        -DTEST_GPU=0
+
+test_c_real_double_2stage_SOURCES = test/C/test.c
+test_c_real_double_2stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_real_double_2stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_REAL \
+        -DTEST_DOUBLE \
+        -DTEST_SOLVER_2STAGE \
+        -DTEST_GPU=0
+
+test_c_complex_double_1stage_SOURCES = test/C/test.c
+test_c_complex_double_1stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_complex_double_1stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_COMPLEX \
+        -DTEST_DOUBLE \
+        -DTEST_SOLVER_1STAGE \
+        -DTEST_GPU=0
+
+test_c_complex_double_2stage_SOURCES = test/C/test.c
+test_c_complex_double_2stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_complex_double_2stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_COMPLEX \
+        -DTEST_DOUBLE \
+        -DTEST_SOLVER_2STAGE \
+        -DTEST_GPU=0
+
+if WANT_SINGLE_PRECISION_REAL
+test_c_single_real_1stage_SOURCES = test/C/test.c
+test_c_single_real_1stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_single_real_1stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_REAL \
+        -DTEST_SINGLE \
+        -DTEST_SOLVER_1STAGE \
+        -DTEST_GPU=0
+
+test_c_single_real_2stage_SOURCES = test/C/test.c
+test_c_single_real_2stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_single_real_2stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_REAL \
+        -DTEST_SINGLE \
+        -DTEST_SOLVER_2STAGE \
+        -DTEST_GPU=0
+endif
+
+if WANT_SINGLE_PRECISION_COMPLEX
+test_c_single_complex_1stage_SOURCES = test/C/test.c
+test_c_single_complex_1stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_single_complex_1stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_COMPLEX \
+        -DTEST_SINGLE \
+        -DTEST_SOLVER_1STAGE \
+        -DTEST_GPU=0
+
+test_c_single_complex_2stage_SOURCES = test/C/test.c
+test_c_single_complex_2stage_LDADD = $(build_lib) $(FCLIBS)
+test_c_single_complex_2stage_CFLAGS = $(AM_CFLAGS) \
+        -DTEST_COMPLEX \
+        -DTEST_SINGLE \
+        -DTEST_SOLVER_2STAGE \
+        -DTEST_GPU=0
+endif
+
 double_instance@SUFFIX@_SOURCES = test/Fortran/elpa2/double_instance.F90
 double_instance@SUFFIX@_LDADD = $(build_lib)
 double_instance@SUFFIX@_FCFLAGS = $(AM_FCFLAGS) $(FC_MODOUT)private_modules $(FC_MODINC)private_modules
@@ -1051,6 +1127,10 @@ check_SCRIPTS = \
   test_real_double_2stage.sh \
   test_complex_double_1stage.sh \
   test_complex_double_2stage.sh \
+  test_c_real_double_1stage.sh \
+  test_c_real_double_2stage.sh \
+  test_c_complex_double_1stage.sh \
+  test_c_complex_double_2stage.sh \
   double_instance@SUFFIX@.sh \
   real_2stage_banded@SUFFIX@.sh \
   complex_2stage_banded@SUFFIX@.sh \
@@ -1063,6 +1143,8 @@ if WANT_SINGLE_PRECISION_REAL
 check_SCRIPTS += \
   test_single_real_1stage.sh \
   test_single_real_2stage.sh \
+  test_c_single_real_1stage.sh \
+  test_c_single_real_2stage.sh \
   single_real_2stage_banded@SUFFIX@.sh \
   single_real_2stage@SUFFIX@.sh \
   single_real_1stage@SUFFIX@.sh
@@ -1072,6 +1154,8 @@ if WANT_SINGLE_PRECISION_COMPLEX
 check_SCRIPTS += \
   test_single_complex_1stage.sh \
   test_single_complex_2stage.sh \
+  test_c_single_complex_1stage.sh \
+  test_c_single_complex_2stage.sh \
   single_complex_2stage_banded@SUFFIX@.sh \
   single_complex_2stage@SUFFIX@.sh \
   single_complex_1stage@SUFFIX@.sh

elpa/elpa.h

@@ -2,9 +2,14 @@
 #define ELPA_H
 
 #include <limits.h>
+#include <complex.h>
+
+struct elpa_struct;
+typedef struct elpa_struct *elpa_t;
 
 #include <elpa/elpa_constants.h>
 #include <elpa/elpa_generated.h>
+#include <elpa/elpa_generic.h>
 
 const char *elpa_strerr(int elpa_error);
 

elpa/elpa_generic.h

+#pragma once
+
+/**
+ * \todo document elpa_set()
+ */
+#define elpa_set(e, name, value, error) _Generic((value), \
+                int: \
+                  elpa_set_integer, \
+                \
+                double: \
+                  elpa_set_double \
+        )(e, name, value, error)
+
+/**
+ * \todo document elpa_solve()
+ */
+#define elpa_solve(handle, a, ev, q, error) _Generic((a), \
+                double*: \
+                  elpa_solve_real_double, \
+                \
+                float*: \
+                  elpa_solve_real_single, \
+                \
+                double complex*: \
+                  elpa_solve_complex_double, \
+                \
+                float complex*: \
+                  elpa_solve_complex_single \
+        )(handle, a, ev, q, error)

generated_headers.am

@@ -13,9 +13,8 @@ config-f90.h: config.h
 	@echo "Generating $@...";
 	@grep "^#define" $^ > $@ || { rm $@; exit 1; }
 
-elpa/elpa_generated.h: $(top_srcdir)/src/elpa_driver/legacy_interface/elpa_driver_c_interface.F90 \
-                       $(top_srcdir)/src/elpa1/legacy_interface/elpa_1stage_c_interface.F90 \
-                       $(top_srcdir)/src/elpa2/legacy_interface/elpa_2stage_c_interface.F90 | elpa
+elpa/elpa_generated.h: $(top_srcdir)/src/elpa_impl.F90 \
+                       $(top_srcdir)/src/elpa_api.F90 | elpa
 	@rm -f $@
 	$(call extract_interface,!c>)
 

src/elpa_api.F90

@@ -773,9 +773,11 @@ module elpa_api
     !> Parameters
     !> \param   api_version integer, api_version that ELPA should use
     !> \result  error       integer
-    function elpa_init(api_version) result(error)
+    !
+    !c> int elpa_init(int api_version);
+    function elpa_init(api_version) result(error) bind(C, name="elpa_init")
       use elpa_utilities, only : error_unit
-      integer, intent(in) :: api_version
+      integer, intent(in), value :: api_version
       integer             :: error
 
       if (earliest_api_version <= api_version .and. api_version <= current_api_version) then
@@ -800,7 +802,9 @@ module elpa_api
     end function
 
     !> \brief subroutine to uninit the ELPA library. Does nothing at the moment. Might do sth. later
-    subroutine elpa_uninit()
+    !
+    !c> void elpa_uninit(void);
+    subroutine elpa_uninit() bind(C, name="elpa_uninit")
     end subroutine
 
     !> \brief helper function for error strings: NOT public to the user

src/elpa_impl.F90

@@ -156,6 +156,29 @@ module elpa_impl
       endif
     end function
 
+
+    !c> elpa_t elpa_allocate();
+    function elpa_impl_allocate_for_c(error) result(ptr) bind(C, name="elpa_allocate")
+      integer(kind=c_int) :: error
+      type(c_ptr) :: ptr
+      type(elpa_impl_t), pointer :: obj
+
+      obj => elpa_impl_allocate(error)
+      ptr = c_loc(obj)
+    end function
+
+
+    !c> void elpa_deallocate(elpa_t handle);
+    subroutine elpa_impl_deallocate_for_c(handle) bind(C, name="elpa_deallocate")
+      type(c_ptr), value :: handle
+      type(elpa_impl_t), pointer :: self
+
+      call c_f_pointer(handle, self)
+      call self%destroy()
+      deallocate(self)
+    end subroutine
+
+
     !> \brief function to setup an ELPA object and to store the MPI communicators internally
     !> Parameters
     !> \param   self       class(elpa_impl_t), the allocated ELPA object
@@ -166,8 +189,8 @@ module elpa_impl
       integer :: error, error2
       integer :: mpi_comm_rows, mpi_comm_cols, mpierr
 
+#ifdef WITH_MPI
       error = ELPA_ERROR
-
       if (self%is_set("mpi_comm_parent") == 1 .and. &
           self%is_set("process_row") == 1 .and. &
           self%is_set("process_col") == 1) then
@@ -188,6 +211,9 @@ module elpa_impl
       if (self%is_set("mpi_comm_rows") == 1 .and. self%is_set("mpi_comm_cols") == 1) then
         error = ELPA_OK
       endif
+#else
+      error = ELPA_OK
+#endif
 
       if (self%get("timings") == 1) then
         call self%timer%enable()
@@ -195,6 +221,18 @@ module elpa_impl
 
     end function
 
+
+    !c> int elpa_setup(elpa_t handle);
+    function elpa_setup_for_c(handle) result(error) bind(C, name="elpa_setup")
+      type(c_ptr), intent(in), value :: handle
+      type(elpa_impl_t), pointer :: self
+      integer(kind=c_int) :: error
+
+      call c_f_pointer(handle, self)
+      error = self%setup()
+    end function
+
+
     !> \brief subroutine to set an integer key/value pair
     !> Parameters
     !> \param   self       class(elpa_impl_t) the allocated ELPA object
@@ -222,6 +260,22 @@ module elpa_impl
       end if
     end subroutine
 
+
+    !c> void elpa_set_integer(elpa_t handle, const char *name, int value, int *error);
+    subroutine elpa_set_integer_for_c(handle, name_p, value, error) bind(C, name="elpa_set_integer")
+      type(c_ptr), intent(in), value :: handle
+      type(elpa_impl_t), pointer :: self
+      type(c_ptr), intent(in), value :: name_p
+      character(len=elpa_strlen_c(name_p)), pointer :: name
+      integer(kind=c_int), intent(in), value :: value
+      integer(kind=c_int), optional, intent(in) :: error
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(name_p, name)
+      call elpa_set_integer(self, name, value, error)
+    end subroutine
+
+
     !> \brief function to get an integer key/value pair
     !> Parameters
     !> \param   self       class(elpa_impl_t) the allocated ELPA object
@@ -331,6 +385,21 @@ module elpa_impl
     end subroutine
 
 
+    !c> void elpa_set_double(elpa_t handle, const char *name, double value, int *error);
+    subroutine elpa_set_double_for_c(handle, name_p, value, error) bind(C, name="elpa_set_double")
+      type(c_ptr), intent(in), value :: handle
+      type(elpa_impl_t), pointer :: self
+      type(c_ptr), intent(in), value :: name_p
+      character(len=elpa_strlen_c(name_p)), pointer :: name
+      real(kind=c_double), intent(in), value :: value
+      integer(kind=c_int), optional, intent(in) :: error
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(name_p, name)
+      call elpa_set_double(self, name, value, error)
+    end subroutine
+
+
     function elpa_get_double(self, name, error) result(value)
       use iso_c_binding
       use elpa_generated_fortran_interfaces
@@ -426,6 +495,22 @@ module elpa_impl
       endif
     end subroutine
 
+    !c> void elpa_solve_real_double(elpa_t handle, double *a, double *ev, double *q, int *error);
+    subroutine elpa_solve_real_double_for_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_real_double")
+      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
+      integer(kind=c_int), optional, intent(in) :: error
+
+      real(kind=c_double), pointer :: a(:, :), q(:, :), ev(:)
+      type(elpa_impl_t), pointer  :: self
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
+      call c_f_pointer(ev_p, ev, [self%na])
+      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
+
+      call elpa_solve_real_double(self, a, ev, q, error)
+    end subroutine
+
 
     subroutine elpa_solve_real_single(self, a, ev, q, error)
       use elpa2_impl
@@ -473,6 +558,23 @@ module elpa_impl
     end subroutine
 
 
+    !c> void elpa_solve_real_single(elpa_t handle, float *a, float *ev, float *q, int *error);
+    subroutine elpa_solve_real_single_for_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_real_single")
+      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
+      integer(kind=c_int), optional, intent(in) :: error
+
+      real(kind=c_float), pointer :: a(:, :), q(:, :), ev(:)
+      type(elpa_impl_t), pointer  :: self
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
+      call c_f_pointer(ev_p, ev, [self%na])
+      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
+
+      call elpa_solve_real_single(self, a, ev, q, error)
+    end subroutine
+
+
     subroutine elpa_solve_complex_double(self, a, ev, q, error)
       use elpa2_impl
       use elpa1_impl
@@ -514,6 +616,24 @@ module elpa_impl
     end subroutine
 
 
+    !c> void elpa_solve_complex_double(elpa_t handle, double complex *a, double *ev, double complex *q, int *error);
+    subroutine elpa_solve_complex_double_for_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_complex_double")
+      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
+      integer(kind=c_int), optional, intent(in) :: error
+
+      complex(kind=c_double_complex), pointer :: a(:, :), q(:, :)
+      real(kind=c_double), pointer :: ev(:)
+      type(elpa_impl_t), pointer  :: self
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
+      call c_f_pointer(ev_p, ev, [self%na])
+      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
+
+      call elpa_solve_complex_double(self, a, ev, q, error)
+    end subroutine
+
+
     subroutine elpa_solve_complex_single(self, a, ev, q, error)
       use elpa2_impl
       use elpa1_impl
@@ -561,6 +681,24 @@ module elpa_impl
     end subroutine
 
 
+    !c> void elpa_solve_complex_single(elpa_t handle, float complex *a, float *ev, float complex *q, int *error);
+    subroutine elpa_solve_complex_single_for_c(handle, a_p, ev_p, q_p, error) bind(C, name="elpa_solve_complex_single")
+      type(c_ptr), intent(in), value :: handle, a_p, ev_p, q_p
+      integer(kind=c_int), optional, intent(in) :: error
+
+      complex(kind=c_float_complex), pointer :: a(:, :), q(:, :)
+      real(kind=c_float), pointer :: ev(:)
+      type(elpa_impl_t), pointer  :: self
+
+      call c_f_pointer(handle, self)
+      call c_f_pointer(a_p, a, [self%local_nrows, self%local_ncols])
+      call c_f_pointer(ev_p, ev, [self%na])
+      call c_f_pointer(q_p, q, [self%local_nrows, self%local_ncols])
+
+      call elpa_solve_complex_single(self, a, ev, q, error)
+    end subroutine
+
+
     subroutine elpa_multiply_at_b_double (self,uplo_a, uplo_c, na, ncb, a, lda, ldaCols, b, ldb, ldbCols, &
                                           c, ldc, ldcCols, error)
       use iso_c_binding

test/C/test.c

+/*   This file is part of ELPA.
+
+     The ELPA library was originally created by the ELPA consortium,
+     consisting of the following organizations:
+
+     - Max Planck Computing and Data Facility (MPCDF), formerly known as
+       Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
+     - Bergische Universität Wuppertal, Lehrstuhl für angewandte
+       Informatik,
+     - Technische Universität München, Lehrstuhl für Informatik mit
+       Schwerpunkt Wissenschaftliches Rechnen ,
+     - Fritz-Haber-Institut, Berlin, Abt. Theorie,
+     - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
+       Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
+       and
+     - IBM Deutschland GmbH
+
+
+     More information can be found here:
+     http://elpa.mpcdf.mpg.de/
+
+     ELPA is free software: you can redistribute it and/or modify
+     it under the terms of the version 3 of the license of the
+     GNU Lesser General Public License as published by the Free
+     Software Foundation.
+
+     ELPA is distributed in the hope that it will be useful,
+     but WITHOUT ANY WARRANTY; without even the implied warranty of
+     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+     GNU Lesser General Public License for more details.
+
+     You should have received a copy of the GNU Lesser General Public License
+     along with ELPA.  If not, see <http://www.gnu.org/licenses/>
+
+     ELPA reflects a substantial effort on the part of the original
+     ELPA consortium, and we ask you to respect the spirit of the
+     license that we chose: i.e., please contribute any changes you
+     may have back to the original ELPA library distribution, and keep
+     any derivatives of ELPA under the same license that we chose for
+     the original distribution, the GNU Lesser General Public License. */
+#include <stdio.h>
+#include <stdlib.h>
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+#include <math.h>
+
+#include <elpa/elpa.h>
+#include <assert.h>
+
+#include "config.h"
+#include "test/shared/generated.h"
+
+#if !(defined(TEST_REAL) ^ defined(TEST_COMPLEX))
+error: define exactly one of TEST_REAL or TEST_COMPLEX
+#endif
+
+#if !(defined(TEST_SINGLE) ^ defined(TEST_DOUBLE))
+error: define exactly one of TEST_SINGLE or TEST_DOUBLE
+#endif
+
+#if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE))
+error: define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE
+#endif
+
+#ifdef TEST_SINGLE
+#  define EV_TYPE float
+#  ifdef TEST_REAL
+#    define MATRIX_TYPE float
+#  else
+#    define MATRIX_TYPE float complex
+#  endif
+#else
+#  define EV_TYPE double
+#  ifdef TEST_REAL
+#    define MATRIX_TYPE double
+#  else
+#    define MATRIX_TYPE double complex
+#  endif
+#endif
+
+#define assert_elpa_ok(x) assert(x == ELPA_OK)
+
+
+int main(int argc, char** argv) {
+   /* matrix dimensions */
+   int na, nev, nblk;
+
+   /* mpi */
+   int myid, nprocs;
+   int na_cols, na_rows;
+   int np_cols, np_rows;
+   int my_prow, my_pcol;
+   int mpi_comm;
+
+   /* blacs */
+   int my_blacs_ctxt, sc_desc[9], info, nprow, npcol;
+
+   /* The Matrix */
+   MATRIX_TYPE *a, *as, *z;
+   EV_TYPE *ev;
+
+   int error, status;
+
+   elpa_t handle;
+
+
+#ifdef WITH_MPI
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+#else
+   nprocs = 1;
+   myid = 0;
+#endif
+
+   na = 1000;
+   nev = 500;
+   nblk = 16;
+
+   for (np_cols = (int) sqrt((double) nprocs); np_cols > 1; np_cols--) {
+     if (nprocs % np_cols == 0) {
+       break;
+     }
+   }
+
+   np_rows = nprocs/np_cols;
+
+   /* set up blacs */
+   /* convert communicators before */
+#ifdef WITH_MPI
+   mpi_comm = MPI_Comm_c2f(MPI_COMM_WORLD);
+#else
+   mpi_comm = 0;
+#endif
+   set_up_blacsgrid_f(mpi_comm, &my_blacs_ctxt, &np_rows, &np_cols, &nprow, &npcol, &my_prow, &my_pcol);
+   set_up_blacs_descriptor_f(na, nblk, my_prow, my_pcol, np_rows, np_cols, &na_rows, &na_cols, sc_desc, my_blacs_ctxt, &info);
+
+   /* allocate the matrices needed for elpa */
+   a  = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
+   z  = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
+   as = calloc(na_rows*na_cols, sizeof(MATRIX_TYPE));
+   ev = calloc(na, sizeof(EV_TYPE));
+
+#ifdef TEST_REAL
+#ifdef TEST_DOUBLE
+   prepare_matrix_real_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
+#else
+   prepare_matrix_real_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
+#endif
+#else
+#ifdef TEST_DOUBLE
+   prepare_matrix_complex_double_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
+#else
+   prepare_matrix_complex_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
+#endif
+#endif
+
+   if (elpa_init(CURRENT_API_VERSION) != ELPA_OK) {
+     fprintf(stderr, "Error: ELPA API version not supported");
+     exit(1);
+   }
+
+   handle = elpa_allocate(&error);
+   assert_elpa_ok(error);
+
+   /* Set parameters */
+   elpa_set(handle, "na", na, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "nev", nev, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "local_nrows", na_rows, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "local_ncols", na_cols, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "nblk", nblk, &error);
+   assert_elpa_ok(error);
+
+#ifdef WITH_MPI
+   elpa_set(handle, "mpi_comm_parent", MPI_COMM_WORLD, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "process_row", my_prow, &error);
+   assert_elpa_ok(error);
+
+   elpa_set(handle, "process_col", my_pcol, &error);
+   assert_elpa_ok(error);
+#endif
+
+   /* Setup */
+   assert_elpa_ok(elpa_setup(handle));
+
+   /* Set tunables */
+#ifdef TEST_SOLVER_1STAGE
+   elpa_set(handle, "solver", ELPA_SOLVER_1STAGE, &error);
+#else
+   elpa_set(handle, "solver", ELPA_SOLVER_2STAGE, &error);
+#endif
+   assert_elpa_ok(error);
+