Test "eigenvalues" and "solve_tridiagonal"

The routines "eigenvalues" and "solve_tridiagonal" are now
also tested directly with the new API (and not only via the
legacy API -> new API mapping)

Todo: adapt test generator script to contain some logic

generate_automake_test_programs.py

@@ -19,7 +19,14 @@ gpu_flag = {
         1 : "-DTEST_GPU=1",
 }
 
-for g, p, d, s in product(sorted(gpu_flag.keys()),
+test_type_flag = {
+        "eigenvectors" : "-D__EIGENVECTORS",
+        "eigenvalues"  : "-D__EIGENVALUES",
+        "solve_tridiagonal"  : "-D__SOLVE_TRIDIAGONAL",
+}
+
+for g, t, p, d, s in product(sorted(gpu_flag.keys()),
+                             sorted(test_type_flag.keys()),
                              sorted(prec_flag.keys()),
                              sorted(domain_flag.keys()),
                              sorted(solver_flag.keys())):
@@ -45,7 +52,7 @@ for g, p, d, s in product(sorted(gpu_flag.keys()),
                 raise Exception("Oh no!")
             endifs += 1
 
-        name = "test_{0}_{1}_{2}{3}{4}".format(d, p, s, "" if kernel == "nokernel" else "_" + kernel, "_gpu" if g else "")
+        name = "test_{0}_{1}_{2}_{3}{4}{5}".format(d, p, t, s, "" if kernel == "nokernel" else "_" + kernel, "_gpu" if g else "")
         print("noinst_PROGRAMS += " + name)
         print("check_SCRIPTS += " + name + ".sh")
         print(name + "_SOURCES = test/Fortran/test.F90")
@@ -54,6 +61,7 @@ for g, p, d, s in product(sorted(gpu_flag.keys()),
         print("  " + " \\\n  ".join([
             domain_flag[d],
             prec_flag[p],
+            test_type_flag[t],
             solver_flag[s],
             gpu_flag[g]] + extra_flags))
 

test/Fortran/test.F90

@@ -132,6 +132,17 @@ program test
    MATRIX_TYPE, allocatable :: z(:,:)
    ! eigenvalues
    EV_TYPE, allocatable :: ev(:)
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+   EV_TYPE, allocatable :: d(:), sd(:), ev_analytic(:), ds(:), sds(:)
+   EV_TYPE              :: diagonalELement, subdiagonalElement, tmp, maxerr
+#ifdef TEST_DOUBLE
+   EV_TYPE, parameter   :: pi = 3.141592653589793238462643383279_rk8
+#else
+   EV_TYPE, parameter   :: pi = 3.1415926535897932_rk4
+#endif
+   integer              :: loctmp ,rowLocal, colLocal, j,ii
+#endif
+
 
    integer :: error, status
 
@@ -142,6 +153,9 @@ program test
 #endif
    integer :: kernel
 
+#if defined(TEST_COMPLEX) && defined(__SOLVE_TRIDIAGONAL)
+   stop 77
+#endif
    call read_input_parameters_traditional(na, nev, nblk, write_to_file)
    call setup_mpi(myid, nprocs)
 
@@ -170,11 +184,54 @@ program test
    allocate(z (na_rows,na_cols))
    allocate(ev(na))
 
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+   allocate(d (na), ds(na))
+   allocate(sd (na), sds(na))
+   allocate(ev_analytic(na))
+#endif
+
    a(:,:) = 0.0
    z(:,:) = 0.0
    ev(:) = 0.0
 
+#ifdef __EIGENVECTORS
    call prepare_matrix(na, myid, sc_desc, a, z, as)
+#endif
+
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+   ! set up simple toeplitz matrix
+#ifdef TEST_DOUBLE
+   diagonalElement = 0.45_rk8
+   subdiagonalElement =  0.78_rk8
+#else
+   diagonalElement = 0.45_rk4
+   subdiagonalElement =  0.78_rk4
+#endif
+
+   d(:) = diagonalElement
+   sd(:) = subdiagonalElement
+
+   ! set up the diagonal and subdiagonals (for general solver test)
+   do ii=1, na ! for diagonal elements
+     if (map_global_array_index_to_local_index(ii, ii, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
+       a(rowLocal,colLocal) = diagonalElement
+     endif
+   enddo
+   do ii=1, na-1
+     if (map_global_array_index_to_local_index(ii, ii+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
+       a(rowLocal,colLocal) = subdiagonalElement
+     endif
+   enddo
+
+   do ii=2, na
+     if (map_global_array_index_to_local_index(ii, ii-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
+       a(rowLocal,colLocal) = subdiagonalElement
+     endif
+   enddo
+   ds = d
+   sds = sd
+   as = a
+#endif
 
    if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
      print *, "ELPA API version not supported"
@@ -238,35 +295,140 @@ program test
 
      call e%timer_start(elpa_int_value_to_string(KERNEL_KEY, kernel))
 #else
+
+#ifdef __EIGENVECTORS
      call e%timer_start("e%eigenvectors()")
+#endif
+#ifdef __EIGENVALUES
+     call e%timer_start("e%eigenvalues()")
+#endif
+#ifdef __SOLVE_TRIDIAGONAL
+     call e%timer_start("e%solve_tridiagonal()")
+#endif
 #endif
 
      ! The actual solve step
+#ifdef __EIGENVECTORS
      call e%eigenvectors(a, ev, z, error)
+#endif
+#ifdef __EIGENVALUES
+     call e%eigenvalues(a, ev, error)
+#endif
+#if defined(__SOLVE_TRIDIAGONAL) && !defined(TEST_COMPLEX)
+     call e%solve_tridiagonal(d, sd, z, error)
+     ev(:) = d(:)
+#endif
+
      assert_elpa_ok(error)
 
 #ifdef TEST_SOLVER_2STAGE
      call e%timer_stop(elpa_int_value_to_string(KERNEL_KEY, kernel))
 #else
+#ifdef __EIGENVECTORS
      call e%timer_stop("e%eigenvectors()")
+#endif
+#ifdef __EIGENVALUES
+     call e%timer_stop("e%eigenvalues()")
+#endif
+#ifdef __SOLVE_TRIDIAGONAL
+     call e%timer_stop("e%solve_tridiagonal()")
+#endif
 #endif
 
      if (myid .eq. 0) then
 #ifdef TEST_SOLVER_2STAGE
        call e%print_times(elpa_int_value_to_string(KERNEL_KEY, kernel))
 #else
+#ifdef __EIGENVECTORS
        call e%print_times("e%eigenvectors()")
+#endif
+#ifdef __EIGENVALUES
+       call e%print_times("e%eigenvalues()")
+#endif
+#ifdef __SOLVE_TRIDIAGONAL
+     call e%print_times("e%solve_tridiagonal()")
+#endif
 #endif
      endif
+
+#ifdef __EIGENVECTORS
      status = check_correctness(na, nev, as, z, ev, sc_desc, myid)
      if (status /= 0) then
        if (myid == 0) print *, "Result incorrect!"
        call exit(status)
      endif
      if (myid == 0) print *, ""
+#endif
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+     status = 0
+     ! analytic solution
+     do ii=1, na
+#ifdef TEST_DOUBLE
+       ev_analytic(ii) = diagonalElement + 2.0 * subdiagonalElement *cos( pi*real(ii,kind=rk8)/ real(na+1,kind=rk8) )
+#else
+       ev_analytic(ii) = diagonalElement + 2.0 * subdiagonalElement *cos( pi*real(ii,kind=rk4)/ real(na+1,kind=rk4) )
+#endif
+     enddo
+
+     ! sort analytic solution:
+
+     ! this hack is neihter elegant, nor optimized: for huge matrixes it might be expensive
+     ! a proper sorting algorithmus might be implemented here
+
+     tmp    = minval(ev_analytic)
+     loctmp = minloc(ev_analytic, 1)
+
+     ev_analytic(loctmp) = ev_analytic(1)
+     ev_analytic(1) = tmp
+
+     do ii=2, na
+       tmp = ev_analytic(ii)
+       do j= ii, na
+         if (ev_analytic(j) .lt. tmp) then
+           tmp    = ev_analytic(j)
+           loctmp = j
+         endif
+       enddo
+       ev_analytic(loctmp) = ev_analytic(ii)
+       ev_analytic(ii) = tmp
+     enddo
+
+     ! compute a simple error max of eigenvalues
+     maxerr = 0.0
+     maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)
+
+#ifdef TEST_DOUBLE
+     if (maxerr .gt. 8.e-13) then
+#else
+     if (maxerr .gt. 8.e-4) then
+#endif
+       status = 1
+       if (myid .eq. 0) then
+         print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
+       endif
+     endif
+
+     if (status /= 0) then
+       call exit(status)
+     endif
+#ifdef __SOLVE_TRIDIAGONAL
+     ! check eigenvectors
+     status = check_correctness(na, nev, as, z, ev, sc_desc, myid)
+     if (status /= 0) then
+       if (myid == 0) print *, "Result incorrect!"
+       call exit(status)
+     endif
+     if (myid == 0) print *, ""
+#endif
+
+#endif
 
 #ifdef TEST_ALL_KERNELS
      a(:,:) = as(:,:)
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+     d = ds
+     sd = sds
+#endif
    end do
 #endif
 
@@ -278,6 +440,12 @@ program test
    deallocate(z)
    deallocate(ev)
 
+#ifdef __EIGENVALUES
+   deallocate(d, ds)
+   deallocate(sd, sds)
+   deallocate(ev_analytic)
+#endif
+
 #ifdef WITH_MPI
    call blacs_gridexit(my_blacs_ctxt)
    call mpi_finalize(mpierr)
@@ -285,4 +453,60 @@ program test
 
    call exit(status)
 
+#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
+   contains
+
+     !Processor col for global col number
+     pure function pcol(global_col, nblk, np_cols) result(local_col)
+       implicit none
+       integer(kind=c_int), intent(in) :: global_col, nblk, np_cols
+       integer(kind=c_int)             :: local_col
+       local_col = MOD((global_col-1)/nblk,np_cols)
+     end function
+
+     !Processor row for global row number
+     pure function prow(global_row, nblk, np_rows) result(local_row)
+       implicit none
+       integer(kind=c_int), intent(in) :: global_row, nblk, np_rows
+       integer(kind=c_int)             :: local_row
+       local_row = MOD((global_row-1)/nblk,np_rows)
+     end function
+
+     function map_global_array_index_to_local_index(iGLobal, jGlobal, iLocal, jLocal , nblk, np_rows, np_cols, my_prow, my_pcol) &
+       result(possible)
+       implicit none
+
+       integer(kind=c_int)              :: pi, pj, li, lj, xi, xj
+       integer(kind=c_int), intent(in)  :: iGlobal, jGlobal, nblk, np_rows, np_cols, my_prow, my_pcol
+       integer(kind=c_int), intent(out) :: iLocal, jLocal
+       logical                       :: possible
+
+       possible = .true.
+       iLocal = 0
+       jLocal = 0
+
+       pi = prow(iGlobal, nblk, np_rows)
+
+       if (my_prow .ne. pi) then
+         possible = .false.
+         return
+       endif
+
+       pj = pcol(jGlobal, nblk, np_cols)
+
+       if (my_pcol .ne. pj) then
+         possible = .false.
+         return
+       endif
+       li = (iGlobal-1)/(np_rows*nblk) ! block number for rows
+       lj = (jGlobal-1)/(np_cols*nblk) ! block number for columns
+
+       xi = mod( (iGlobal-1),nblk)+1   ! offset in block li
+       xj = mod( (jGlobal-1),nblk)+1   ! offset in block lj
+
+       iLocal = li * nblk + xi
+       jLocal = lj * nblk + xj
+
+     end function
+#endif
 end program