Restore the number of OpenMP threads at the end of ELPA

This solves issue #66

Makefile.am

@@ -40,6 +40,7 @@ libelpa@SUFFIX@_private_la_SOURCES = \
   src/helpers/mod_precision.F90 \
   src/helpers/mod_mpi.F90 \
   src/helpers/mod_mpi_stubs.F90 \
+  src/helpers/mod_omp.F90 \
   src/elpa_generated_fortran_interfaces.F90 \
   src/elpa2/mod_redist_band.F90 \
   src/elpa2/mod_pack_unpack_cpu.F90 \

src/elpa1/elpa1_template.F90

@@ -66,9 +66,8 @@ function elpa_solve_evp_&
    use elpa_abstract_impl
    use elpa_mpi
    use elpa1_compute
-#ifdef WITH_OPENMP
-   use omp_lib
-#endif
+   use elpa_omp
+
    implicit none
 #include "../general/precision_kinds.F90"
    class(elpa_abstract_impl_t), intent(inout) :: obj
@@ -122,7 +121,11 @@ function elpa_solve_evp_&
    &")
 
 #ifdef WITH_OPENMP
-   !nrThreads = omp_get_max_threads()
+   ! store the number of OpenMP threads used in the calling function
+   ! restore this at the end of ELPA 2
+   omp_threads_caller = omp_get_max_threads()
+
+   ! check the number of threads that ELPA should use internally
    call obj%get("omp_threads",nrThreads,error)
    call omp_set_num_threads(nrThreads)
 #else
@@ -156,6 +159,13 @@ function elpa_solve_evp_&
      if (.not.(obj%eigenvalues_only)) then
        q(1,1) = ONE
      endif
+
+     ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+     ! store the number of OpenMP threads used in the calling function
+     ! restore this at the end of ELPA 2
+     call omp_set_num_threads(omp_threads_caller)
+#endif
      call obj%timer%stop("elpa_solve_evp_&
      &MATH_DATATYPE&
      &_1stage_&
@@ -411,6 +421,13 @@ function elpa_solve_evp_&
      endif
    endif
 
+   ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+   ! store the number of OpenMP threads used in the calling function
+   ! restore this at the end of ELPA 2
+   call omp_set_num_threads(omp_threads_caller)
+#endif
+
    call obj%timer%stop("elpa_solve_evp_&
    &MATH_DATATYPE&
    &_1stage_&

src/elpa1/elpa1_tridiag_template.F90

@@ -102,9 +102,8 @@ call prmat(na,useGpu,a_mat,a_dev,lda,matrixCols,nblk,my_prow,my_pcol,np_rows,np_
       use precision
       use elpa_abstract_impl
       use matrix_plot
-#ifdef WITH_OPENMP
-      use omp_lib
-#endif
+      use elpa_omp
+
       implicit none
 #include "../general/precision_kinds.F90"
       class(elpa_abstract_impl_t), intent(inout) :: obj

src/elpa1/elpa_cholesky_template.F90

@@ -48,9 +48,8 @@
      use elpa_mpi
      use precision
      use elpa_abstract_impl
-#ifdef WITH_OPENMP
-     use omp_lib
-#endif
+     use elpa_omp
+
      implicit none
 #include "../general/precision_kinds.F90"
       class(elpa_abstract_impl_t), intent(inout) :: obj
@@ -80,7 +79,11 @@
       &")
 
 #ifdef WITH_OPENMP
-      !nrThreads=omp_get_max_threads()
+      ! store the number of OpenMP threads used in the calling function
+      ! restore this at the end of ELPA 2
+      omp_threads_caller = omp_get_max_threads()
+
+      ! check the number of threads that ELPA should use internally
       call obj%get("omp_threads",nrThreads,error)
       call omp_set_num_threads(nrThreads)
 #else
@@ -333,6 +336,14 @@
           a(l_row1:l_rows,l_col1) = 0
         endif
       enddo
+
+      ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+      ! store the number of OpenMP threads used in the calling function
+      ! restore this at the end of ELPA 2
+      call omp_set_num_threads(omp_threads_caller)
+#endif
+      
       call obj%timer%stop("elpa_cholesky_&
       &MATH_DATATYPE&
       &_&

src/elpa1/elpa_solve_tridi_impl_public.F90

@@ -63,9 +63,7 @@
                          &_impl
       use precision
       use elpa_abstract_impl
-#ifdef WITH_OPENMP
-      use omp_lib
-#endif
+      use elpa_omp
 
       implicit none
       class(elpa_abstract_impl_t), intent(inout) :: obj
@@ -95,7 +93,12 @@
       matrixCols = obj%local_ncols
 
 #ifdef WITH_OPENMP
-      !nrThreads=omp_get_max_threads()
+      ! store the number of OpenMP threads used in the calling function
+      ! restore this at the end of ELPA 2 
+      omp_threads_caller = omp_get_max_threads()
+
+      ! check the number of threads that ELPA should use internally
+
       call obj%get("omp_threads",nrThreads,error)
 #else
       nrThreads=1
@@ -130,6 +133,15 @@
                mpi_comm_rows, mpi_comm_cols,.false., wantDebug, success, &
                nrThreads)
 
+
+      ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+      ! store the number of OpenMP threads used in the calling function
+      ! restore this at the end of ELPA 2
+      call omp_set_num_threads(omp_threads_caller)
+#endif
+
+
       call obj%timer%stop("elpa_solve_tridi_public_&
       &MATH_DATATYPE&
       &_&

src/elpa1/legacy_interface/elpa_solve_tridi.F90

@@ -62,9 +62,8 @@
       &_impl
       use elpa
       use elpa_abstract_impl
-#ifdef WITH_OPENMP
-      use omp_lib
-#endif
+      use elpa_omp
+
       implicit none
       integer(kind=ik)            :: na, nev, ldq, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols
       real(kind=REAL_DATATYPE)    :: d(na), e(na)
@@ -148,6 +147,12 @@
       endif
 
 #ifdef WITH_OPENMP
+      ! store the number of OpenMP threads used in the calling function
+      ! restore this at the end of ELPA 2
+      omp_threads_caller = omp_get_max_threads()
+      
+      ! check the number of threads that ELPA should use internally
+
       nrThreads = omp_get_max_threads()
       call obj%set("omp_threads", nrThreads, error)
 #else
@@ -167,6 +172,12 @@
       if (error /= ELPA_OK) then
         print *, "Cannot run solve_tridi"
         success = .false.
+        ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+        ! store the number of OpenMP threads used in the calling function
+        ! restore this at the end of ELPA 2
+        call omp_set_num_threads(omp_threads_caller)
+#endif
         return
       else
         success = .true.
@@ -184,6 +195,13 @@
 !        stop
       endif
 
+     ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+     ! store the number of OpenMP threads used in the calling function
+     ! restore this at the end of ELPA 2
+     call omp_set_num_threads(omp_threads_caller)
+#endif
+
      !call timer%stop("elpa_solve_tridi_&
      !&PRECISION&
      !&_legacy_interface")

src/elpa2/elpa2_template.F90

@@ -63,9 +63,8 @@
    use elpa_mpi
    use cuda_functions
    use mod_check_for_gpu
-#ifdef WITH_OPENMP
-   use omp_lib
-#endif
+   use elpa_omp
+
    use iso_c_binding
    implicit none
 #include "../general/precision_kinds.F90"
@@ -150,7 +149,11 @@
 
 
 #ifdef WITH_OPENMP
-    !nrThreads = omp_get_max_threads()
+    ! store the number of OpenMP threads used in the calling function
+    ! restore this at the end of ELPA 2
+    omp_threads_caller = omp_get_max_threads()
+
+    ! check the number of threads that ELPA should use internally
     call obj%get("omp_threads",nrThreads,error)
     call omp_set_num_threads(nrThreads)
 #else
@@ -209,6 +212,14 @@
      if (.not.(obj%eigenvalues_only)) then
        q(1,1) = ONE
      endif
+
+     ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+     ! store the number of OpenMP threads used in the calling function
+     ! restore this at the end of ELPA 2
+     call omp_set_num_threads(omp_threads_caller)
+#endif
+
      call obj%timer%stop("elpa_solve_evp_&
      &MATH_DATATYPE&
      &_2stage_&
@@ -779,6 +790,13 @@
        endif
      endif
 
+     ! restore original OpenMP settings
+#ifdef WITH_OPENMP
+    ! store the number of OpenMP threads used in the calling function
+    ! restore this at the end of ELPA 2
+    call omp_set_num_threads(omp_threads_caller)
+#endif
+
      call obj%timer%stop("elpa_solve_evp_&
      &MATH_DATATYPE&
      &_2stage_&

src/helpers/mod_omp.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.
+!
+! Author: Andreas Marek, MPCDF
+
+#include "config-f90.h"
+
+!> \brief Fortran module which exports the MPI functions to ELPA
+module elpa_omp
+#ifdef WITH_OPENMP
+   use omp_lib
+#endif
+   use iso_c_binding
+   use precision
+   implicit none
+   public
+
+   integer(kind=ik) :: omp_threads_caller
+
+
+end module
+

test/C/test_autotune.c

@@ -178,9 +178,9 @@ int main(int argc, char** argv) {
    handle = elpa_allocate();
 #else
    handle = elpa_allocate(&error);
-#endif
    assert_elpa_ok(error);
 #endif
+   assert_elpa_ok(error);
 
    /* Set parameters */
    elpa_set(handle, "na", na, &error);

test/C/test_multiple_objs.c

@@ -234,7 +234,6 @@ int main(int argc, char** argv) {
    elpa_handle_2 = elpa_allocate();
 #else
    elpa_handle_2 = elpa_allocate(&error);
-#endif
    assert_elpa_ok(error);
 #endif
 

test/Fortran/test.F90

@@ -180,7 +180,7 @@ program test
                                   do_test_hermitian_multiply
 
 #ifdef WITH_OPENMP
-   integer                    :: max_threads
+   integer                    :: max_threads, threads_caller
 #endif
 
 #ifdef SPLIT_COMM_MYSELF
@@ -532,6 +532,14 @@ program test
      do_test_cholesky = .false.
    endif
 
+
+#ifdef WITH_OPENMP
+   threads_caller = omp_get_max_threads()
+   if (myid == 0) then
+     print *,"The calling program uses ",threads_caller," threads"
+   endif
+#endif
+
    e => elpa_allocate(error)
    assert_elpa_ok(error)
 
@@ -789,6 +797,15 @@ program test
      endif
 #endif
 
+
+#ifdef WITH_OPENMP
+     if (threads_caller .ne. omp_get_max_threads()) then
+       if (myid .eq. 0) then
+         print *, " ERROR! the number of OpenMP threads has not been restored correctly"
+       endif
+       status = 1
+     endif
+#endif
      if (myid == 0) then
        print *, ""
      endif