test.F90 27.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
Andreas Marek committed
105
106
107
108
109
110
111
#ifdef HAVE_64BIT_INTEGER_SUPPORT
#define TEST_INT_TYPE integer(kind=c_int64_t)
#define INT_TYPE c_int64_t
#else
#define TEST_INT_TYPE integer(kind=c_int32_t)
#define INT_TYPE c_int32_t
#endif
112
113
114
#include "assert.h"

program test
Andreas Marek's avatar
Andreas Marek committed
115
116
   use elpa !, only : elpa_allocate, elpa_init, elpa_allocate, elpa_t, ELPA_SOLVER_1STAGE, elpa_deallocate, &
            !       elpa_uninit, 
117

Andreas Marek's avatar
WIP    
Andreas Marek committed
118
   !use test_util
Andreas Marek's avatar
Andreas Marek committed
119
120
   use test_setup_mpi
   use test_prepare_matrix
121
   use test_read_input_parameters
Andreas Marek's avatar
Andreas Marek committed
122
123
124
   use test_blacs_infrastructure
   use test_check_correctness
   use test_analytic
125
#ifdef WITH_SCALAPACK_TESTS
Andreas Marek's avatar
Andreas Marek committed
126
   use test_scalapack
127
#endif
128

129
#ifdef HAVE_REDIRECT
Andreas Marek's avatar
Andreas Marek committed
130
   use test_redirect
131
132
#endif
#ifdef WITH_OPENMP
Andreas Marek's avatar
Andreas Marek committed
133
   use omp_lib
134
#endif
135
   use precision_for_tests
136

137
138
139
   implicit none

   ! matrix dimensions
Andreas Marek's avatar
Andreas Marek committed
140
   TEST_INT_TYPE      :: na, nev, nblk
141
142

   ! mpi
Andreas Marek's avatar
Andreas Marek committed
143
144
145
146
147
   TEST_INT_TYPE      :: myid, nprocs
   TEST_INT_TYPE      :: na_cols, na_rows  ! local matrix size
   TEST_INT_TYPE      :: np_cols, np_rows  ! number of MPI processes per column/row
   TEST_INT_TYPE      :: my_prow, my_pcol  ! local MPI task position (my_prow, my_pcol) in the grid (0..np_cols -1, 0..np_rows -1)
   TEST_INT_TYPE      :: mpierr
148
149

   ! blacs
Andreas Marek's avatar
Andreas Marek committed
150
   TEST_INT_TYPE     :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
151
152

   ! The Matrix
153
   MATRIX_TYPE, allocatable    :: a(:,:), as(:,:)
154
#if defined(TEST_HERMITIAN_MULTIPLY)
155
   MATRIX_TYPE, allocatable    :: b(:,:), c(:,:)
156
157
158
#endif
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
   MATRIX_TYPE, allocatable    :: b(:,:), bs(:,:)
159
#endif
160
   ! eigenvectors
161
   MATRIX_TYPE, allocatable    :: z(:,:)
162
   ! eigenvalues
163
   EV_TYPE, allocatable        :: ev(:)
164

165
   logical                     :: check_all_evals, skip_check_correctness
166

167
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
168
169
   EV_TYPE, allocatable        :: d(:), sd(:), ds(:), sds(:)
   EV_TYPE                     :: diagonalELement, subdiagonalElement
170
#endif
171

Andreas Marek's avatar
Andreas Marek committed
172
173
   TEST_INT_TYPE               :: status
   integer(kind=c_int)         :: error_elpa
174

175
176
   type(output_t)              :: write_to_file
   class(elpa_t), pointer      :: e
177
#ifdef TEST_ALL_KERNELS
Andreas Marek's avatar
Andreas Marek committed
178
   TEST_INT_TYPE      :: i
179
180
181
#endif
#ifdef TEST_ALL_LAYOUTS
   character(len=1), parameter :: layouts(2) = [ 'C', 'R' ]
Andreas Marek's avatar
Andreas Marek committed
182
   TEST_INT_TYPE      :: i_layout
183
#endif
Andreas Marek's avatar
Andreas Marek committed
184
   integer(kind=c_int)         :: kernel
185
   character(len=1)            :: layout
186
187
   logical                     :: do_test_numeric_residual, do_test_numeric_residual_generalized, &
                                  do_test_analytic_eigenvalues, &
Andreas Marek's avatar
Andreas Marek committed
188
189
190
191
192
                                  do_test_analytic_eigenvalues_eigenvectors,   &
                                  do_test_frank_eigenvalues,  &
                                  do_test_toeplitz_eigenvalues, do_test_cholesky,   &
                                  do_test_hermitian_multiply

193
#ifdef WITH_OPENMP
Andreas Marek's avatar
Andreas Marek committed
194
   TEST_INT_TYPE      :: max_threads, threads_caller
195
196
#endif

197
#ifdef SPLIT_COMM_MYSELF
Andreas Marek's avatar
Andreas Marek committed
198
   TEST_INT_TYPE      :: mpi_comm_rows, mpi_comm_cols, mpi_string_length, mpierr2
199
200
201
   character(len=MPI_MAX_ERROR_STRING) :: mpierr_string
#endif

202
   call read_input_parameters_traditional(na, nev, nblk, write_to_file, skip_check_correctness)
Andreas Marek's avatar
Andreas Marek committed
203
   call setup_mpi(myid, nprocs)
204
205
206
207
#ifdef HAVE_REDIRECT
#ifdef WITH_MPI
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
     call redirect_stdout(myid)
208
#endif
209
#endif
210

211
   check_all_evals = .true.
212

Andreas Marek's avatar
Andreas Marek committed
213
214

   do_test_numeric_residual = .false.
215
   do_test_numeric_residual_generalized = .false.
Andreas Marek's avatar
Andreas Marek committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
   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
231
232
233
234
235
   if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
     print *, "ELPA API version not supported"
     stop 1
   endif

236
237
238
239
240
   if (myid == 0) then
     print '((a,i0))', 'Program ' // TEST_CASE
     print *, ""
   endif

241
#ifdef TEST_ALL_LAYOUTS
242
   do i_layout = 1, size(layouts)               ! layouts
243
     layout = layouts(i_layout)
244
     do np_cols = 1, nprocs                     ! factors
245
246
247
248
249
       if (mod(nprocs,np_cols) /= 0 ) then
         cycle
       endif
#else
   layout = 'C'
250
251
252
   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
      if(mod(nprocs,np_cols) == 0 ) exit
   enddo
253
#endif
254
255

   np_rows = nprocs/np_cols
256
   assert(nprocs == np_rows * np_cols)
257

258
259
260
261
   if (myid == 0) then
     print '((a,i0))', 'Matrix size: ', na
     print '((a,i0))', 'Num eigenvectors: ', nev
     print '((a,i0))', 'Blocksize: ', nblk
262
#ifdef WITH_MPI
263
264
     print '((a,i0))', 'Num MPI proc: ', nprocs
     print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
265
     print '(a)',      'Process layout: ' // layout
266
#endif
267
268
269
     print *,''
   endif

Andreas Marek's avatar
Andreas Marek committed
270
#if TEST_QR_DECOMPOSITION == 1
271
272
273
274
275
276

#if TEST_GPU == 1
#ifdef WITH_MPI
     call mpi_finalize(mpierr)
#endif
     stop 77
Andreas Marek's avatar
Andreas Marek committed
277
#endif /* TEST_GPU */
278
279
280
281
282
283
284
285
286
287
288
289
290
   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
291
292
#endif /* TEST_QR_DECOMPOSITION */

293

Andreas Marek's avatar
Andreas Marek committed
294
295
   call set_up_blacsgrid(mpi_comm_world, np_rows, np_cols, layout, &
                         my_blacs_ctxt, my_prow, my_pcol)
296

Andreas Marek's avatar
Andreas Marek committed
297
298
   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)
299

Pavel Kus's avatar
Pavel Kus committed
300
301
   allocate(a (na_rows,na_cols))
   allocate(as(na_rows,na_cols))
302
303
304
   allocate(z (na_rows,na_cols))
   allocate(ev(na))

305
306
307
308
309
#ifdef TEST_HERMITIAN_MULTIPLY
   allocate(b (na_rows,na_cols))
   allocate(c (na_rows,na_cols))
#endif

310
311
312
313
314
#ifdef TEST_GENERALIZED_EIGENPROBLEM
   allocate(b (na_rows,na_cols))
   allocate(bs (na_rows,na_cols))
#endif

315
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
316
317
318
319
   allocate(d (na), ds(na))
   allocate(sd (na), sds(na))
#endif

320
321
322
323
   a(:,:) = 0.0
   z(:,:) = 0.0
   ev(:) = 0.0

Andreas Marek's avatar
Andreas Marek committed
324
325
326
327
328
#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
329
   ! RANDOM + TEST_CHOLESKY: wee need SPD matrix
Andreas Marek's avatar
Andreas Marek committed
330
331
332
333
   ! 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
Andreas Marek committed
334
    call prepare_matrix_random(na, myid, sc_desc, a, z, as)
Andreas Marek's avatar
Andreas Marek committed
335
336
#else /* TEST_EIGENVECTORS */
    if (nev .ge. 1) then
Andreas Marek's avatar
Andreas Marek committed
337
      call prepare_matrix_random(na, myid, sc_desc, a, z, as)
Andreas Marek's avatar
Andreas Marek committed
338
#ifndef TEST_HERMITIAN_MULTIPLY
Andreas Marek's avatar
Andreas Marek committed
339
      do_test_numeric_residual = .true.
Andreas Marek's avatar
Andreas Marek committed
340
#endif
341
   else
Andreas Marek's avatar
Andreas Marek committed
342
343
344
345
346
     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)
347
#endif
Andreas Marek's avatar
Andreas Marek committed
348
349
350
     stop 77
   endif
#endif /* TEST_EIGENVECTORS */
351
352
353
354
    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
355
#endif /* (TEST_MATRIX_RANDOM) */
356
357

#if defined(TEST_MATRIX_RANDOM) && defined(TEST_CHOLESKY)
Andreas Marek's avatar
Andreas Marek committed
358
359
     call prepare_matrix_random_spd(na, myid, sc_desc, a, z, as, &
                 nblk, np_rows, np_cols, my_prow, my_pcol)
360
361
362
363
364
365
    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 */

366
#if defined(TEST_MATRIX_RANDOM) && defined(TEST_GENERALIZED_EIGENPROBLEM)
Andreas Marek's avatar
Andreas Marek committed
367
368
369
   ! 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)
370
371
372
373
374
375
376
377
    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 */

378
#if defined(TEST_MATRIX_RANDOM) && (defined(TEST_SOLVE_TRIDIAGONAL) || defined(TEST_EIGENVALUES))
Andreas Marek's avatar
Andreas Marek committed
379
#error "Random matrix is not allowed in this configuration"
380
#endif
381

Andreas Marek's avatar
Andreas Marek committed
382
383
384
385
386
387
#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
388

Andreas Marek's avatar
Andreas Marek committed
389
   call prepare_matrix_analytic(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
390
   as(:,:) = a
391

Andreas Marek's avatar
Andreas Marek committed
392
393
394
395
   do_test_numeric_residual = .false.
   do_test_analytic_eigenvalues_eigenvectors = .false.
#ifndef TEST_HERMITIAN_MULTIPLY
   do_test_analytic_eigenvalues = .true.
396
#endif
Andreas Marek's avatar
Andreas Marek committed
397
398
399
#if defined(TEST_EIGENVECTORS)
   if (nev .ge. 1) then
     do_test_analytic_eigenvalues_eigenvectors = .true.
Pavel Kus's avatar
Pavel Kus committed
400
     do_test_analytic_eigenvalues = .false.
Andreas Marek's avatar
Andreas Marek committed
401
402
403
   else
     do_test_analytic_eigenvalues_eigenvectors = .false.
   endif
404
#endif
Andreas Marek's avatar
Andreas Marek committed
405
406
   do_test_frank_eigenvalues = .false.
   do_test_toeplitz_eigenvalues = .false.
407
#endif /* TEST_MATRIX_ANALYTIC */
Andreas Marek's avatar
Andreas Marek committed
408
409
410
#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
411

Andreas Marek's avatar
Andreas Marek committed
412
413
#if defined(TEST_MATRIX_TOEPLITZ)
   ! The Toeplitz matrix works in each test
Andreas Marek's avatar
Andreas Marek committed
414
415
416
417
418
419
420
#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
421

422
! actually we test cholesky for diagonal matrix only
423
424
#if defined(TEST_CHOLESKY)
#ifdef TEST_SINGLE
425
426
  diagonalElement = (2.546_c_float, 0.0_c_float)
  subdiagonalElement =  (0.0_c_float, 0.0_c_float)
427
#else
428
429
  diagonalElement = (2.546_c_double, 0.0_c_double)
  subdiagonalElement =  (0.0_c_double, 0.0_c_double)
430
#endif
Andreas Marek's avatar
Andreas Marek committed
431
432
#endif /* TEST_CHOLESKY */

Andreas Marek's avatar
Andreas Marek committed
433
434
435
   call prepare_matrix_toeplitz(na, diagonalElement, subdiagonalElement, &
                                d, sd, ds, sds, a, as, nblk, np_rows, &
                                np_cols, my_prow, my_pcol)
436
437


Andreas Marek's avatar
Andreas Marek committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
   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
455

Andreas Marek's avatar
Andreas Marek committed
456
457
458
459
460
461
462
463
464
465
466
467
#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
Andreas Marek committed
468
    call prepare_matrix_frank(na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
469
470
471
472
473
474
475
476
477
478
479

    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
Andreas Marek committed
480
      call prepare_matrix_frank(na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
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
527
528
529
530
531
532

    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 */
533

534
535
536
537
538
539
540
541
542
543
544
545
! 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

546
547
548
549
550
551
552
553

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

Andreas Marek's avatar
Andreas Marek committed
554
555
   e => elpa_allocate(error_elpa)
   assert_elpa_ok(error_elpa)
556

Andreas Marek's avatar
Andreas Marek committed
557
558
559
560
561
562
563
564
565
566
   call e%set("na", int(na,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("nev", int(nev,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("local_nrows", int(na_rows,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("local_ncols", int(na_cols,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("nblk", int(nblk,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
567
568

#ifdef WITH_MPI
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
#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

Andreas Marek's avatar
Andreas Marek committed
584
585
586
587
588
589
   call e%set("mpi_comm_parent", int(MPI_COMM_WORLD,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("mpi_comm_rows", int(mpi_comm_rows,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("mpi_comm_cols", int(mpi_comm_cols,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
590
#else
Andreas Marek's avatar
Andreas Marek committed
591
592
593
594
595
596
   call e%set("mpi_comm_parent", int(MPI_COMM_WORLD,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("process_row", int(my_prow,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
   call e%set("process_col", int(my_pcol,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
597
#endif
598
#endif
599
#ifdef TEST_GENERALIZED_EIGENPROBLEM
Andreas Marek's avatar
Andreas Marek committed
600
601
   call e%set("blacs_context", int(my_blacs_ctxt,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
602
#endif
Andreas Marek's avatar
Andreas Marek committed
603
   call e%set("timings", 1_ik, error_elpa)
604
605
606
   assert_elpa_ok(e%setup())

#ifdef TEST_SOLVER_1STAGE
Andreas Marek's avatar
Andreas Marek committed
607
   call e%set("solver", ELPA_SOLVER_1STAGE, error_elpa)
608
#else
Andreas Marek's avatar
Andreas Marek committed
609
   call e%set("solver", ELPA_SOLVER_2STAGE, error_elpa)
610
#endif
Andreas Marek's avatar
Andreas Marek committed
611
   assert_elpa_ok(error_elpa)
612

Andreas Marek's avatar
Andreas Marek committed
613
614
   call e%set("gpu", TEST_GPU, error_elpa)
   assert_elpa_ok(error_elpa)
615

Andreas Marek's avatar
Andreas Marek committed
616
#if TEST_QR_DECOMPOSITION == 1
Andreas Marek's avatar
Andreas Marek committed
617
618
   call e%set("qr", 1_ik, error_elpa)
   assert_elpa_ok(error_elpa)
619
620
#endif

621
622
#ifdef WITH_OPENMP
   max_threads=omp_get_max_threads()
Andreas Marek's avatar
Andreas Marek committed
623
624
   call e%set("omp_threads", int(max_threads,kind=c_int), error_elpa)
   assert_elpa_ok(error_elpa)
625
626
#endif

627
628
   if (myid == 0) print *, ""

629
#ifdef TEST_ALL_KERNELS
630
   do i = 0, elpa_option_cardinality(KERNEL_KEY)  ! kernels
Andreas Marek's avatar
Andreas Marek committed
631
     if (TEST_GPU .eq. 0) then
Andreas Marek's avatar
Andreas Marek committed
632
       kernel = elpa_option_enumerate(KERNEL_KEY, int(i,kind=c_int))
Andreas Marek's avatar
Andreas Marek committed
633
634
635
       if (kernel .eq. ELPA_2STAGE_REAL_GPU) continue
       if (kernel .eq. ELPA_2STAGE_COMPLEX_GPU) continue
     endif
636
#endif
637
#ifdef TEST_KERNEL
638
     kernel = TEST_KERNEL
639
#endif
640
641

#ifdef TEST_SOLVER_2STAGE
Andreas Marek's avatar
Andreas Marek committed
642
     call e%set(KERNEL_KEY, kernel, error_elpa)
643
#ifdef TEST_KERNEL
Andreas Marek's avatar
Andreas Marek committed
644
     assert_elpa_ok(error_elpa)
645
#else
Andreas Marek's avatar
Andreas Marek committed
646
     if (error_elpa /= ELPA_OK) then
647
648
       cycle
     endif
649
     ! actually used kernel might be different if forced via environment variables
Andreas Marek's avatar
Andreas Marek committed
650
651
     call e%get(KERNEL_KEY, kernel, error_elpa)
     assert_elpa_ok(error_elpa)
652
653
654
655
#endif
     if (myid == 0) then
       print *, elpa_int_value_to_string(KERNEL_KEY, kernel) // " kernel"
     endif
656
657
#endif

Pavel Kus's avatar
Pavel Kus committed
658
659

! print all parameters
Andreas Marek's avatar
Andreas Marek committed
660
661
     call e%print_settings(error_elpa)
     assert_elpa_ok(error_elpa)
Pavel Kus's avatar
Pavel Kus committed
662

663
#ifdef TEST_ALL_KERNELS
664
     call e%timer_start(elpa_int_value_to_string(KERNEL_KEY, kernel))
665
#endif
666

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

#ifdef TEST_EIGENVALUES
     call e%timer_start("e%eigenvalues()")
Andreas Marek's avatar
Andreas Marek committed
691
     call e%eigenvalues(a, ev, error_elpa)
692
     call e%timer_stop("e%eigenvalues()")
693
#endif
694
695
696

#if defined(TEST_SOLVE_TRIDIAGONAL)
     call e%timer_start("e%solve_tridiagonal()")
Andreas Marek's avatar
Andreas Marek committed
697
     call e%solve_tridiagonal(d, sd, z, error_elpa)
698
     call e%timer_stop("e%solve_tridiagonal()")
699
700
701
     ev(:) = d(:)
#endif

702
703
#if defined(TEST_CHOLESKY)
     call e%timer_start("e%cholesky()")
Andreas Marek's avatar
Andreas Marek committed
704
705
     call e%cholesky(a, error_elpa)
     assert_elpa_ok(error_elpa)
706
707
708
     call e%timer_stop("e%cholesky()")
#endif

709
710
#if defined(TEST_HERMITIAN_MULTIPLY)
     call e%timer_start("e%hermitian_multiply()")
Andreas Marek's avatar
Andreas Marek committed
711
712
713
     call e%hermitian_multiply('F','F', int(na,kind=c_int), a, b, int(na_rows,kind=c_int), &
                               int(na_cols,kind=c_int), c, int(na_rows,kind=c_int),        &
                               int(na_cols,kind=c_int), error_elpa)
714
715
     call e%timer_stop("e%hermitian_multiply()")
#endif
Pavel Kus's avatar
Pavel Kus committed
716

717
718
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
     call e%timer_start("e%generalized_eigenvectors()")
719
720
721
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
     call e%timer_start("is_already_decomposed=.false.")
#endif
Andreas Marek's avatar
Andreas Marek committed
722
     call e%generalized_eigenvectors(a, b, ev, z, .false., error_elpa)
723
724
725
726
#if defined(TEST_GENERALIZED_DECOMP_EIGENPROBLEM)
     call e%timer_stop("is_already_decomposed=.false.")
     a = as
     call e%timer_start("is_already_decomposed=.true.")
Andreas Marek's avatar
Andreas Marek committed
727
     call e%generalized_eigenvectors(a, b, ev, z, .true., error_elpa)
728
729
     call e%timer_stop("is_already_decomposed=.true.")
#endif
730
731
732
     call e%timer_stop("e%generalized_eigenvectors()")
#endif

Andreas Marek's avatar
Andreas Marek committed
733
     assert_elpa_ok(error_elpa)
734

735
#ifdef TEST_ALL_KERNELS
736
737
738
     call e%timer_stop(elpa_int_value_to_string(KERNEL_KEY, kernel))
#endif

739
     if (myid .eq. 0) then
740
#ifdef TEST_ALL_KERNELS
741
       call e%print_times(elpa_int_value_to_string(KERNEL_KEY, kernel))
742
743
#else /* TEST_ALL_KERNELS */

Andreas Marek's avatar
Andreas Marek committed
744
745
746
747
#if defined(TEST_EIGENVECTORS)
#if TEST_QR_DECOMPOSITION == 1
       call e%print_times("e%eigenvectors_qr()")
#else
748
       call e%print_times("e%eigenvectors()")
749
#endif
Andreas Marek's avatar
Andreas Marek committed
750
#endif
751
#ifdef TEST_EIGENVALUES
752
753
       call e%print_times("e%eigenvalues()")
#endif
754
755
#ifdef TEST_SOLVE_TRIDIAGONAL
       call e%print_times("e%solve_tridiagonal()")
756
#endif
757
758
#ifdef TEST_CHOLESKY
       call e%print_times("e%cholesky()")
759
#endif
760
761
762
#ifdef TEST_HERMITIAN_MULTIPLY
       call e%print_times("e%hermitian_multiply()")
#endif
763
764
765
#ifdef TEST_GENERALIZED_EIGENPROBLEM
      call e%print_times("e%generalized_eigenvectors()")
#endif
766
#endif /* TEST_ALL_KERNELS */
767
     endif
768

Andreas Marek's avatar
Andreas Marek committed
769
     if (do_test_analytic_eigenvalues) then
Andreas Marek's avatar
Andreas Marek committed
770
       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
771
772
773
774
       call check_status(status, myid)
     endif

     if (do_test_analytic_eigenvalues_eigenvectors) then
Andreas Marek's avatar
Andreas Marek committed
775
       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
776
777
       call check_status(status, myid)
     endif
778

Andreas Marek's avatar
Andreas Marek committed
779
     if(do_test_numeric_residual) then
Andreas Marek's avatar
Andreas Marek committed
780
       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
781
782
783
784
       call check_status(status, myid)
     endif

     if (do_test_frank_eigenvalues) then
Andreas Marek's avatar
Andreas Marek committed
785
       status = check_correctness_eigenvalues_frank(na, ev, z, myid)
Andreas Marek's avatar
Andreas Marek committed
786
       call check_status(status, myid)
787
     endif
788

Andreas Marek's avatar
Andreas Marek committed
789
     if (do_test_toeplitz_eigenvalues) then
790
#if defined(TEST_EIGENVALUES) || defined(TEST_SOLVE_TRIDIAGONAL)
Andreas Marek's avatar
Andreas Marek committed
791
792
       status = check_correctness_eigenvalues_toeplitz(na, diagonalElement, &
                                                       subdiagonalElement, ev, z, myid)
Andreas Marek's avatar
Andreas Marek committed
793
       call check_status(status, myid)
794
#endif
Andreas Marek's avatar
Andreas Marek committed
795
     endif
796

Andreas Marek's avatar
Andreas Marek committed
797
     if (do_test_cholesky) then
Andreas Marek's avatar
Andreas Marek committed
798
       status = check_correctness_cholesky(na, a, as, na_rows, sc_desc, myid )
Andreas Marek's avatar
Andreas Marek committed
799
800
       call check_status(status, myid)
     endif
801

Andreas Marek's avatar
Andreas Marek committed
802
803
#ifdef TEST_HERMITIAN_MULTIPLY
     if (do_test_hermitian_multiply) then
Andreas Marek's avatar
Andreas Marek committed
804
       status = check_correctness_hermitian_multiply(na, a, b, c, na_rows, sc_desc, myid )
Andreas Marek's avatar
Andreas Marek committed
805
806
       call check_status(status, myid)
     endif
807
#endif
808

809
810
#ifdef TEST_GENERALIZED_EIGENPROBLEM
     if(do_test_numeric_residual_generalized) then
Andreas Marek's avatar
Andreas Marek committed
811
812
       status = check_correctness_evp_numeric_residuals(na, nev, as, z, ev, sc_desc, nblk, myid, np_rows,np_cols, my_prow, &
       my_pcol, bs)
813
814
815
816
       call check_status(status, myid)
     endif
#endif

817
818
819
820
821
822
823
824
825

#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
826
827
828
829
     if (myid == 0) then
       print *, ""
     endif

830
831
#ifdef TEST_ALL_KERNELS
     a(:,:) = as(:,:)
832
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
833
834
835
     d = ds
     sd = sds
#endif
836
   end do ! kernels
837
#endif
Andreas Marek's avatar
Andreas Marek committed
838

Andreas Marek's avatar
Andreas Marek committed
839
840
   call elpa_deallocate(e, error_elpa)
   assert_elpa_ok(error_elpa)
841
842
843
844
845

   deallocate(a)
   deallocate(as)
   deallocate(z)
   deallocate(ev)
846
847
848
849
#ifdef TEST_HERMITIAN_MULTIPLY
   deallocate(b)
   deallocate(c)
#endif
850
#if defined(TEST_MATRIX_TOEPLITZ) || defined(TEST_MATRIX_FRANK)
851
852
853
   deallocate(d, ds)
   deallocate(sd, sds)
#endif
854
855
856
#if defined(TEST_GENERALIZED_EIGENPROBLEM)
  deallocate(b, bs)
#endif
857

858
859
860
861
#ifdef TEST_ALL_LAYOUTS
   end do ! factors
   end do ! layouts
#endif
Andreas Marek's avatar
Andreas Marek committed
862
863
   call elpa_uninit(error_elpa)
   assert_elpa_ok(error_elpa)
864

865
866
867
868
869
870
#ifdef WITH_MPI
   call blacs_gridexit(my_blacs_ctxt)
   call mpi_finalize(mpierr)
#endif
   call exit(status)

871
872
873
874
   contains

     subroutine check_status(status, myid)
       implicit none
Andreas Marek's avatar
Andreas Marek committed
875
876
       TEST_INT_TYPE, intent(in) :: status, myid
       TEST_INT_TYPE             :: mpierr
877
878
879
880
881
882
883
884
885
       if (status /= 0) then
         if (myid == 0) print *, "Result incorrect!"
#ifdef WITH_MPI
         call mpi_finalize(mpierr)
#endif
         call exit(status)
       endif
     end subroutine

886
end program