test.F90 13.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
!    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"

! Define one of TEST_REAL or TEST_COMPLEX
! Define one of TEST_SINGLE or TEST_DOUBLE
! Define one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE
! Define TEST_GPU \in [0, 1]
49
! Define either TEST_ALL_KERNELS or a TEST_KERNEL \in [any valid kernel]
50

51
52
#if !(defined(TEST_REAL) ^ defined(TEST_COMPLEX))
error: define exactly one of TEST_REAL or TEST_COMPLEX
53
54
55
56
57
58
#endif

#if !(defined(TEST_SINGLE) ^ defined(TEST_DOUBLE))
error: define exactly one of TEST_SINGLE or TEST_DOUBLE
#endif

Pavel Kus's avatar
Pavel Kus committed
59
60
#if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE) ^ defined(TEST_SCALAPACK_ALL))
error: define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE or TEST_SCALAPACK_ALL 
61
62
#endif

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#ifdef TEST_SOLVER_1STAGE
#ifdef TEST_ALL_KERNELS
error: TEST_ALL_KERNELS cannot be defined for TEST_SOLVER_1STAGE
#endif
#ifdef TEST_KERNEL
error: TEST_KERNEL cannot be defined for TEST_SOLVER_1STAGE
#endif
#endif

#ifdef TEST_SOLVER_2STAGE
#if !(defined(TEST_KERNEL) ^ defined(TEST_ALL_KERNELS))
error: define either TEST_ALL_KERNELS or a valid TEST_KERNEL
#endif
#endif


79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#ifdef TEST_SINGLE
#  define EV_TYPE real(kind=C_FLOAT)
#  ifdef TEST_REAL
#    define MATRIX_TYPE real(kind=C_FLOAT)
#  else
#    define MATRIX_TYPE complex(kind=C_FLOAT_COMPLEX)
#  endif
#else
#  define EV_TYPE real(kind=C_DOUBLE)
#  ifdef TEST_REAL
#    define MATRIX_TYPE real(kind=C_DOUBLE)
#  else
#    define MATRIX_TYPE complex(kind=C_DOUBLE_COMPLEX)
#  endif
#endif

95
96
97
98
99
100
101
#ifdef TEST_REAL
#define KERNEL_KEY "real_kernel"
#endif
#ifdef TEST_COMPLEX
#define KERNEL_KEY "complex_kernel"
#endif

102
103
104
105
#include "assert.h"

program test
   use elpa
106
107
108
109
110
111
112

   use test_util
   use test_setup_mpi
   use test_prepare_matrix
   use test_read_input_parameters
   use test_blacs_infrastructure
   use test_check_correctness
Pavel Kus's avatar
Pavel Kus committed
113
   use test_analytic
114
#ifdef WITH_SCALAPACK_TESTS
Pavel Kus's avatar
Pavel Kus committed
115
   use test_scalapack
116
#endif
117

118
119
120
#ifdef HAVE_REDIRECT
   use test_redirect
#endif
121
122
123
124
125
126
127
128
129
130
131
132
133
   implicit none

   ! matrix dimensions
   integer :: na, nev, nblk

   ! mpi
   integer :: myid, nprocs
   integer :: na_cols, na_rows  ! local matrix size
   integer :: np_cols, np_rows  ! number of MPI processes per column/row
   integer :: my_prow, my_pcol  ! local MPI task position (my_prow, my_pcol) in the grid (0..np_cols -1, 0..np_rows -1)
   integer :: mpierr

   ! blacs
134
   integer :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
135
136
137
138
139
140
141

   ! The Matrix
   MATRIX_TYPE, allocatable :: a(:,:), as(:,:)
   ! eigenvectors
   MATRIX_TYPE, allocatable :: z(:,:)
   ! eigenvalues
   EV_TYPE, allocatable :: ev(:)
142
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_CHOLESKY)
143
   EV_TYPE, allocatable :: d(:), sd(:), ev_analytic(:), ds(:), sds(:)
Andreas Marek's avatar
Andreas Marek committed
144
   EV_TYPE              :: diagonalELement, subdiagonalElement
145
146
#endif

147

148
   integer :: error, status
149
150
151

   type(output_t) :: write_to_file
   class(elpa_t), pointer :: e
152
#ifdef TEST_ALL_KERNELS
153
   integer :: i
154
155
156
157
#endif
#ifdef TEST_ALL_LAYOUTS
   character(len=1), parameter :: layouts(2) = [ 'C', 'R' ]
   integer :: i_layout
158
#endif
159
   integer :: kernel
160
   character(len=1) :: layout
161

162
   call read_input_parameters_traditional(na, nev, nblk, write_to_file)
163
   call setup_mpi(myid, nprocs)
164
165
166
167
#ifdef HAVE_REDIRECT
#ifdef WITH_MPI
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
     call redirect_stdout(myid)
168
#endif
169
#endif
170

171
172
173
174
175
   if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
     print *, "ELPA API version not supported"
     stop 1
   endif

176
177
178
179
180
   if (myid == 0) then
     print '((a,i0))', 'Program ' // TEST_CASE
     print *, ""
   endif

181
#ifdef TEST_ALL_LAYOUTS
182
   do i_layout = 1, size(layouts)               ! layouts
183
     layout = layouts(i_layout)
184
     do np_cols = 1, nprocs                     ! factors
185
186
187
188
189
       if (mod(nprocs,np_cols) /= 0 ) then
         cycle
       endif
#else
   layout = 'C'
190
191
192
   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
      if(mod(nprocs,np_cols) == 0 ) exit
   enddo
193
#endif
194
195

   np_rows = nprocs/np_cols
196
   assert(nprocs == np_rows * np_cols)
197

198
199
200
201
   if (myid == 0) then
     print '((a,i0))', 'Matrix size: ', na
     print '((a,i0))', 'Num eigenvectors: ', nev
     print '((a,i0))', 'Blocksize: ', nblk
202
#ifdef WITH_MPI
203
204
     print '((a,i0))', 'Num MPI proc: ', nprocs
     print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
205
     print '(a)',      'Process layout: ' // layout
206
#endif
207
208
209
     print *,''
   endif

210
211
   call set_up_blacsgrid(mpi_comm_world, np_rows, np_cols, layout, &
                         my_blacs_ctxt, my_prow, my_pcol)
212
213
214
215

   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)

Pavel Kus's avatar
Pavel Kus committed
216
217
   allocate(a (na_rows,na_cols))
   allocate(as(na_rows,na_cols))
218
219
220
   allocate(z (na_rows,na_cols))
   allocate(ev(na))

221
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_CHOLESKY)
222
223
224
225
226
   allocate(d (na), ds(na))
   allocate(sd (na), sds(na))
   allocate(ev_analytic(na))
#endif

227
228
229
230
   a(:,:) = 0.0
   z(:,:) = 0.0
   ev(:) = 0.0

231
#ifdef TEST_EIGENVECTORS
Pavel Kus's avatar
Pavel Kus committed
232
233
#ifdef TEST_MATRIX_ANALYTIC
   call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
234
235
236
237
238
239
240
241
242
     as(:,:) = a
#else
   if (nev .ge. 1) then
     call prepare_matrix(na, myid, sc_desc, a, z, as)
   else
     ! zero eigenvectors and not analytic test => toeplitz matrix
#ifdef TEST_SINGLE
     diagonalElement = 0.45_c_float
     subdiagonalElement =  0.78_c_float
Pavel Kus's avatar
Pavel Kus committed
243
#else
244
245
246
247
248
249
250
     diagonalElement = 0.45_c_double
     subdiagonalElement =  0.78_c_double
#endif
     call prepare_toeplitz_matrix(na, diagonalElement, subdiagonalElement, &
                                  d, sd, ds, sds, a, as, nblk, np_rows, &
                                  np_cols, my_prow, my_pcol)
   endif
251
#endif
Pavel Kus's avatar
Pavel Kus committed
252
#endif
253

254
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL)
Andreas Marek's avatar
Andreas Marek committed
255
256
257
258
259
260
261
262
263
264
265

#ifdef TEST_SINGLE
   diagonalElement = 0.45_c_float
   subdiagonalElement =  0.78_c_float
#else
   diagonalElement = 0.45_c_double
   subdiagonalElement =  0.78_c_double
#endif
   call prepare_toeplitz_matrix(na, diagonalElement, subdiagonalElement, &
                                d, sd, ds, sds, a, as, nblk, np_rows, &
                                np_cols, my_prow, my_pcol)
266
#endif /* EIGENVALUES OR TRIDIAGONAL */
267

268
269
270
271
272
273
274
275
276
277
278
279
280
281
#if defined(TEST_CHOLESKY)

#ifdef TEST_SINGLE
   diagonalElement = 2.546_c_float
   subdiagonalElement =  0.0_c_float
#else
   diagonalElement = 2.546_c_double
   subdiagonalElement =  0.0_c_double
#endif
   call prepare_toeplitz_matrix(na, diagonalElement, subdiagonalElement, &
                                d, sd, ds, sds, a, as, nblk, np_rows, &
                                np_cols, my_prow, my_pcol)
#endif /* TEST_CHOLESKY */

282
283
   e => elpa_allocate()

284
285
286
287
288
289
290
291
292
293
   call e%set("na", na, error)
   assert_elpa_ok(error)
   call e%set("nev", nev, error)
   assert_elpa_ok(error)
   call e%set("local_nrows", na_rows, error)
   assert_elpa_ok(error)
   call e%set("local_ncols", na_cols, error)
   assert_elpa_ok(error)
   call e%set("nblk", nblk, error)
   assert_elpa_ok(error)
294
295

#ifdef WITH_MPI
296
297
298
299
300
301
   call e%set("mpi_comm_parent", MPI_COMM_WORLD, error)
   assert_elpa_ok(error)
   call e%set("process_row", my_prow, error)
   assert_elpa_ok(error)
   call e%set("process_col", my_pcol, error)
   assert_elpa_ok(error)
302
#endif
303

Andreas Marek's avatar
Andreas Marek committed
304
305
   call e%set("timings",1)

306
307
308
309
310
311
312
   assert_elpa_ok(e%setup())

#ifdef TEST_SOLVER_1STAGE
   call e%set("solver", ELPA_SOLVER_1STAGE)
#else
   call e%set("solver", ELPA_SOLVER_2STAGE)
#endif
313
   assert_elpa_ok(error)
314

315
316
   call e%set("gpu", TEST_GPU, error)
   assert_elpa_ok(error)
317

318
319
   if (myid == 0) print *, ""

320
#ifdef TEST_ALL_KERNELS
321
   do i = 0, elpa_option_cardinality(KERNEL_KEY)  ! kernels
322
     kernel = elpa_option_enumerate(KERNEL_KEY, i)
323
#endif
324
#ifdef TEST_KERNEL
325
     kernel = TEST_KERNEL
326
#endif
327
328

#ifdef TEST_SOLVER_2STAGE
329
     call e%set(KERNEL_KEY, kernel, error)
330
331
332
#ifdef TEST_KERNEL
     assert_elpa_ok(error)
#else
333
334
335
     if (error /= ELPA_OK) then
       cycle
     endif
336
337
338
339
340
341
     ! actually used kernel might be different if forced via environment variables
     call e%get(KERNEL_KEY, kernel)
#endif
     if (myid == 0) then
       print *, elpa_int_value_to_string(KERNEL_KEY, kernel) // " kernel"
     endif
342
343
#endif

344
#ifdef TEST_ALL_KERNELS
345
     call e%timer_start(elpa_int_value_to_string(KERNEL_KEY, kernel))
346
#endif
347

348
     ! The actual solve step
349
350
#ifdef TEST_EIGENVECTORS
     call e%timer_start("e%eigenvectors()")
Pavel Kus's avatar
Pavel Kus committed
351
352
353
#ifdef TEST_SCALAPACK_ALL
     call solve_scalapack_all(na, a, sc_desc, ev, z)
#else
354
     call e%eigenvectors(a, ev, z, error)
Pavel Kus's avatar
Pavel Kus committed
355
#endif
356
     call e%timer_stop("e%eigenvectors()")
357
#endif /* TEST_EIGENVECTORS */
358
359
360

#ifdef TEST_EIGENVALUES
     call e%timer_start("e%eigenvalues()")
361
     call e%eigenvalues(a, ev, error)
362
     call e%timer_stop("e%eigenvalues()")
363
#endif
364
365
366

#if defined(TEST_SOLVE_TRIDIAGONAL)
     call e%timer_start("e%solve_tridiagonal()")
367
     call e%solve_tridiagonal(d, sd, z, error)
368
     call e%timer_stop("e%solve_tridiagonal()")
369
370
371
     ev(:) = d(:)
#endif

372
373
374
375
376
377
#if defined(TEST_CHOLESKY)
     call e%timer_start("e%cholesky()")
     call e%cholesky(a, error)
     call e%timer_stop("e%cholesky()")
#endif

Pavel Kus's avatar
Pavel Kus committed
378

379
380
     assert_elpa_ok(error)

381
#ifdef TEST_ALL_KERNELS
382
383
384
     call e%timer_stop(elpa_int_value_to_string(KERNEL_KEY, kernel))
#endif

385
     if (myid .eq. 0) then
386
#ifdef TEST_ALL_KERNELS
387
       call e%print_times(elpa_int_value_to_string(KERNEL_KEY, kernel))
388
389
#else /* TEST_ALL_KERNELS */

390
#ifdef TEST_EIGENVECTORS
391
       call e%print_times("e%eigenvectors()")
392
#endif
393
#ifdef TEST_EIGENVALUES
394
395
       call e%print_times("e%eigenvalues()")
#endif
396
397
#ifdef TEST_SOLVE_TRIDIAGONAL
       call e%print_times("e%solve_tridiagonal()")
398
#endif
399
400
#ifdef TEST_CHOLESKY
       call e%print_times("e%cholesky()")
401
#endif
402
#endif /* TEST_ALL_KERNELS */
403
     endif
404

405
#ifdef TEST_EIGENVECTORS
Pavel Kus's avatar
Pavel Kus committed
406
407
408
#ifdef TEST_MATRIX_ANALYTIC
     status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
#else
409
410
411
412
413
414
415
     if (nev .ge. 1) then
       status = check_correctness(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
     else
       ! zero eigenvectors and no analytic test => toeplitz
       status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &
         subdiagonalElement, ev, z, myid)
     endif
Andreas Marek's avatar
Andreas Marek committed
416
#endif
417
     call check_status(status, myid)
418
419
#endif

420
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL)
Andreas Marek's avatar
Andreas Marek committed
421
422
     status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &
         subdiagonalElement, ev, z, myid)
423
     call check_status(status, myid)
424

425
#ifdef TEST_SOLVE_TRIDIAGONAL
426
     ! check eigenvectors
Pavel Kus's avatar
Pavel Kus committed
427
     status = check_correctness(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
428
     call check_status(status, myid)
429
430
#endif
#endif
431

432
433
434
435
436
#if defined(TEST_CHOLESKY)
     status = check_correctness_cholesky(na, a, as, na_rows, sc_desc, myid )
     call check_status(status, myid)
#endif

437
438
439
440
     if (myid == 0) then
       print *, ""
     endif

441
442
#ifdef TEST_ALL_KERNELS
     a(:,:) = as(:,:)
443
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_CHOLESKY)
444
445
446
     d = ds
     sd = sds
#endif
447
   end do ! kernels
448
#endif
Andreas Marek's avatar
Andreas Marek committed
449

450
451
452
453
454
455
456
   call elpa_deallocate(e)

   deallocate(a)
   deallocate(as)
   deallocate(z)
   deallocate(ev)

457
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVECTORS) || defined(TEST_CHOLESKY)
458
459
460
461
462
   deallocate(d, ds)
   deallocate(sd, sds)
   deallocate(ev_analytic)
#endif

463
464
465
466
467
468
469
#ifdef TEST_ALL_LAYOUTS
   end do ! factors
   end do ! layouts
#endif

   call elpa_uninit()

470
471
472
473
474
475
476
#ifdef WITH_MPI
   call blacs_gridexit(my_blacs_ctxt)
   call mpi_finalize(mpierr)
#endif

   call exit(status)

477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
   contains

     subroutine check_status(status, myid)
       implicit none
       integer, intent(in) :: status, myid
       integer :: mpierr
       if (status /= 0) then
         if (myid == 0) print *, "Result incorrect!"
#ifdef WITH_MPI
         call mpi_finalize(mpierr)
#endif
         call exit(status)
       endif
     end subroutine

492
end program