Move resort_ev to a module

Makefile.am

@@ -59,6 +59,7 @@ libelpa@SUFFIX@_private_la_SOURCES = \
   src/general/mod_elpa_skewsymmetric_blas.F90 \
   src/solve_tridi/mod_global_product.F90 \
   src/solve_tridi/mod_global_gather.F90 \
+  src/solve_tridi/mod_resort_ev.F90 \
   src/elpa_index.c
 
 libelpa@SUFFIX@_private_la_SOURCES += src/elpa_c_interface.c
@@ -681,6 +682,7 @@ EXTRA_DIST = \
   src/elpa1/elpa1_compute_template.F90 \
   src/solve_tridi/global_product_template.F90 \
   src/solve_tridi/global_gather_template.F90 \
+  src/solve_tridi/resort_ev_template.F90 \
   src/elpa1/elpa1_merge_systems_real_template.F90 \
   src/elpa1/elpa1_solve_tridi_real_template.F90 \
   src/elpa1/elpa1_template.F90 \

src/elpa1/elpa1_merge_systems_real_template.F90

@@ -65,6 +65,7 @@
       use elpa_blas_interfaces
       use global_product
       use global_gather
+      use resort_ev
 #ifdef WITH_OPENMP
       use omp_lib
 #endif
@@ -266,7 +267,8 @@
         ! Rearrange eigenvectors
         call resort_ev_&
         &PRECISION &
-                       (obj, idx, na)
+                       (obj, idx, na, na, p_col_out, q, ldq, matrixCols, l_rows, l_rqe, &
+                        l_rqs, mpi_comm_cols, p_col, l_col, l_col_out)
 
         call obj%timer%stop("merge_systems" // PRECISION_SUFFIX)
 
@@ -437,7 +439,9 @@
         enddo
         call resort_ev_&
         &PRECISION&
-        &(obj, idxq1, na)
+        &(obj, idxq1, na, na, p_col_out, q, ldq, matrixCols, l_rows, l_rqe, &
+          l_rqs, mpi_comm_cols, p_col, l_col, l_col_out)
+
       else if (na1>2) then
 
         ! Solve secular equation
@@ -928,6 +932,7 @@
         end subroutine add_tmp_&
         &PRECISION
 
+#if 0
         subroutine resort_ev_&
         &PRECISION&
         &(obj, idx_ev, nLength)
@@ -1006,6 +1011,7 @@
           check_deallocate("resort_ev: qtmp",istat, errorMessage)
         end subroutine resort_ev_&
         &PRECISION
+#endif
 
         subroutine transform_columns_&
         &PRECISION&

src/solve_tridi/mod_resort_ev.F90

+
+#include "config-f90.h"
+module resort_ev
+  use precision
+  implicit none
+  private
+
+  public :: resort_ev_double
+#if defined(WANT_SINGLE_PRECISION_REAL) || defined(WANT_SINGLE_PRECISION_COMPLEX)
+  public :: resort_ev_single
+#endif
+
+  contains
+
+! real double precision first
+#define DOUBLE_PRECISION_REAL
+#define REALCASE
+#define DOUBLE_PRECISION
+#include "../general/precision_macros.h"
+#include "./resort_ev_template.F90"
+#undef DOUBLE_PRECISION_REAL
+#undef REALCASE
+#undef DOUBLE_PRECISION
+
+#ifdef WANT_SINGLE_PRECISION_REAL
+! real single precision first
+#define SINGLE_PRECISION_REAL
+#define REALCASE
+#define SINGLE_PRECISION
+#include "../general/precision_macros.h"
+#include "./resort_ev_template.F90"
+#undef SINGLE_PRECISION_REAL
+#undef REALCASE
+#undef SINGLE_PRECISION
+#endif
+
+end module

src/solve_tridi/resort_ev_template.F90

+#include "../general/error_checking.inc"
+
+subroutine resort_ev_&
+&PRECISION&
+  &(obj, idx_ev, nLength, na, p_col_out, q, ldq, matrixCols, l_rows, l_rqe, l_rqs, &
+    mpi_comm_cols, p_col, l_col, l_col_out)
+    use precision
+#ifdef WITH_OPENMP
+    use elpa_omp
+#endif
+    use elpa_mpi
+    use ELPA_utilities
+    use elpa_abstract_impl
+    implicit none
+    class(elpa_abstract_impl_t), intent(inout) :: obj
+    integer(kind=ik), intent(in)               :: nLength, na
+    integer(kind=ik), intent(in)               :: ldq, matrixCols, l_rows, l_rqe, l_rqs
+    integer(kind=ik), intent(in)               :: mpi_comm_cols
+    integer(kind=ik), intent(in)               :: p_col(na), l_col(na), l_col_out(na)
+#ifdef WITH_MPI
+    integer(kind=MPI_KIND)                     :: mpierrMPI, my_pcolMPI
+    integer(kind=ik)                           :: mpierr
+#endif
+    integer(kind=ik)                           :: my_pcol
+#ifdef USE_ASSUMED_SIZE
+    real(kind=REAL_DATATYPE), intent(inout)    :: q(ldq,*)
+#else
+    real(kind=REAL_DATATYPE), intent(inout)    :: q(ldq,matrixCols)
+#endif
+    integer(kind=ik), intent(in)               :: p_col_out(na)
+    integer(kind=ik)                           :: idx_ev(nLength)
+    integer(kind=ik)                           :: i, nc, pc1, pc2, lc1, lc2, l_cols_out
+
+    real(kind=REAL_DATATYPE), allocatable      :: qtmp(:,:)
+    integer(kind=ik)                           :: istat
+    character(200)                             :: errorMessage
+
+    if (l_rows==0) return ! My processor column has no work to do
+
+    ! Resorts eigenvectors so that q_new(:,i) = q_old(:,idx_ev(i))
+
+    l_cols_out = COUNT(p_col_out(1:na)==my_pcol)
+    allocate(qtmp(l_rows,l_cols_out), stat=istat, errmsg=errorMessage)
+    check_allocate("resort_ev: qtmp",istat, errorMessage)
+
+    nc = 0
+
+    do i=1,na
+
+      pc1 = p_col(idx_ev(i))
+      lc1 = l_col(idx_ev(i))
+      pc2 = p_col_out(i)
+
+      if (pc2<0) cycle ! This column is not needed in output
+
+      if (pc2==my_pcol) nc = nc+1 ! Counter for output columns
+
+      if (pc1==my_pcol) then
+        if (pc2==my_pcol) then
+          ! send and recieve column are local
+          qtmp(1:l_rows,nc) = q(l_rqs:l_rqe,lc1)
+        else
+#ifdef WITH_MPI
+          call obj%timer%start("mpi_communication")
+          call mpi_send(q(l_rqs,lc1), int(l_rows,kind=MPI_KIND), MPI_REAL_PRECISION, pc2, int(mod(i,4096),kind=MPI_KIND), &
+                              int(mpi_comm_cols,kind=MPI_KIND), mpierr)
+          call obj%timer%stop("mpi_communication")
+#else /* WITH_MPI */
+
+#endif /* WITH_MPI */
+        endif
+      else if (pc2==my_pcol) then
+#ifdef WITH_MPI
+        call obj%timer%start("mpi_communication")
+        call mpi_recv(qtmp(1,nc), int(l_rows,kind=MPI_KIND), MPI_REAL_PRECISION, pc1, int(mod(i,4096),kind=MPI_KIND), &
+                      int(mpi_comm_cols,kind=MPI_KIND), MPI_STATUS_IGNORE, mpierr)
+        call obj%timer%stop("mpi_communication")
+#else /* WITH_MPI */
+        qtmp(1:l_rows,nc) = q(l_rqs:l_rqe,lc1)
+#endif /* WITH_MPI */
+       endif
+     enddo
+
+     ! Insert qtmp into (output) q
+
+     nc = 0
+
+     do i=1,na
+
+       pc2 = p_col_out(i)
+       lc2 = l_col_out(i)
+
+       if (pc2==my_pcol) then
+         nc = nc+1
+         q(l_rqs:l_rqe,lc2) = qtmp(1:l_rows,nc)
+       endif
+     enddo
+
+     deallocate(qtmp, stat=istat, errmsg=errorMessage)
+     check_deallocate("resort_ev: qtmp",istat, errorMessage)
+   end subroutine resort_ev_&
+        &PRECISION
+