Commit 9750a05e authored by Andreas Marek's avatar Andreas Marek
Browse files

Remove obsolete elpa_tests.F90

parent 35b8b870
......@@ -286,9 +286,7 @@ popd
%files tools
%attr(0755,root,root) %{_bindir}/elpa2_print_kernels
%attr(0755,root,root) %{_bindir}/elpa_tests
%attr(0644,root,root) %{_mandir}/man1/elpa2_print_kernels.1.gz
%attr(0644,root,root) %{_mandir}/man1/elpa_tests.1.gz
%files devel
%defattr(-,root,root)
......@@ -313,7 +311,6 @@ popd
%files -n %{name}_openmp-tools
%defattr(-,root,root)
%attr(0755,root,root) %{_bindir}/elpa2_print_kernels_openmp
%attr(0755,root,root) %{_bindir}/elpa_tests_openmp
%files -n %{name}_openmp-devel
......
.TH "elpa_tests" 1 "Thur Mar 17 2017" "ELPA" \" -*- nroff -*-
.ad l
.nh
.SH NAME
elpa_tests \- Provide all tests for the ELPA library\&.
.SH SYNOPSIS
.br
elpa_tests [--help] [datatype={real|complex}] [na=number] [nev=number] [nblk=size of block cyclic distribution] [--output_eigenvalues] [--output_eigenvectors] [--real-kernel=name_of_kernel] [--complex-kernel=name_of_kernel] [--use-gpu={0|1}] [--use-qr={0,1}] [--tests={all|solve-tridi|1stage|2stage|cholesky|invert-triangular|transpose-mulitply}]
.br
.SH "Description"
.PP
Provide in a configurable way all test for the ELPA library.
.br
It is possible to run all test which have been implemented in ELPA
If no options are set, default options are choses (and printed).
.SH "Options"
.PP
.br
.RI "--help print help information"
.br
.RI "dataype={real|complex} choose whether real or complex test cases are run. If not given, the default is real case."
.br
.RI "na=integer number choose the size of the matrix. If not given, the default na=4000 is set."
.br
.RI "nev=integer number choose number of eigenvalues/eigenvectors to compute. If not given, the default nev=1500 is set."
.br.
.RI "nblk=integer number set the size of the block cyclic distribution. If not given, the default nblk=16 is set."
.br
.RI "--output-eigenvalues if set, the computed eigenvalues will be stored in files. (default no)."
.br
.RI "--output-eigenvectors if set, the computed eigenvectors will be stored in files. (default no)."
.br
.RI "--real-kernel=string if given, use only this kernel for the real 2stage step. Available kernels can be querried with elpa2_print_kernels."
.br
.RI "--complex-kernel=string if given, use only this kernel for the complex 2stage step Available kernels can be querried with elpa2_print_kernels"
.br
.RI "--use-gpu={0|1}] switch on GPU usage. Fails, if GPUs are not available. Default no"
.br
.RI "--use-qr={0,1} use QR-decomposition in real 2stage step. Default no"
.br
.RI "--tests={all|solve-tridi|1stage|2stage| if given, run only specified tests. Default is all"
.RI " cholesky|invert-triangular|"
.RI " transpose-mulitply}"
.SH "Author"
A. Marek, MPCDF
.SH "Reporting bugs"
Report bugs to the ELPA mail elpa-library@mpcdf.mpg.de
.SH "SEE ALSO"
\fBelpa2_print_kernels\fP(1) \fBelpa_get_communicators\fP(3) \fBelpa_solve_evp_real_double\fP(3) \fBelpa_solve_evp_real_single\fP(3) \fBelpa_solve_evp_complex_double\fP(3) \fBelpa_solve_evp_complex_single\fP(3) \fBelpa_solve_evp_real_1stage_double\fP(3) \fBelpa_solve_evp_real_1stage_single\fP(3) \fBelpa_solve_evp_real_2stage_double\fP(3) \fBelpa_solve_evp_real_2stage_single\fP(3) \fBelpa_solve_evp_complex_2stage_double\fP(3) \fBelpa_solve_evp_complex_2stage_single\fP(3)
! 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"
program test_all_real
!-------------------------------------------------------------------------------
! Standard eigenvalueGproblem - 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 elpa1
use elpa2
use elpa_utilities, only : error_unit, map_global_array_index_to_local_index
use elpa2_utilities
use test_util
use test_read_input_parameters
use test_check_correctness
use test_setup_mpi
use test_blacs_infrastructure
use test_prepare_matrix
#ifdef HAVE_REDIRECT
use test_redirect
#endif
use test_output_type
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, external :: numroc
logical :: wantDebug
real(kind=rk8), allocatable :: ev(:)
real(kind=rk8), allocatable, target :: a_real(:,:), z_real(:,:), tmp1_real(:,:), tmp2_real(:,:), as_real(:,:)
complex(kind=ck8), allocatable, target :: a_complex(:,:), z_complex(:,:), tmp1_complex(:,:), tmp2_complex(:,:), as_complex(:,:)
real(kind=rk8), allocatable, target :: b_real(:,:), bs_real(:,:), c_real(:,:)
complex(kind=ck8), allocatable, target :: b_complex(:,:), bs_complex(:,:), c_complex(:,:)
complex(kind=ck8), parameter :: CZERO = (0.0_rk8,0.0_rk8), CONE = (1.0_rk8,0.0_rk8)
real(kind=rk8), allocatable, target :: d_real(:), e_real(:), ev_analytic_real(:)
real(kind=rk8), target :: diagonalELement_real, subdiagonalElement_real
complex(kind=ck8) :: diagonalElement_complex
real(kind=rk8), target :: tmp
real(kind=rk8) :: norm, normmax
real(kind=rk8), parameter :: pi = 3.141592653589793238462643383279_rk8
integer(kind=ik) :: STATUS
#ifdef WITH_OPENMP
integer(kind=ik) :: omp_get_max_threads, required_mpi_thread_level, &
provided_mpi_thread_level
#endif
type(input_options_t) :: input_options
logical :: success
character(len=8) :: task_suffix
integer(kind=ik) :: j, this_kernel
logical :: this_gpu
logical :: this_qr
integer(kind=ik) :: loctmp ,rowLocal, colLocal
real(kind=rk8) :: tStart, tEnd
real(kind=rk8) :: maxerr
logical :: gpuAvailable
#ifdef WITH_MPI
real(kind=rk8) :: pzlange, pdlange
#else
real(kind=rk8) :: zlange, dlange
#endif
!-------------------------------------------------------------------------------
success = .true.
call read_input_parameters(input_options)
na = input_options%na
nev = input_options%nev
nblk = input_options%nblk
if (input_options%justHelpMessage) then
call EXIT(0)
endif
!-------------------------------------------------------------------------------
! MPI Initialization
call setup_mpi(myid, nprocs)
STATUS = 0
!#define DATATYPE REAL !!! check here
!#include "elpa_print_headers.F90"
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 *
if (input_options%datatype .eq. 1) then
print '(a)','Standard eigenvalue problem - REAL version'
else if (input_options%datatype .eq. 2) then
print '(a)','Standard eigenvalue problem - COMPLEX version'
endif
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
if (input_options%datatype .eq. 1) then
allocate(a_real (na_rows,na_cols))
allocate(z_real (na_rows,na_cols))
allocate(as_real(na_rows,na_cols))
if (input_options%doInvertTrm .or. input_options%doTransposeMultiply) then
allocate(b_real(na_rows,na_cols))
allocate(bs_real(na_rows,na_cols))
endif
if (input_options%doTransposeMultiply) then
allocate(c_real(na_rows,na_cols))
endif
endif
if (input_options%datatype .eq. 2) then
allocate(a_complex (na_rows,na_cols))
allocate(z_complex (na_rows,na_cols))
allocate(as_complex(na_rows,na_cols))
if (input_options%doInvertTrm .or. input_options%doTransposeMultiply) then
allocate(b_complex(na_rows,na_cols))
allocate(bs_complex(na_rows,na_cols))
endif
if (input_options%doTransposeMultiply) then
allocate(c_complex(na_rows,na_cols))
endif
endif
if (input_options%datatype .eq. 1) then
allocate(d_real (na))
allocate(e_real (na))
allocate(ev_analytic_real(na))
endif
allocate(ev(na))
if (input_options%datatype .eq. 1) then
allocate(tmp1_real(na_rows,na_cols))
allocate(tmp2_real(na_rows,na_cols))
endif
if (input_options%datatype .eq. 2) then
allocate(tmp1_complex(na_rows,na_cols))
allocate(tmp2_complex(na_rows,na_cols))
endif
if (input_options%datatype .eq. 1) then
call prepare_matrix_random(na, myid, sc_desc, a_real, z_real, as_real)
if (input_options%doInvertTrm) then
b_real(:,:) = a_real(:,:)
bs_real(:,:) = a_real(:,:)
endif
endif
if (input_options%datatype .eq. 2) then
call prepare_matrix_random(na, myid, sc_desc, a_complex, z_complex, as_complex)
if (input_options%doInvertTrm) then
b_complex(:,:) = a_complex(:,:)
bs_complex(:,:) = a_complex(:,:)
endif
endif
if (input_options%doSolveTridi) then
! first the toeplitz test only in real case
! changeable numbers here would be nice
if (input_options%datatype .eq. 1) then
diagonalElement_real = 0.45_rk8
subdiagonalElement_real = 0.78_rk8
d_real(:) = diagonalElement_real
e_real(:) = subdiagonalElement_real
! set up the diagonal and subdiagonals (for general solver test)
do i=1, na ! for diagonal elements
if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a_real(rowLocal,colLocal) = diagonalElement_real
endif
enddo
do i=1, na-1
if (map_global_array_index_to_local_index(i, i+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a_real(rowLocal,colLocal) = subdiagonalElement_real
endif
enddo
do i=2, na
if (map_global_array_index_to_local_index(i, i-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a_real(rowLocal,colLocal) = subdiagonalElement_real
endif
enddo
as_real = a_real
!-------------------------------------------------------------------------------
! Calculate eigenvalues/eigenvectors
if (myid==0) then
print '(a)','| Entering elpa_solve_tridi ... '
print *
end if
#ifdef WITH_MPI
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
#endif
success = elpa_solve_tridi_double(na, nev, d_real, e_real, a_real, na_rows, nblk, na_cols, mpi_comm_rows, &
mpi_comm_cols, wantDebug)
if (.not.(success)) then
write(error_unit,*) "elpa_solve_tridi produced an error! Aborting..."
#ifdef WITH_MPI
call MPI_ABORT(mpi_comm_world, 1, mpierr)
#else
call exit(1)
#endif
endif
if (myid==0) then
print '(a)','| elpa_solve_tridi complete.'
print *
end if
ev = d_real
! analytic solution
do i=1, na
ev_analytic_real(i) = diagonalElement_real + 2.0_rk8 * subdiagonalElement_real * &
cos( pi*real(i,kind=rk8)/ real(na+1,kind=rk8) )
enddo
! sort analytic solution:
! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
! a proper sorting algorithmus might be implemented here
tmp = minval(ev_analytic_real)
loctmp = minloc(ev_analytic_real, 1)
ev_analytic_real(loctmp) = ev_analytic_real(1)
ev_analytic_real(1) = tmp
do i=2, na
tmp = ev_analytic_real(i)
do j= i, na
if (ev_analytic_real(j) .lt. tmp) then
tmp = ev_analytic_real(j)
loctmp = j
endif
enddo
ev_analytic_real(loctmp) = ev_analytic_real(i)
ev_analytic_real(i) = tmp
enddo
! compute a simple error max of eigenvalues
maxerr = 0.0_rk8
maxerr = maxval( (d_real(:) - ev_analytic_real(:))/ev_analytic_real(:) , 1)
if (maxerr .gt. 8.e-13_rk8) then
if (myid .eq. 0) then
print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
endif
status = 1
endif
endif ! datatype == 1 for toeplitz
endif ! doSolve-tridi
if (input_options%doCholesky) then
if (input_options%datatype .eq. 1) then
a_real(:,:) = 0.0_rk8
diagonalElement_real = 2.546_rk8
do i = 1, na
if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a_real(rowLocal,colLocal) = diagonalElement_real * abs(cos( pi*real(i,kind=rk8)/ real(na+1,kind=rk8) ))
endif
enddo
as_real(:,:) = a_real(:,:)
if (myid==0) then
print '(a)','| Compute real cholesky decomposition ... '
print *
end if
#ifdef WITH_MPI
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
#endif
success = elpa_cholesky_real_double(na, a_real, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, .true.)
if (.not.(success)) then
write(error_unit,*) "elpa_cholseky_real produced an error! Aborting..."
#ifdef WITH_MPI
call MPI_ABORT(mpi_comm_world, 1, mpierr)
#else
call exit(1)
#endif
endif
if (myid==0) then
print '(a)','| Real cholesky decomposition complete.'
print *
end if
tmp1_real(:,:) = 0.0_rk8
#ifdef WITH_MPI
call pdtran(na, na, 1.0_rk8, a_real, 1, 1, sc_desc, 0.0_rk8, tmp1_real, 1, 1, sc_desc)
#else
tmp1_real = transpose(a_real)
#endif
! tmp2 = a * a**T
#ifdef WITH_MPI
call pdgemm("N","N", na, na, na, 1.0_rk8, a_real, 1, 1, sc_desc, tmp1_real, 1, 1, &
sc_desc, 0.0_rk8, tmp2_real, 1, 1, sc_desc)
#else
call dgemm("N","N", na, na, na, 1.0_rk8, a_real, na, tmp1_real, na, 0.0_rk8, tmp2_real, na)
#endif
! compare tmp2 with original matrix
tmp2_real(:,:) = tmp2_real(:,:) - as_real(:,:)
#ifdef WITH_MPI
norm = pdlange("M",na, na, tmp2_real, 1, 1, sc_desc, tmp1_real)
#else
norm = dlange("M", na, na, tmp2_real, na_rows, tmp1_real)
#endif
#ifdef WITH_MPI
call mpi_allreduce(norm,normmax,1,MPI_REAL8,MPI_MAX,MPI_COMM_WORLD,mpierr)
#else
normmax = norm
#endif
if (myid .eq. 0) then
print *," Maximum error of result: ", normmax
endif
if (normmax .gt. 5e-12_rk8) then
status = 1
endif
endif ! real case
if (input_options%datatype .eq. 2) then
a_complex(:,:) = CONE - CONE
diagonalElement_complex = (2.546_rk8, 0.0_rk8)
do i = 1, na
if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a_complex(rowLocal,colLocal) = diagonalElement_complex * abs(cos( pi*real(i,kind=rk8)/ real(na+1,kind=rk8) ))
endif
enddo
as_complex(:,:) = a_complex(:,:)
if (myid==0) then
print '(a)','| Compute complex cholesky decomposition ... '
print *
end if
#ifdef WITH_MPI
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
#endif
success = elpa_cholesky_complex_double(na, a_complex, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, .true.)
if (.not.(success)) then
write(error_unit,*) " elpa_cholesky_complex produced an error! Aborting..."
#ifdef WITH_MPI
call MPI_ABORT(mpi_comm_world, 1, mpierr)
#else
call exit(1)
#endif
endif
if (myid==0) then
print '(a)','| Solve cholesky decomposition complete.'
print *
end if
tmp1_complex(:,:) = 0.0_ck8
! tmp1 = a**H
#ifdef WITH_MPI
call pztranc(na, na, CONE, a_complex, 1, 1, sc_desc, CZERO, tmp1_complex, 1, 1, sc_desc)
#else
tmp1_complex = transpose(conjg(a_complex))
#endif
! tmp2 = a * a**H
#ifdef WITH_MPI
call pzgemm("N","N", na, na, na, CONE, a_complex, 1, 1, sc_desc, tmp1_complex, 1, 1, &
sc_desc, CZERO, tmp2_complex, 1, 1, sc_desc)
#else
call zgemm("N","N", na, na, na, CONE, a_complex, na, tmp1_complex, na, CZERO, tmp2_complex, na)
#endif
! compare tmp2 with c
tmp2_complex(:,:) = tmp2_complex(:,:) - as_complex(:,:)
#ifdef WITH_MPI
norm = pzlange("M",na, na, tmp2_complex, 1, 1, sc_desc,tmp1_complex)
#else
norm = zlange("M",na, na, tmp2_complex, na_rows, tmp1_complex)
#endif
#ifdef WITH_MPI
call mpi_allreduce(norm,normmax,1,MPI_REAL8,MPI_MAX,MPI_COMM_WORLD,mpierr)
#else
normmax = norm
#endif
if (myid .eq. 0) then
print *," Maximum error of result: ", normmax
endif