Commit 4d35ddc8 authored by Andreas Marek's avatar Andreas Marek

More cleanup of test programs

parent 11cae8ed
/* 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"
#include <stdio.h>
#include <stdlib.h>
#ifdef WITH_MPI
#include <mpi.h>
#endif
#include <math.h>
#include <elpa/elpa_legacy.h>
#include <test/shared/generated.h>
#include <complex.h>
int main(int argc, char** argv) {
int myid;
int nprocs;
#ifndef WITH_MPI
int MPI_COMM_WORLD;
#endif
int na, nev, nblk;
int status;
int np_cols, np_rows, np_colsStart;
int my_blacs_ctxt, nprow, npcol, my_prow, my_pcol;
int mpierr;
int my_mpi_comm_world;
int mpi_comm_rows, mpi_comm_cols;
int info, *sc_desc;
int na_rows, na_cols;
float startVal;
complex *a, *z, *as;
float *ev;
int useGPU;
int success;
int i;
int THIS_COMPLEX_ELPA_KERNEL_API;
#ifdef WITH_MPI
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
#else
nprocs = 1;
myid =0;
MPI_COMM_WORLD=1;
#endif
na = 1000;
nev = 500;
nblk = 16;
if (myid == 0) {
printf("This is the c version of an ELPA test-programm\n");
printf("\n");
printf("It will call the ELPA complex solver for a matrix\n");
printf("of matrix size %d. It will compute %d eigenvalues\n",na,nev);
printf("and uses a blocksize of %d\n",nblk);
printf("\n");
printf("This is an example program with much less functionality\n");
printf("as it's Fortran counterpart. It's only purpose is to show how \n");
printf("to evoke ELPA1 from a c programm\n");
printf("\n");
}
status = 0;
startVal = sqrt((float) nprocs);
np_colsStart = (int) round(startVal);
for (np_cols=np_colsStart;np_cols>1;np_cols--){
if (nprocs %np_cols ==0){
break;
}
}
np_rows = nprocs/np_cols;
if (myid == 0) {
printf("\n");
printf("Number of processor rows %d, cols %d, total %d \n",np_rows,np_cols,nprocs);
}
/* set up blacs */
/* convert communicators before */
#ifdef WITH_MPI
my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
#else
my_mpi_comm_world = 1;
#endif
set_up_blacsgrid_f(my_mpi_comm_world, &my_blacs_ctxt, &np_rows, &np_cols, &nprow, &npcol, &my_prow, &my_pcol);
if (myid == 0) {
printf("\n");
printf("Past BLACS_Gridinfo...\n");
printf("\n");
}
/* get the ELPA row and col communicators. */
/* These are NOT usable in C without calling the MPI_Comm_f2c function on them !! */
#ifdef WITH_MPI
my_mpi_comm_world = MPI_Comm_c2f(MPI_COMM_WORLD);
#endif
mpierr = elpa_get_communicators(my_mpi_comm_world, my_prow, my_pcol, &mpi_comm_rows, &mpi_comm_cols);
if (myid == 0) {
printf("\n");
printf("Past split communicator setup for rows and columns...\n");
printf("\n");
}
sc_desc = malloc(9*sizeof(int));
set_up_blacs_descriptor_f(na, nblk, my_prow, my_pcol, np_rows, np_cols, &na_rows, &na_cols, sc_desc, my_blacs_ctxt, &info);
if (myid == 0) {
printf("\n");
printf("Past scalapack descriptor setup...\n");
printf("\n");
}
/* allocate the matrices needed for elpa */
if (myid == 0) {
printf("\n");
printf("Allocating matrices with na_rows=%d and na_cols=%d\n",na_rows, na_cols);
printf("\n");
}
a = malloc(na_rows*na_cols*sizeof(complex));
z = malloc(na_rows*na_cols*sizeof(complex));
as = malloc(na_rows*na_cols*sizeof(complex));
ev = malloc(na*sizeof(float));
prepare_matrix_complex_single_f(na, myid, na_rows, na_cols, sc_desc, a, z, as);
if (myid == 0) {
printf("\n");
printf("Entering ELPA 1stage complex solver\n");
printf("\n");
}
#ifdef WITH_MPI
mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
useGPU = 0;
THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
success = elpa_solve_evp_complex_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_COMPLEX_ELPA_KERNEL_API, useGPU,"1stage");
if (success != 1) {
printf("error in ELPA solve \n");
#ifdef WITH_MPI
mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
#endif
}
if (myid == 0) {
printf("\n");
printf("1stage ELPA complex solver complete\n");
printf("\n");
}
for (i=0;i<na_rows*na_cols;i++){
a[i] = as[i];
z[i] = as[i];
}
if (myid == 0) {
printf("\n");
printf("Entering ELPA 2stage complex solver\n");
printf("\n");
}
#ifdef WITH_MPI
mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
useGPU = 0;
THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
success = elpa_solve_evp_complex_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_COMPLEX_ELPA_KERNEL_API, useGPU, "2stage");
if (success != 1) {
printf("error in ELPA solve \n");
#ifdef WITH_MPI
mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
#endif
}
if (myid == 0) {
printf("\n");
printf("2stage ELPA complex solver complete\n");
printf("\n");
}
for (i=0;i<na_rows*na_cols;i++){
a[i] = as[i];
z[i] = as[i];
}
if (myid == 0) {
printf("\n");
printf("Entering auto-chosen ELPA complex solver\n");
printf("\n");
}
#ifdef WITH_MPI
mpierr = MPI_Barrier(MPI_COMM_WORLD);
#endif
useGPU = 0;
THIS_COMPLEX_ELPA_KERNEL_API = ELPA_2STAGE_COMPLEX_DEFAULT;
success = elpa_solve_evp_complex_single(na, nev, a, na_rows, ev, z, na_rows, nblk, na_cols, mpi_comm_rows, mpi_comm_cols, my_mpi_comm_world, THIS_COMPLEX_ELPA_KERNEL_API, useGPU, "auto");
if (success != 1) {
printf("error in ELPA solve \n");
#ifdef WITH_MPI
mpierr = MPI_Abort(MPI_COMM_WORLD, 99);
#endif
}
if (myid == 0) {
printf("\n");
printf("Auto-chosen ELPA complex solver complete\n");
printf("\n");
}
/* check the results */
status = check_correctness_complex_single_f(na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid);
if (status !=0){
printf("The computed EVs are not correct !\n");
}
if (status ==0){
if (myid == 0) {
printf("All ok!\n");
}
}
free(sc_desc);
free(a);
free(z);
free(as);
free(ev);
#ifdef WITH_MPI
MPI_Finalize();
#endif
return 0;
}
......@@ -154,8 +154,6 @@ program test_complex2_double_banded
STATUS = 0
#include "../elpa_print_headers.X90"
#ifdef HAVE_DETAILED_TIMINGS
! initialise the timing functionality
......
......@@ -82,7 +82,6 @@ program test_real2_double_banded
use precision
use elpa
use mod_check_for_gpu, only : check_for_gpu
#ifdef WITH_OPENMP
use test_util
#endif
......@@ -127,7 +126,6 @@ program test_real2_double_banded
#endif
integer(kind=ik) :: success
integer(kind=ik) :: numberOfDevices
logical :: gpuAvailable
type(output_t) :: write_to_file
character(len=8) :: task_suffix
integer(kind=ik) :: j
......@@ -136,7 +134,6 @@ program test_real2_double_banded
class(elpa_t), pointer :: e
#define DOUBLE_PRECISION_REAL 1
gpuAvailable = .false.
call read_input_parameters_traditional(na, nev, nblk, write_to_file)
......@@ -144,12 +141,9 @@ program test_real2_double_banded
! 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
......
......@@ -154,8 +154,6 @@ program test_complex2_single_banded
STATUS = 0
#include "../elpa_print_headers.X90"
#ifdef HAVE_DETAILED_TIMINGS
! initialise the timing functionality
......
......@@ -149,7 +149,6 @@ program test_real2_single_banded
STATUS = 0
#define REALCASE
#include "../elpa_print_headers.X90"
#ifdef HAVE_DETAILED_TIMINGS
......
! 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 1 real case library.
!> This program can read a matrix from an ascii
!> file and computes then the Eigenvectors.
!> 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".
!>
program read_real
!-------------------------------------------------------------------------------
! Standard eigenvalue problem - REAL version
!
! This program demonstrates the use of the ELPA module
! together with standard scalapack routines
!-------------------------------------------------------------------------------
use precision
use elpa1
use elpa_utilities, only : error_unit
#ifdef WITH_OPENMP
use test_util
#endif
#ifdef HAVE_REDIRECT
use redirect
#endif
#ifdef HAVE_DETAILED_TIMINGS
use timings
#endif
#ifdef HAVE_MPI_MODULE
use mpi
implicit none
#else
implicit none
include 'mpif.h'
#endif
!-------------------------------------------------------------------------------
! Please set system size parameters below!
! nblk: Blocking factor in block cyclic distribution
!-------------------------------------------------------------------------------
integer(kind=ik), parameter :: nblk = 16
!-------------------------------------------------------------------------------
! Local Variables
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, lenarg
integer(kind=ik), external :: numroc
real(kind=rk) :: err, errmax
real(kind=rk), allocatable :: a(:,:), z(:,:), tmp1(:,:), tmp2(:,:), as(:,:), ev(:)
character*256 :: filename
#ifdef WITH_OPENMP
integer(kind=iK) :: omp_get_max_threads, required_mpi_thread_level, provided_mpi_thread_level
#endif
!-------------------------------------------------------------------------------
! MPI Initialization
#ifndef WITH_OPENMP
call mpi_init(mpierr)
#else
required_mpi_thread_level = MPI_THREAD_MULTIPLE
call mpi_init_thread(required_mpi_thread_level, &
provided_mpi_thread_level, mpierr)
if (required_mpi_thread_level .ne. provided_mpi_thread_level) then
write(error_unit,*) "MPI ERROR: MPI_THREAD_MULTIPLE is not provided on this system"
write(error_unit,*) " only ", mpi_thread_level_name(provided_mpi_thread_level), " is available"
call exit(77)
endif
#endif
call mpi_comm_rank(mpi_comm_world,myid,mpierr)
call mpi_comm_size(mpi_comm_world,nprocs,mpierr)
#ifdef HAVE_REDIRECT
if (check_redirect_environment_variable()) then
if (myid .eq. 0) then
print *," "
print *,"Redirection of mpi processes is used"
print *," "
if (create_directories() .ne. 1) then
write(error_unit,*) "Unable to create directory for stdout and stderr!"
stop 1
endif
endif
call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
call redirect_stdout(myid)
endif
#endif
#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")
#endif
!-------------------------------------------------------------------------------
! Get the name of the input file containing the matrix and open input file
! Please note:
! get_command_argument is a FORTRAN 2003 intrinsic which may not be implemented
! for every Fortran compiler!!!
if(myid==0) then
call get_command_argument(1,filename,lenarg,info)
if(info/=0) then
write(error_unit,*) 'Usage: test_real matrix_file'
call mpi_abort(mpi_comm_world,1,mpierr)
endif
open(10,file=filename,action='READ',status='OLD',iostat=info)
if(info/=0) then
write(error_unit,*) 'Error: Unable to open ',trim(filename)
call mpi_abort(mpi_comm_world,1,mpierr)
endif
endif
call mpi_barrier(mpi_comm_world, mpierr) ! Just for safety
!-------------------------------------------------------------------------------
! 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))','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).
my_blacs_ctxt = mpi_comm_world
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in elpa_get_communicators
call elpa_get_communicators(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols)
! Read matrix size
if(myid==0) read(10,*) na
call mpi_bcast(na, 1, mpi_integer, 0, mpi_comm_world, mpierr)
! Quick check for plausibility
if(na<=0 .or. na>10000000) then
if(myid==0) write(error_unit,*) 'Illegal value for matrix size: ',na
call mpi_finalize(mpierr)
stop 1
endif
if(myid==0) print *,'Matrix size: ',na
! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that.
na_rows = numroc(na, nblk, my_prow, 0, np_rows)
na_cols = numroc(na, nblk, my_pcol, 0, np_cols)
! Set up a scalapack descriptor for the checks below.
! For ELPA the following restrictions hold: