Playing with an interface

Makefile.am

@@ -17,6 +17,8 @@ libelpa@SUFFIX@_public_la_SOURCES = \
 	src/elpa1_auxiliary.F90 \
         src/elpa1_utilities.F90 \
         src/elpa2_utilities.F90 \
+	src/elpa_init.F90 \
+	src/elpa_t.F90 \
         src/elpa_utilities.F90
 
 # internal parts
@@ -322,6 +324,7 @@ dist_files_DATA = \
   test/Fortran/test_invert_trm_real.F90 \
   test/Fortran/test_cholesky_complex.F90 \
   test/Fortran/test_invert_trm_complex.F90 \
+  test/Fortran/test_new_interface.F90 \
   src/elpa2_print_kernels.F90
 
 #end needed
@@ -363,6 +366,7 @@ noinst_PROGRAMS = \
   elpa2_test_real_c_version@SUFFIX@ \
   elpa2_test_complex_c_version@SUFFIX@ \
   elpa_driver_real_c_version@SUFFIX@ \
+  elpa_test_new_interface@SUFFIX@ \
   elpa_driver_complex_c_version@SUFFIX@
 
 if WANT_SINGLE_PRECISION_COMPLEX
@@ -435,6 +439,11 @@ libelpatest@SUFFIX@_la_SOURCES += \
 	test/shared/redir.c \
 	test/shared/redirect.F90
 endif
+elpa_test_new_interface@SUFFIX@_SOURCES = test/Fortran/test_new_interface.F90
+elpa_test_new_interface@SUFFIX@_LDADD = $(build_lib) $(FCLIBS)
+elpa_test_new_interface@SUFFIX@_FCFLAGS = $(AM_FCFLAGS) @FC_MODOUT@private_modules @FC_MODINC@private_modules
+EXTRA_elpa_test_new_interface@SUFFIX@_DEPENDENCIES = test/Fortran/elpa_print_headers.X90
+
 
 elpa1_test_real_c_version@SUFFIX@_SOURCES = test/C/elpa1_test_real_c_version.c
 elpa1_test_real_c_version@SUFFIX@_LDADD = $(build_lib) $(FCLIBS)
@@ -753,6 +762,7 @@ check_SCRIPTS = \
   elpa2_test_real_c_version@SUFFIX@.sh \
   elpa2_test_complex_c_version@SUFFIX@.sh \
   elpa_driver_real_c_version@SUFFIX@.sh \
+  elpa_test_new_interface@SUFFIX@.sh \
   elpa_driver_complex_c_version@SUFFIX@.sh
 
 if WANT_SINGLE_PRECISION_REAL

src/elpa_init.F90

+module init_elpa
+
+  private
+  public :: elpa_init, initDone
+
+  logical :: initDone = .false.
+
+  contains
+
+  subroutine elpa_init()
+
+    implicit none
+
+    ! must be done by all task using ELPA !!!
+
+    initDone = .true.
+
+  end subroutine
+
+
+end module init_elpa

src/elpa_t.F90

+module elpa_type
+ use iso_c_binding
+
+ private
+
+ public :: elpa_create, elpa_t
+
+ type :: elpa_t
+   integer(kind=c_int) :: mpi_comm_rows, mpi_comm_cols, mpi_comm_global
+   integer(kind=c_int) :: na, nev, local_nrows, local_ncols, nblk
+
+   integer(kind=c_int) :: real_kernel, complex_kernel
+
+   integer(kind=c_int) :: useQR, useGPU
+   character(6)        :: solver
+   character(8)        :: timings
+   contains
+     generic, public :: set_option => elpa_set_option_string, elpa_set_option_integer
+     procedure, private :: elpa_set_option_string, elpa_set_option_integer
+     generic, public :: get_option => elpa_get_option_string, elpa_get_option_integer
+     procedure, private :: elpa_get_option_string, elpa_get_option_integer
+
+     procedure :: get_communicators => get_communicators
+     procedure :: solve_real_double => elpa_solve_real_double
+
+ end type elpa_t
+
+ contains
+
+
+
+   function elpa_create(na, nev, local_nrows, local_ncols, nblk) result(elpa)
+     use precision
+     use init_elpa
+     use elpa2_utilities, only : DEFAULT_REAL_ELPA_KERNEL, DEFAULT_COMPLEX_ELPA_KERNEL
+     implicit none
+
+      integer(kind=ik), intent(in) :: na, nev, local_nrows, local_ncols, nblk
+      type(elpa_t) :: elpa
+
+      ! check whether init has ever been called
+      if (.not.(initDone))  then
+        print *,"ERROR: you must call elpa_init() once before creating instances of ELPA"
+        stop
+      endif
+
+      elpa%na          = na
+      elpa%nev         = nev
+      elpa%local_nrows = local_nrows
+      elpa%local_ncols = local_ncols
+      elpa%nblk        = nblk
+
+      ! some default values
+      elpa%solver         = "2stage"
+      elpa%real_kernel    = DEFAULT_REAL_ELPA_KERNEL
+      elpa%complex_kernel = DEFAULT_COMPLEX_ELPA_KERNEL
+
+      elpa%useQR          = 0
+      elpa%useGPU         = 0
+      elpa%timings        = "none"
+
+    end function
+
+    function elpa_set_option_string(self, keyword, value) result(success)
+      use iso_c_binding
+      use elpa1, only : elpa_print_times
+      implicit none
+      class(elpa_t)            :: self
+      character(*), intent(in) :: keyword
+      character(*), intent(in) :: value
+      integer(kind=c_int)      :: success
+
+      success = 0
+
+      if (trim(keyword) .eq. "solver") then
+        if (trim(value) .eq. "1stage") then
+          self%solver = "1stage"
+          success = 1
+        else if (trim(value) .eq. "2stage") then
+          self%solver = "2stage"
+          success = 1
+        else if (trim(value) .eq. "auto") then
+          self%solver = "auto "
+          success = 1
+        else
+          print *," not allowed key/value pair: ",trim(keyword),"/",trim(value)
+          success = 0
+        endif
+      else if (trim(keyword) .eq. "timings") then
+        if (trim(value) .eq. "balanced") then
+          elpa_print_times = .true.
+          success = 1
+        else if (trim(value) .eq. "detailed") then
+          print *,"detailed timings not yet implemented"
+          elpa_print_times = .false.
+          success = 1
+        else if (trim(value) .eq. "none") then
+          elpa_print_times = .false.
+          success = 1
+        else
+          print *," not allowed key/value pair: ",trim(keyword),"/",trim(value)
+          success = 0
+        endif
+      else
+        print *," not allowed key/value pair: ",trim(keyword),"/",trim(value)
+        success = 0
+      endif
+
+    end function elpa_set_option_string
+
+    function elpa_set_option_integer(self, keyword, value) result(success)
+      use iso_c_binding
+      use elpa2_utilities, only : check_allowed_real_kernels, check_allowed_complex_kernels
+      implicit none
+      class(elpa_t)                   :: self
+      character(*), intent(in)        :: keyword
+      integer(kind=c_int), intent(in) :: value
+      integer(kind=c_int)             :: success
+
+      success = 0
+
+      if (trim(keyword) .eq. "real_kernel") then
+        if (.not.(check_allowed_real_kernels(value))) then
+          self%real_kernel = value
+          success = 1
+        else
+          print *,"Setting this real_kernel is not possible"
+          success = 0
+        endif
+      else if (trim(keyword) .eq. "complex_kernel" ) then
+        if (.not.(check_allowed_complex_kernels(value))) then
+          self%complex_kernel = value
+          success = 1
+        else
+          print *,"Setting this complex_kernel is not possible"
+          success = 0
+        endif
+      else if (trim(keyword) .eq. "use_qr") then
+        if (value .eq. 1) then
+          self%useQr = 1
+          success = 1
+        else if (value .eq. 0) then
+          self%useQr = 0
+          success = 1
+        else
+          print *," not allowed key/value pair: ",trim(keyword),"/",value
+          success = 0
+        endif
+      else if (trim(keyword) .eq. "use_gpu") then
+        if (value .eq. 1) then
+          self%useGPU = 1
+          success = 1
+        else if (value .eq. 0) then
+          self%useGPU = 0
+          success = 1
+        else
+          print *," not allowed key/value pair: ",trim(keyword),"/",value
+          success = 0
+        endif
+      else
+        print *," not allowed key/value pair: ",trim(keyword),"/",value
+        success = 0
+      endif
+
+    end function elpa_set_option_integer
+
+    function elpa_get_option_string(self, keyword, value) result(success)
+      use iso_c_binding
+      use elpa1, only : elpa_print_times
+      implicit none
+      class(elpa_t)               :: self
+      character(*), intent(in)    :: keyword
+      character(*), intent(inout) :: value
+      integer(kind=c_int)         :: success
+
+      success = 0
+
+      if (trim(keyword) .eq. "solver") then
+        value = trim(self%solver)
+        success = 1
+      else if (trim(keyword) .eq. "timings") then
+        if (elpa_print_times) then
+          value = "balanced"
+          success = 1
+        else
+          ! detailed not yet implemented
+          success = 1
+        endif
+      else
+        print *," not allowed key/value pair: ",trim(keyword),"/",trim(value)
+        success = 0
+      endif
+
+    end function elpa_get_option_string
+
+    function elpa_get_option_integer(self, keyword, value) result(success)
+      use iso_c_binding
+      implicit none
+      class(elpa_t)                      :: self
+      character(*), intent(in)           :: keyword
+      integer(kind=c_int), intent(inout) :: value
+      integer(kind=c_int)                :: success
+
+      success = 0
+
+      if (trim(keyword) .eq. "real_kernel") then
+        value = self%real_kernel
+        success = 1
+      else if (trim(keyword) .eq. "complex_kernel" ) then
+        value = self%complex_kernel
+        success = 1
+      else if (trim(keyword) .eq. "use_qr") then
+        value = self%useQr
+        success = 1
+      else if (trim(keyword) .eq. "use_gpu") then
+        value =  self%useGPU
+        success = 1
+      else
+        print *," not allowed key/value pair: ",trim(keyword),"/",value
+        success = 0
+      endif
+
+    end function elpa_get_option_integer
+
+    function get_communicators(self, mpi_comm_global, my_prow, my_pcol, mpi_comm_rows, mpi_comm_cols) result(mpierr)
+      use iso_c_binding
+      use elpa_mpi
+      use elpa1, only : elpa_get_communicators
+      implicit none
+      class(elpa_t)                    :: self
+
+      integer(kind=c_int), intent(in)  :: mpi_comm_global, my_prow, my_pcol
+      integer(kind=c_int), intent(out) :: mpi_comm_rows, mpi_comm_cols
+
+      integer(kind=c_int)              :: mpierr
+
+      mpierr = elpa_get_communicators(mpi_comm_global, my_prow, my_pcol, mpi_comm_rows, mpi_comm_cols)
+
+      self%mpi_comm_rows   = mpi_comm_rows
+      self%mpi_comm_cols   = mpi_comm_cols
+      self%mpi_comm_global = mpi_comm_global
+    end function
+
+    function elpa_solve_real_double(self, a, ev, q) result(success)
+      use elpa
+
+      use iso_c_binding
+      implicit none
+      class(elpa_t)                    :: self
+
+      real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols), &
+                             ev(self%na)
+      integer(kind=c_int) :: success
+
+      logical :: successDummy
+
+      successDummy =  elpa_solve_evp_real_double(self%na, self%nev, a, self%local_nrows, ev, q,  &
+                                                 self%local_nrows,  self%nblk, self%local_ncols, &
+                                                 self%mpi_comm_rows, self%mpi_comm_cols,         &
+                                                 self%mpi_comm_global, method=trim(self%solver))
+
+      if (successDummy) then
+        success = 1
+      else
+        success = 0
+      endif
+    end function
+
+
+end module
+

test/Fortran/test_new_interface.F90

+!    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 "config-f90.h"
+!>
+!> Fortran test programm to demonstrates the use of
+!> ELPA 2 real case library.
+!> If "HAVE_REDIRECT" was defined at build time
+!> the stdout and stderr output of each MPI task
+!> can be redirected to files if the environment
+!> variable "REDIRECT_ELPA_TEST_OUTPUT" is set
+!> to "true".
+!>
+!> By calling executable [arg1] [arg2] [arg3] [arg4]
+!> one can define the size (arg1), the number of
+!> Eigenvectors to compute (arg2), and the blocking (arg3).
+!> If these values are not set default values (4000, 1500, 16)
+!> are choosen.
+!> If these values are set the 4th argument can be
+!> "output", which specifies that the EV's are written to
+!> an ascii file.
+!>
+!> The real ELPA 2 kernel is set as the default kernel.
+!> However, this can be overriden by setting
+!> the environment variable "REAL_ELPA_KERNEL" to an
+!> appropiate value.
+!>
+program test_real2_double_precision
+
+!-------------------------------------------------------------------------------
+! Standard eigenvalue problem - REAL version
+!
+! This program demonstrates the use of the ELPA module
+! together with standard scalapack routines
+!
+! Copyright of the original code rests with the authors inside the ELPA
+! consortium. The copyright of any additional modifications shall rest
+! with their original authors, but shall adhere to the licensing terms
+! distributed along with the original code in the file "COPYING".
+!
+!-------------------------------------------------------------------------------
+   use precision
+   use ELPA1
+   use ELPA2
+
+   use elpa2_utilities
+
+   use mod_check_for_gpu, only : check_for_gpu
+   use elpa_utilities, only : error_unit
+#ifdef WITH_OPENMP
+   use test_util
+#endif
+
+   use mod_read_input_parameters
+   use mod_check_correctness
+   use mod_setup_mpi
+   use mod_blacs_infrastructure
+   use mod_prepare_matrix
+   use elpa_mpi
+#ifdef HAVE_REDIRECT
+  use redirect
+#endif
+#ifdef HAVE_DETAILED_TIMINGS
+ use timings
+#endif
+ use output_types
+
+
+  use elpa_type
+  use init_elpa
+
+
+   implicit none
+
+   !-------------------------------------------------------------------------------
+   ! Please set system size parameters below!
+   ! na:   System size
+   ! nev:  Number of eigenvectors to be calculated
+   ! nblk: Blocking factor in block cyclic distribution
+   !-------------------------------------------------------------------------------
+
+   integer(kind=ik)           :: nblk
+   integer(kind=ik)           :: na, nev
+
+   integer(kind=ik)           :: np_rows, np_cols, na_rows, na_cols
+
+   integer(kind=ik)           :: myid, nprocs, my_prow, my_pcol, mpi_comm_rows, mpi_comm_cols
+   integer(kind=ik)           :: i, mpierr, my_blacs_ctxt, sc_desc(9), info, nprow, npcol
+
+   integer(kind=ik), external :: numroc
+
+   real(kind=rk8), allocatable :: a(:,:), z(:,:), tmp1(:,:), tmp2(:,:), as(:,:), ev(:)
+
+   integer(kind=ik)           :: iseed(4096) ! Random seed, size should be sufficient for every generator
+   integer(kind=ik)           :: STATUS
+#ifdef WITH_OPENMP
+   integer(kind=ik)           :: omp_get_max_threads,  required_mpi_thread_level, provided_mpi_thread_level
+#endif
+   logical                    :: successELPA, success
+   integer(kind=ik)           :: success_test
+   integer(kind=ik)           :: numberOfDevices
+   logical                    :: gpuAvailable
+   type(output_t)             :: write_to_file
+   character(len=8)           :: task_suffix
+   integer(kind=ik)           :: j
+
+
+
+  type(elpa_t) :: instance1
+  integer(kind=ik) :: use_qr
+  character(7)     :: solver
+
+
+
+#define DOUBLE_PRECISION_REAL 1
+
+   successELPA   = .true.
+   gpuAvailable  = .false.
+
+   call read_input_parameters(na, nev, nblk, write_to_file)
+
+   !-------------------------------------------------------------------------------
+   !  MPI Initialization
+   call setup_mpi(myid, nprocs)
+
+   gpuAvailable = check_for_gpu(myid, numberOfDevices)
+
+   STATUS = 0
+
+#define REALCASE
+#include "elpa_print_headers.X90"
+
+#ifdef HAVE_DETAILED_TIMINGS
+
+   ! initialise the timing functionality
+
+#ifdef HAVE_LIBPAPI
+   call timer%measure_flops(.true.)
+#endif
+
+   call timer%measure_allocated_memory(.true.)
+   call timer%measure_virtual_memory(.true.)
+   call timer%measure_max_allocated_memory(.true.)
+
+   call timer%set_print_options(&
+#ifdef HAVE_LIBPAPI
+                print_flop_count=.true., &
+                print_flop_rate=.true., &
+#endif
+                print_allocated_memory = .true. , &
+                print_virtual_memory=.true., &
+                print_max_allocated_memory=.true.)
+
+
+  call timer%enable()
+
+  call timer%start("program: test_real2_double_precision")
+#endif
+
+   !-------------------------------------------------------------------------------
+   ! Selection of number of processor rows/columns
+   ! We try to set up the grid square-like, i.e. start the search for possible
+   ! divisors of nprocs with a number next to the square root of nprocs
+   ! and decrement it until a divisor is found.
+
+   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
+      if(mod(nprocs,np_cols) == 0 ) exit
+   enddo
+   ! at the end of the above loop, nprocs is always divisible by np_cols
+
+   np_rows = nprocs/np_cols
+
+   if(myid==0) then
+      print *
+      print '(a)','Standard eigenvalue problem - REAL version'
+      print *
+      print '(3(a,i0))','Matrix size=',na,', Number of eigenvectors=',nev,', Block size=',nblk
+      print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
+      print *
+   endif
+
+   !-------------------------------------------------------------------------------
+   ! Set up BLACS context and MPI communicators
+   !
+   ! The BLACS context is only necessary for using Scalapack.
+   !
+   ! For ELPA, the MPI communicators along rows/cols are sufficient,
+   ! and the grid setup may be done in an arbitrary way as long as it is
+   ! consistent (i.e. 0<=my_prow<np_rows, 0<=my_pcol<np_cols and every
+   ! process has a unique (my_prow,my_pcol) pair).
+
+   call set_up_blacsgrid(mpi_comm_world, my_blacs_ctxt, np_rows, np_cols, &
+                         nprow, npcol, my_prow, my_pcol)
+
+   if (myid==0) then
+     print '(a)','| Past BLACS_Gridinfo.'
+   end if
+
+!   ! All ELPA routines need MPI communicators for communicating within
+!   ! rows or columns of processes, these are set in elpa_get_communicators.
+!
+!   mpierr = elpa_get_communicators(mpi_comm_world, my_prow, my_pcol, &
+!                                   mpi_comm_rows, mpi_comm_cols)
+
+   if (myid==0) then
+     print '(a)','| Past split communicator setup for rows and columns.'
+   end if
+
+   call set_up_blacs_descriptor(na ,nblk, my_prow, my_pcol, np_rows, np_cols, &
+                                na_rows, na_cols, sc_desc, my_blacs_ctxt, info)
+
+   if (myid==0) then
+     print '(a)','| Past scalapack descriptor setup.'
+   end if
+
+   !-------------------------------------------------------------------------------
+   ! Allocate matrices and set up a test matrix for the eigenvalue problem
+#ifdef HAVE_DETAILED_TIMINGS
+   call timer%start("set up matrix")
+#endif
+   allocate(a (na_rows,na_cols))
+   allocate(z (na_rows,na_cols))
+   allocate(as(na_rows,na_cols))
+
+   allocate(ev(na))
+
+   call prepare_matrix_double(na, myid, sc_desc, iseed,  a, z, as)
+
+#ifdef HAVE_DETAILED_TIMINGS
+   call timer%stop("set up matrix")
+#endif
+   ! set print flag in elpa1
+   elpa_print_times = .true.
+
+