test.F90 26.2 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

59
60
#if !(defined(TEST_SOLVER_1STAGE) ^ defined(TEST_SOLVER_2STAGE) ^ defined(TEST_SCALAPACK_ALL) ^ defined(TEST_SCALAPACK_PART))
error: define exactly one of TEST_SOLVER_1STAGE or TEST_SOLVER_2STAGE or TEST_SCALAPACK_ALL or TEST_SCALAPACK_PART
61
62
#endif

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#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

78
79
80
#ifdef TEST_GENERALIZED_DECOMP_EIGENPROBLEM
#define TEST_GENERALIZED_EIGENPROBLEM
#endif
81

82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#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

98
99
100
101
102
103
104
#ifdef TEST_REAL
#define KERNEL_KEY "real_kernel"
#endif
#ifdef TEST_COMPLEX
#define KERNEL_KEY "complex_kernel"
#endif

Andreas Marek's avatar
WIP    
Andreas Marek committed
105
106
107

#define TEST_INTEGER integer(kind=C_int64_t)

108
109
110
#include "assert.h"

program test
Andreas Marek's avatar
WIP    
Andreas Marek committed
111
   !use elpa
112

Andreas Marek's avatar
WIP    
Andreas Marek committed
113
114
115
   !use test_util
   !use test_setup_mpi
   !use test_prepare_matrix
116
   use test_read_input_parameters
Andreas Marek's avatar
WIP    
Andreas Marek committed
117
118
119
   !use test_blacs_infrastructure
   !use test_check_correctness
   !use test_analytic
120
#ifdef WITH_SCALAPACK_TESTS
Andreas Marek's avatar
WIP    
Andreas Marek committed
121
   !use test_scalapack
122
#endif
123

124
#ifdef HAVE_REDIRECT
Andreas Marek's avatar
WIP    
Andreas Marek committed
125
   !use test_redirect
126
127
#endif
#ifdef WITH_OPENMP
Andreas Marek's avatar
WIP    
Andreas Marek committed
128
   !use omp_lib
129
#endif
130
   use precision_for_tests
131

132
133
134
   implicit none

   ! matrix dimensions
Andreas Marek's avatar
WIP    
Andreas Marek committed
135
   TEST_INTEGER      :: na, nev, nblk
136
137

   ! mpi
Andreas Marek's avatar
WIP    
Andreas Marek committed
138
139
140
141
142
   TEST_INTEGER      :: myid, nprocs
   TEST_INTEGER      :: na_cols, na_rows  ! local matrix size
   TEST_INTEGER      :: np_cols, np_rows  ! number of MPI processes per column/row
   TEST_INTEGER      :: my_prow, my_pcol  ! local MPI task position (my_prow, my_pcol) in the grid (0..np_cols -1, 0..np_rows -1)
   TEST_INTEGER      :: mpierr
143
144

   ! blacs
Andreas Marek's avatar
WIP    
Andreas Marek committed
145
   TEST_INTEGER     :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
146
147

   ! The Matrix
148
   MATRIX_TYPE, allocatable    :: a(:,:), as(:,:)
149
#if defined(TEST_HERMITIAN_MULTIPLY)
150
   MATRIX_TYPE, allocatable    :: b(:,:), c(:,:)
151
152
153
#endif
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
   MATRIX_TYPE, allocatable    :: b(:,:), bs(:,:)
154
#endif
155
   ! eigenvectors
156
   MATRIX_TYPE, allocatable    :: z(:,:)
157
   ! eigenvalues
158
   EV_TYPE, allocatable        :: ev(:)
159

160
   logical                     :: check_all_evals, skip_check_correctness
161

162
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
163
164
   EV_TYPE, allocatable        :: d(:), sd(:), ds(:), sds(:)
   EV_TYPE                     :: diagonalELement, subdiagonalElement
165
#endif
166

Andreas Marek's avatar
WIP    
Andreas Marek committed
167
   TEST_INTEGER      :: error, status
168

169
170
   type(output_t)              :: write_to_file
   class(elpa_t), pointer      :: e
171
#ifdef TEST_ALL_KERNELS
Andreas Marek's avatar
WIP    
Andreas Marek committed
172
   TEST_INTEGER      :: i
173
174
175
#endif
#ifdef TEST_ALL_LAYOUTS
   character(len=1), parameter :: layouts(2) = [ 'C', 'R' ]
Andreas Marek's avatar
WIP    
Andreas Marek committed
176
   TEST_INTEGER      :: i_layout
177
#endif
Andreas Marek's avatar
WIP    
Andreas Marek committed
178
   TEST_INTEGER      :: kernel
179
   character(len=1)            :: layout
180
181
   logical                     :: do_test_numeric_residual, do_test_numeric_residual_generalized, &
                                  do_test_analytic_eigenvalues, &
Andreas Marek's avatar
Andreas Marek committed
182
183
184
185
186
                                  do_test_analytic_eigenvalues_eigenvectors,   &
                                  do_test_frank_eigenvalues,  &
                                  do_test_toeplitz_eigenvalues, do_test_cholesky,   &
                                  do_test_hermitian_multiply

187
#ifdef WITH_OPENMP
Andreas Marek's avatar
WIP    
Andreas Marek committed
188
   TEST_INTEGER      :: max_threads, threads_caller
189
190
#endif

191
#ifdef SPLIT_COMM_MYSELF
Andreas Marek's avatar
WIP    
Andreas Marek committed
192
   TEST_INTEGER      :: mpi_comm_rows, mpi_comm_cols, mpi_string_length, mpierr2
193
194
195
   character(len=MPI_MAX_ERROR_STRING) :: mpierr_string
#endif

196
   call read_input_parameters_traditional(na, nev, nblk, write_to_file, skip_check_correctness)
Andreas Marek's avatar
WIP    
Andreas Marek committed
197
   !call setup_mpi(myid, nprocs)
198
199
200
201
#ifdef HAVE_REDIRECT
#ifdef WITH_MPI
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
     call redirect_stdout(myid)
202
#endif
203
#endif
204

205
   check_all_evals = .true.
206

Andreas Marek's avatar
Andreas Marek committed
207
208

   do_test_numeric_residual = .false.
209
   do_test_numeric_residual_generalized = .false.
Andreas Marek's avatar
Andreas Marek committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
   do_test_analytic_eigenvalues = .false.
   do_test_analytic_eigenvalues_eigenvectors = .false.
   do_test_frank_eigenvalues = .false.
   do_test_toeplitz_eigenvalues = .false. 

   do_test_cholesky = .false.
#if defined(TEST_CHOLESKY)
   do_test_cholesky = .true.
#endif
   do_test_hermitian_multiply = .false.
#if defined(TEST_HERMITIAN_MULTIPLY)
   do_test_hermitian_multiply = .true.
#endif

   status = 0
225
226
227
228
229
   if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
     print *, "ELPA API version not supported"
     stop 1
   endif

230
231
232
233
234
   if (myid == 0) then
     print '((a,i0))', 'Program ' // TEST_CASE
     print *, ""
   endif

235
#ifdef TEST_ALL_LAYOUTS
236
   do i_layout = 1, size(layouts)               ! layouts
237
     layout = layouts(i_layout)
238
     do np_cols = 1, nprocs                     ! factors
239
240
241
242
243
       if (mod(nprocs,np_cols) /= 0 ) then
         cycle
       endif
#else
   layout = 'C'
244
245
246
   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
      if(mod(nprocs,np_cols) == 0 ) exit
   enddo
247
#endif
248
249

   np_rows = nprocs/np_cols
250
   assert(nprocs == np_rows * np_cols)
251

252
253
254
255
   if (myid == 0) then
     print '((a,i0))', 'Matrix size: ', na
     print '((a,i0))', 'Num eigenvectors: ', nev
     print '((a,i0))', 'Blocksize: ', nblk
256
#ifdef WITH_MPI
257
258
     print '((a,i0))', 'Num MPI proc: ', nprocs
     print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
259
     print '(a)',      'Process layout: ' // layout
260
#endif
261
262
263
     print *,''
   endif

Andreas Marek's avatar
Andreas Marek committed
264
#if TEST_QR_DECOMPOSITION == 1
265
266
267
268
269
270

#if TEST_GPU == 1
#ifdef WITH_MPI
     call mpi_finalize(mpierr)
#endif
     stop 77
Andreas Marek's avatar
Andreas Marek committed
271
#endif /* TEST_GPU */
272
273
274
275
276
277
278
279
280
281
282
283
284
   if (nblk .lt. 64) then
     if (myid .eq. 0) then
       print *,"At the moment QR decomposition need blocksize of at least 64"
     endif
     if ((na .lt. 64) .and. (myid .eq. 0)) then
       print *,"This is why the matrix size must also be at least 64 or only 1 MPI task can be used"
     endif

#ifdef WITH_MPI
     call mpi_finalize(mpierr)
#endif
     stop 77
   endif
Andreas Marek's avatar
Andreas Marek committed
285
286
#endif /* TEST_QR_DECOMPOSITION */

287

Andreas Marek's avatar
WIP    
Andreas Marek committed
288
289
   !call set_up_blacsgrid(mpi_comm_world, np_rows, np_cols, layout, &
   !                      my_blacs_ctxt, my_prow, my_pcol)
290

Andreas Marek's avatar
WIP    
Andreas Marek committed
291
292
   !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)
293

Pavel Kus's avatar
Pavel Kus committed
294
295
   allocate(a (na_rows,na_cols))
   allocate(as(na_rows,na_cols))
296
297
298
   allocate(z (na_rows,na_cols))
   allocate(ev(na))

299
300
301
302
303
#ifdef TEST_HERMITIAN_MULTIPLY
   allocate(b (na_rows,na_cols))
   allocate(c (na_rows,na_cols))
#endif

304
305
306
307
308
#ifdef TEST_GENERALIZED_EIGENPROBLEM
   allocate(b (na_rows,na_cols))
   allocate(bs (na_rows,na_cols))
#endif

309
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
310
311
312
313
   allocate(d (na), ds(na))
   allocate(sd (na), sds(na))
#endif

314
315
316
317
   a(:,:) = 0.0
   z(:,:) = 0.0
   ev(:) = 0.0

Andreas Marek's avatar
Andreas Marek committed
318
319
320
321
322
#if defined(TEST_MATRIX_RANDOM) && !defined(TEST_SOLVE_TRIDIAGONAL) && !defined(TEST_CHOLESKY) && !defined(TEST_EIGENVALUES)
   ! the random matrix can be used in allmost all tests; but for some no
   ! correctness checks have been implemented; do not allow these
   ! combinations
   ! RANDOM + TEST_SOLVE_TRIDIAGONAL: we need a TOEPLITZ MATRIX
323
   ! RANDOM + TEST_CHOLESKY: wee need SPD matrix
Andreas Marek's avatar
Andreas Marek committed
324
325
326
327
   ! RANDOM + TEST_EIGENVALUES: no correctness check known

   ! We also have to take care of special case in TEST_EIGENVECTORS
#if !defined(TEST_EIGENVECTORS)
Andreas Marek's avatar
WIP    
Andreas Marek committed
328
    !call prepare_matrix_random(na, myid, sc_desc, a, z, as)
Andreas Marek's avatar
Andreas Marek committed
329
330
#else /* TEST_EIGENVECTORS */
    if (nev .ge. 1) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
331
     !call prepare_matrix_random(na, myid, sc_desc, a, z, as)
Andreas Marek's avatar
Andreas Marek committed
332
333
334
#ifndef TEST_HERMITIAN_MULTIPLY
    do_test_numeric_residual = .true.
#endif
335
   else
Andreas Marek's avatar
Andreas Marek committed
336
337
338
339
340
     if (myid .eq. 0) then
       print *,"At the moment with the random matrix you need nev >=1"
     endif
#ifdef WITH_MPI
     call mpi_finalize(mpierr)
341
#endif
Andreas Marek's avatar
Andreas Marek committed
342
343
344
     stop 77
   endif
#endif /* TEST_EIGENVECTORS */
345
346
347
348
    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
    do_test_frank_eigenvalues = .false.
    do_test_toeplitz_eigenvalues = .false.
Andreas Marek's avatar
Andreas Marek committed
349
#endif /* (TEST_MATRIX_RANDOM) */
350
351

#if defined(TEST_MATRIX_RANDOM) && defined(TEST_CHOLESKY)
Andreas Marek's avatar
WIP    
Andreas Marek committed
352
353
     !call prepare_matrix_random_spd(na, myid, sc_desc, a, z, as, &
     !            nblk, np_rows, np_cols, my_prow, my_pcol)
354
355
356
357
358
359
    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
    do_test_frank_eigenvalues = .false.
    do_test_toeplitz_eigenvalues = .false.
#endif /* TEST_MATRIX_RANDOM and TEST_CHOLESKY */

360
#if defined(TEST_MATRIX_RANDOM) && defined(TEST_GENERALIZED_EIGENPROBLEM)
Andreas Marek's avatar
WIP    
Andreas Marek committed
361
362
363
   !! call prepare_matrix_random(na, myid, sc_desc, a, z, as)
   ! call prepare_matrix_random_spd(na, myid, sc_desc, b, z, bs, &
   !              nblk, np_rows, np_cols, my_prow, my_pcol)
364
365
366
367
368
369
370
371
    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
    do_test_frank_eigenvalues = .false.
    do_test_toeplitz_eigenvalues = .false.
    do_test_numeric_residual = .false.
    do_test_numeric_residual_generalized = .true.
#endif /* TEST_MATRIX_RANDOM and TEST_GENERALIZED_EIGENPROBLEM */

372
#if defined(TEST_MATRIX_RANDOM) && (defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVALUES))
Andreas Marek's avatar
Andreas Marek committed
373
#error "Random matrix is not allowed in this configuration"
374
#endif
375

Andreas Marek's avatar
Andreas Marek committed
376
377
378
379
380
381
#if defined(TEST_MATRIX_ANALYTIC)  && !defined(TEST_SOLVE_TRIDIAGONAL) && !defined(TEST_CHOLESKY)
   ! the analytic matrix can be used in allmost all tests; but for some no
   ! correctness checks have been implemented; do not allow these
   ! combinations
   ! ANALYTIC + TEST_SOLVE_TRIDIAGONAL: we need a TOEPLITZ MATRIX
   ! ANALTIC  + TEST_CHOLESKY: no correctness check yet implemented
382

Andreas Marek's avatar
WIP    
Andreas Marek committed
383
   !call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
384
   as(:,:) = a
385

Andreas Marek's avatar
Andreas Marek committed
386
387
388
389
   do_test_numeric_residual = .false.
   do_test_analytic_eigenvalues_eigenvectors = .false.
#ifndef TEST_HERMITIAN_MULTIPLY
   do_test_analytic_eigenvalues = .true.
390
#endif
Andreas Marek's avatar
Andreas Marek committed
391
392
393
#if defined(TEST_EIGENVECTORS)
   if (nev .ge. 1) then
     do_test_analytic_eigenvalues_eigenvectors = .true.
Pavel Kus's avatar
Pavel Kus committed
394
     do_test_analytic_eigenvalues = .false.
Andreas Marek's avatar
Andreas Marek committed
395
396
397
   else
     do_test_analytic_eigenvalues_eigenvectors = .false.
   endif
398
#endif
Andreas Marek's avatar
Andreas Marek committed
399
400
   do_test_frank_eigenvalues = .false.
   do_test_toeplitz_eigenvalues = .false.
401
#endif /* TEST_MATRIX_ANALYTIC */
Andreas Marek's avatar
Andreas Marek committed
402
403
404
#if defined(TEST_MATRIX_ANALYTIC) && (defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_CHOLESKY))
#error "Analytic matrix is not allowd in this configuration"
#endif
Andreas Marek's avatar
Andreas Marek committed
405

Andreas Marek's avatar
Andreas Marek committed
406
407
#if defined(TEST_MATRIX_TOEPLITZ)
   ! The Toeplitz matrix works in each test
Andreas Marek's avatar
Andreas Marek committed
408
409
410
411
412
413
414
#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
415

416
! actually we test cholesky for diagonal matrix only
417
418
#if defined(TEST_CHOLESKY)
#ifdef TEST_SINGLE
419
420
  diagonalElement = (2.546_c_float, 0.0_c_float)
  subdiagonalElement =  (0.0_c_float, 0.0_c_float)
421
#else
422
423
  diagonalElement = (2.546_c_double, 0.0_c_double)
  subdiagonalElement =  (0.0_c_double, 0.0_c_double)
424
#endif
Andreas Marek's avatar
Andreas Marek committed
425
426
#endif /* TEST_CHOLESKY */

Andreas Marek's avatar
WIP    
Andreas Marek committed
427
428
429
   !call prepare_matrix_toeplitz(na, diagonalElement, subdiagonalElement, &
   !                             d, sd, ds, sds, a, as, nblk, np_rows, &
   !                             np_cols, my_prow, my_pcol)
430
431


Andreas Marek's avatar
Andreas Marek committed
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
   do_test_numeric_residual = .false.
#if defined(TEST_EIGENVECTORS)
   if (nev .ge. 1) then
     do_test_numeric_residual = .true.
   else
     do_test_numeric_residual = .false.
   endif
#endif

   do_test_analytic_eigenvalues = .false.
   do_test_analytic_eigenvalues_eigenvectors = .false.
   do_test_frank_eigenvalues = .false.
#if defined(TEST_CHOLESKY)
   do_test_toeplitz_eigenvalues = .false.
#else
   do_test_toeplitz_eigenvalues = .true.
#endif
449

Andreas Marek's avatar
Andreas Marek committed
450
451
452
453
454
455
456
457
458
459
460
461
#endif /* TEST_MATRIX_TOEPLITZ */


#if defined(TEST_MATRIX_FRANK) && !defined(TEST_SOLVE_TRIDIAGONAL) && !defined(TEST_CHOLESKY)
   ! the random matrix can be used in allmost all tests; but for some no
   ! correctness checks have been implemented; do not allow these
   ! combinations
   ! FRANK + TEST_SOLVE_TRIDIAGONAL: we need a TOEPLITZ MATRIX
   ! FRANK + TEST_CHOLESKY: no correctness check yet implemented

   ! We also have to take care of special case in TEST_EIGENVECTORS
#if !defined(TEST_EIGENVECTORS)
Andreas Marek's avatar
WIP    
Andreas Marek committed
462
    !call prepare_matrix_frank(na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
463
464
465
466
467
468
469
470
471
472
473

    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
#ifndef TEST_HERMITIAN_MULTIPLY
    do_test_frank_eigenvalues = .true.
#endif
    do_test_toeplitz_eigenvalues = .false.

#else /* TEST_EIGENVECTORS */

    if (nev .ge. 1) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
474
     !call prepare_matrix_frank(na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526

    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
#ifndef TEST_HERMITIAN_MULTIPLY
    do_test_frank_eigenvalues = .true.
#endif
    do_test_toeplitz_eigenvalues = .false.
    do_test_numeric_residual = .false.
   else
    do_test_analytic_eigenvalues = .false.
    do_test_analytic_eigenvalues_eigenvectors = .false.
#ifndef TEST_HERMITIAN_MULTIPLY
    do_test_frank_eigenvalues = .true.
#endif
    do_test_toeplitz_eigenvalues = .false.
    do_test_numeric_residual = .false.

   endif

#endif /* TEST_EIGENVECTORS */
#endif /* (TEST_MATRIX_FRANK) */
#if defined(TEST_MATRIX_FRANK) && (defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_CHOLESKY))
#error "FRANK matrix is not allowed in this configuration"
#endif


#ifdef TEST_HERMITIAN_MULTIPLY
#ifdef TEST_REAL

#ifdef TEST_DOUBLE
   b(:,:) = 2.0_c_double * a(:,:)
   c(:,:) = 0.0_c_double
#else
   b(:,:) = 2.0_c_float * a(:,:)
   c(:,:) = 0.0_c_float
#endif

#endif /* TEST_REAL */

#ifdef TEST_COMPLEX

#ifdef TEST_DOUBLE
   b(:,:) = 2.0_c_double * a(:,:)
   c(:,:) = (0.0_c_double, 0.0_c_double)
#else
   b(:,:) = 2.0_c_float * a(:,:)
   c(:,:) = (0.0_c_float, 0.0_c_float)
#endif

#endif /* TEST_COMPLEX */

#endif /* TEST_HERMITIAN_MULTIPLY */
527

528
529
530
531
532
533
534
535
536
537
538
539
! if the test is used for (repeated) performacne tests, one might want to skip the checking
! of the results, which might be time-consuming and not necessary.
   if(skip_check_correctness) then
     do_test_numeric_residual = .false.
     do_test_numeric_residual_generalized = .false.
     do_test_analytic_eigenvalues = .false.
     do_test_analytic_eigenvalues_eigenvectors = .false.
     do_test_frank_eigenvalues = .false.
     do_test_toeplitz_eigenvalues = .false. 
     do_test_cholesky = .false.
   endif

540
541
542
543
544
545
546
547

#ifdef WITH_OPENMP
   threads_caller = omp_get_max_threads()
   if (myid == 0) then
     print *,"The calling program uses ",threads_caller," threads"
   endif
#endif

548
549
   e => elpa_allocate(error)
   assert_elpa_ok(error)
550

551
552
553
554
555
556
557
558
559
560
   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)
561
562

#ifdef WITH_MPI
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
#ifdef SPLIT_COMM_MYSELF
   call mpi_comm_split(MPI_COMM_WORLD,my_pcol,my_prow,mpi_comm_rows,mpierr)
   if (mpierr .ne. MPI_SUCCESS) then
     call MPI_ERROR_STRING(mpierr,mpierr_string, mpi_string_length, mpierr2)
     write(error_unit,*) "MPI ERROR occured during mpi_comm_split for row communicator: ", trim(mpierr_string)
     stop 1
   endif

   call mpi_comm_split(MPI_COMM_WORLD,my_prow,my_pcol,mpi_comm_cols, mpierr)
   if (mpierr .ne. MPI_SUCCESS) then
     call MPI_ERROR_STRING(mpierr,mpierr_string, mpi_string_length, mpierr2)
     write(error_unit,*) "MPI ERROR occured during mpi_comm_split for col communicator: ", trim(mpierr_string)
     stop 1
   endif

578
579
   call e%set("mpi_comm_parent", MPI_COMM_WORLD, error)
   assert_elpa_ok(error)
580
581
582
583
584
   call e%set("mpi_comm_rows", mpi_comm_rows, error)
   assert_elpa_ok(error)
   call e%set("mpi_comm_cols", mpi_comm_cols, error)
   assert_elpa_ok(error)
#else
585
586
587
588
589
590
   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)
591
#endif
592
#endif
593
594
595
#ifdef TEST_GENERALIZED_EIGENPROBLEM
   call e%set("blacs_context", my_blacs_ctxt, error)
   assert_elpa_ok(error)
596
#endif
Andreas Marek's avatar
Andreas Marek committed
597
   call e%set("timings",1,error)
598
599
600
   assert_elpa_ok(e%setup())

#ifdef TEST_SOLVER_1STAGE
Andreas Marek's avatar
Andreas Marek committed
601
   call e%set("solver", ELPA_SOLVER_1STAGE,error)
602
#else
Andreas Marek's avatar
Andreas Marek committed
603
   call e%set("solver", ELPA_SOLVER_2STAGE,error)
604
#endif
605
   assert_elpa_ok(error)
606

607
608
   call e%set("gpu", TEST_GPU, error)
   assert_elpa_ok(error)
609

Andreas Marek's avatar
Andreas Marek committed
610
#if TEST_QR_DECOMPOSITION == 1
611
612
613
614
   call e%set("qr", 1, error)
   assert_elpa_ok(error)
#endif

615
616
617
618
619
620
#ifdef WITH_OPENMP
   max_threads=omp_get_max_threads()
   call e%set("omp_threads", max_threads, error)
   assert_elpa_ok(error)
#endif

621
622
   if (myid == 0) print *, ""

623
#ifdef TEST_ALL_KERNELS
624
   do i = 0, elpa_option_cardinality(KERNEL_KEY)  ! kernels
Andreas Marek's avatar
Andreas Marek committed
625
626
627
628
629
     if (TEST_GPU .eq. 0) then
       kernel = elpa_option_enumerate(KERNEL_KEY, i)
       if (kernel .eq. ELPA_2STAGE_REAL_GPU) continue
       if (kernel .eq. ELPA_2STAGE_COMPLEX_GPU) continue
     endif
630
#endif
631
#ifdef TEST_KERNEL
632
     kernel = TEST_KERNEL
633
#endif
634
635

#ifdef TEST_SOLVER_2STAGE
636
     call e%set(KERNEL_KEY, kernel, error)
637
638
639
#ifdef TEST_KERNEL
     assert_elpa_ok(error)
#else
640
641
642
     if (error /= ELPA_OK) then
       cycle
     endif
643
     ! actually used kernel might be different if forced via environment variables
Andreas Marek's avatar
Andreas Marek committed
644
     call e%get(KERNEL_KEY, kernel, error)
645
     assert_elpa_ok(error)
646
647
648
649
#endif
     if (myid == 0) then
       print *, elpa_int_value_to_string(KERNEL_KEY, kernel) // " kernel"
     endif
650
651
#endif

Pavel Kus's avatar
Pavel Kus committed
652
653

! print all parameters
654
655
     call e%print_settings(error)
     assert_elpa_ok(error)
Pavel Kus's avatar
Pavel Kus committed
656

657
#ifdef TEST_ALL_KERNELS
658
     call e%timer_start(elpa_int_value_to_string(KERNEL_KEY, kernel))
659
#endif
660

661
     ! The actual solve step
Andreas Marek's avatar
Andreas Marek committed
662
663
664
665
#if defined(TEST_EIGENVECTORS)
#if TEST_QR_DECOMPOSITION == 1
     call e%timer_start("e%eigenvectors_qr()")
#else
666
     call e%timer_start("e%eigenvectors()")
Andreas Marek's avatar
Andreas Marek committed
667
#endif
Pavel Kus's avatar
Pavel Kus committed
668
669
#ifdef TEST_SCALAPACK_ALL
     call solve_scalapack_all(na, a, sc_desc, ev, z)
670
671
#elif TEST_SCALAPACK_PART
     call solve_scalapack_part(na, a, sc_desc, nev, ev, z)
672
     check_all_evals = .false. ! scalapack does not compute all eigenvectors
Pavel Kus's avatar
Pavel Kus committed
673
#else
674
     call e%eigenvectors(a, ev, z, error)
Pavel Kus's avatar
Pavel Kus committed
675
#endif
Andreas Marek's avatar
Andreas Marek committed
676
677
678
#if TEST_QR_DECOMPOSITION == 1
     call e%timer_stop("e%eigenvectors_qr()")
#else
679
     call e%timer_stop("e%eigenvectors()")
Andreas Marek's avatar
Andreas Marek committed
680
681
#endif
#endif /* TEST_EIGENVECTORS  */
682
683
684

#ifdef TEST_EIGENVALUES
     call e%timer_start("e%eigenvalues()")
685
     call e%eigenvalues(a, ev, error)
686
     call e%timer_stop("e%eigenvalues()")
687
#endif
688
689
690

#if defined(TEST_SOLVE_TRIDIAGONAL)
     call e%timer_start("e%solve_tridiagonal()")
691
     call e%solve_tridiagonal(d, sd, z, error)
692
     call e%timer_stop("e%solve_tridiagonal()")
693
694
695
     ev(:) = d(:)
#endif

696
697
698
#if defined(TEST_CHOLESKY)
     call e%timer_start("e%cholesky()")
     call e%cholesky(a, error)
699
     assert_elpa_ok(error)
700
701
702
     call e%timer_stop("e%cholesky()")
#endif

703
704
705
706
707
#if defined(TEST_HERMITIAN_MULTIPLY)
     call e%timer_start("e%hermitian_multiply()")
     call e%hermitian_multiply('F','F', na, a, b, na_rows, na_cols, c, na_rows, na_cols, error)
     call e%timer_stop("e%hermitian_multiply()")
#endif
Pavel Kus's avatar
Pavel Kus committed
708

709
710
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
     call e%timer_start("e%generalized_eigenvectors()")
711
712
713
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
     call e%timer_start("is_already_decomposed=.false.")
#endif
714
     call e%generalized_eigenvectors(a, b, ev, z, .false., error)
715
716
717
718
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
     call e%timer_stop("is_already_decomposed=.false.")
     a = as
     call e%timer_start("is_already_decomposed=.true.")
719
     call e%generalized_eigenvectors(a, b, ev, z, .true., error)
720
721
     call e%timer_stop("is_already_decomposed=.true.")
#endif
722
723
724
     call e%timer_stop("e%generalized_eigenvectors()")
#endif

725
726
     assert_elpa_ok(error)

727
#ifdef TEST_ALL_KERNELS
728
729
730
     call e%timer_stop(elpa_int_value_to_string(KERNEL_KEY, kernel))
#endif

731
     if (myid .eq. 0) then
732
#ifdef TEST_ALL_KERNELS
733
       call e%print_times(elpa_int_value_to_string(KERNEL_KEY, kernel))
734
735
#else /* TEST_ALL_KERNELS */

Andreas Marek's avatar
Andreas Marek committed
736
737
738
739
#if defined(TEST_EIGENVECTORS)
#if TEST_QR_DECOMPOSITION == 1
       call e%print_times("e%eigenvectors_qr()")
#else
740
       call e%print_times("e%eigenvectors()")
741
#endif
Andreas Marek's avatar
Andreas Marek committed
742
#endif
743
#ifdef TEST_EIGENVALUES
744
745
       call e%print_times("e%eigenvalues()")
#endif
746
747
#ifdef TEST_SOLVE_TRIDIAGONAL
       call e%print_times("e%solve_tridiagonal()")
748
#endif
749
750
#ifdef TEST_CHOLESKY
       call e%print_times("e%cholesky()")
751
#endif
752
753
754
#ifdef TEST_HERMITIAN_MULTIPLY
       call e%print_times("e%hermitian_multiply()")
#endif
755
756
757
#ifdef TEST_GENERALIZED_EIGENPROBLEM
      call e%print_times("e%generalized_eigenvectors()")
#endif
758
#endif /* TEST_ALL_KERNELS */
759
     endif
760

Andreas Marek's avatar
Andreas Marek committed
761
     if (do_test_analytic_eigenvalues) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
762
       !status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals, .false.)
Andreas Marek's avatar
Andreas Marek committed
763
764
765
766
       call check_status(status, myid)
     endif

     if (do_test_analytic_eigenvalues_eigenvectors) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
767
       !status = check_correctness_analytic(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals, .true.)
Andreas Marek's avatar
Andreas Marek committed
768
769
       call check_status(status, myid)
     endif
770

Andreas Marek's avatar
Andreas Marek committed
771
     if(do_test_numeric_residual) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
772
       !status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
773
774
775
776
       call check_status(status, myid)
     endif

     if (do_test_frank_eigenvalues) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
777
       !status = check_correctness_eigenvalues_frank(na, ev, z, myid)
Andreas Marek's avatar
Andreas Marek committed
778
       call check_status(status, myid)
779
     endif
780

Andreas Marek's avatar
Andreas Marek committed
781
     if (do_test_toeplitz_eigenvalues) then
782
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL)
Andreas Marek's avatar
WIP    
Andreas Marek committed
783
784
       !status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &
       !                                                subdiagonalElement, ev, z, myid)
Andreas Marek's avatar
Andreas Marek committed
785
       call check_status(status, myid)
786
#endif
Andreas Marek's avatar
Andreas Marek committed
787
     endif
788

Andreas Marek's avatar
Andreas Marek committed
789
     if (do_test_cholesky) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
790
       !status = check_correctness_cholesky(na, a, as, na_rows, sc_desc, myid )
Andreas Marek's avatar
Andreas Marek committed
791
792
       call check_status(status, myid)
     endif
793

Andreas Marek's avatar
Andreas Marek committed
794
795
#ifdef TEST_HERMITIAN_MULTIPLY
     if (do_test_hermitian_multiply) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
796
       !status = check_correctness_hermitian_multiply(na, a, b, c, na_rows, sc_desc, myid )
Andreas Marek's avatar
Andreas Marek committed
797
798
       call check_status(status, myid)
     endif
799
#endif
800

801
802
#ifdef TEST_GENERALIZED_EIGENPROBLEM
     if(do_test_numeric_residual_generalized) then
Andreas Marek's avatar
WIP    
Andreas Marek committed
803
804
       !status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, &
       !my_pcol, bs)
805
806
807
808
       call check_status(status, myid)
     endif
#endif

809
810
811
812
813
814
815
816
817

#ifdef WITH_OPENMP
     if (threads_caller .ne. omp_get_max_threads()) then
       if (myid .eq. 0) then
         print *, " ERROR! the number of OpenMP threads has not been restored correctly"
       endif
       status = 1
     endif
#endif
818
819
820
821
     if (myid == 0) then
       print *, ""
     endif

822
823
#ifdef TEST_ALL_KERNELS
     a(:,:) = as(:,:)
824
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
825
826
827
     d = ds
     sd = sds
#endif
828
   end do ! kernels
829
#endif
Andreas Marek's avatar
Andreas Marek committed
830

831
832
   call elpa_deallocate(e, error)
   assert_elpa_ok(error)
833
834
835
836
837

   deallocate(a)
   deallocate(as)
   deallocate(z)
   deallocate(ev)
838
839
840
841
#ifdef TEST_HERMITIAN_MULTIPLY
   deallocate(b)
   deallocate(c)
#endif
842
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
843
844
845
   deallocate(d, ds)
   deallocate(sd, sds)
#endif
846
847
848
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
  deallocate(b, bs)
#endif
849

850
851
852
853
#ifdef TEST_ALL_LAYOUTS
   end do ! factors
   end do ! layouts
#endif
854
855
   call elpa_uninit(error)
   assert_elpa_ok(error)
856

857
858
859
860
861
862
#ifdef WITH_MPI
   call blacs_gridexit(my_blacs_ctxt)
   call mpi_finalize(mpierr)
#endif
   call exit(status)

863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
   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

878
end program