test_skewsymmetric.F90 13.6 KB
Newer Older
Andreas Marek's avatar
Andreas Marek committed
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
!    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
48
49
! Define TEST_NVIDIA_GPU \in [0, 1]
! Define TEST_AMD_GPU \in [0, 1]
Andreas Marek's avatar
Andreas Marek committed
50
51
52
53
54
55
56
57
58
59
60
61
! Define either TEST_ALL_KERNELS or a TEST_KERNEL \in [any valid kernel]

#if !(defined(TEST_REAL) ^ defined(TEST_COMPLEX))
error: define exactly one of TEST_REAL or TEST_COMPLEX
#endif

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

#ifdef TEST_SINGLE
#  define EV_TYPE real(kind=C_FLOAT)
Andreas Marek's avatar
Andreas Marek committed
62
63
#  define EV_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
#  define MATRIX_TYPE_COMPLEX complex(kind=C_FLOAT_COMPLEX)
Andreas Marek's avatar
Andreas Marek committed
64
65
66
67
68
69
#  ifdef TEST_REAL
#    define MATRIX_TYPE real(kind=C_FLOAT)
#  else
#    define MATRIX_TYPE complex(kind=C_FLOAT_COMPLEX)
#  endif
#else
Andreas Marek's avatar
Andreas Marek committed
70
71
#  define MATRIX_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
#  define EV_TYPE_COMPLEX complex(kind=C_DOUBLE_COMPLEX)
Andreas Marek's avatar
Andreas Marek committed
72
73
74
75
76
77
78
79
80
81
82
83
84
#  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

#ifdef TEST_REAL
#  define AUTOTUNE_DOMAIN ELPA_AUTOTUNE_DOMAIN_REAL
#else
#  define AUTOTUNE_DOMAIN ELPA_AUTOTUNE_DOMAIN_COMPLEX
#endif
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#ifdef HAVE_64BIT_INTEGER_MATH_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
#ifdef HAVE_64BIT_INTEGER_MPI_SUPPORT
#define TEST_INT_MPI_TYPE integer(kind=c_int64_t)
#define INT_MPI_TYPE c_int64_t
#else
#define TEST_INT_MPI_TYPE integer(kind=c_int32_t)
#define INT_MPI_TYPE c_int32_t
#endif
99
100
101
102
103
104
105

#define TEST_GPU 0
#if (TEST_NVIDIA_GPU == 1) || (TEST_AMD_GPU == 1)
#undef TEST_GPU
#define TEST_GPU 1
#endif

Andreas Marek's avatar
Andreas Marek committed
106
107
108
109
110
#include "assert.h"

program test
   use elpa

111
   !use test_util
Andreas Marek's avatar
Andreas Marek committed
112
113
114
115
116
   use test_setup_mpi
   use test_prepare_matrix
   use test_read_input_parameters
   use test_blacs_infrastructure
   use test_check_correctness
117
   use precision_for_tests
Andreas Marek's avatar
Andreas Marek committed
118
119
120
121
122
123
124
125
   use iso_fortran_env

#ifdef HAVE_REDIRECT
   use test_redirect
#endif
   implicit none

   ! matrix dimensions
126
   TEST_INT_TYPE                          :: na, nev, nblk
Andreas Marek's avatar
Andreas Marek committed
127
128

   ! mpi
129
130
131
132
133
   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_MPI_TYPE                      :: mpierr
Andreas Marek's avatar
Andreas Marek committed
134
135

   ! blacs
136
   character(len=1)                 :: layout
137
   TEST_INT_TYPE                          :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
Andreas Marek's avatar
Andreas Marek committed
138
139

   ! The Matrix
140
141
   MATRIX_TYPE, allocatable         :: a_skewsymmetric(:,:), as_skewsymmetric(:,:)
   MATRIX_TYPE_COMPLEX, allocatable :: a_complex(:,:), as_complex(:,:)
Andreas Marek's avatar
Andreas Marek committed
142
   ! eigenvectors
143
144
   MATRIX_TYPE, allocatable         :: z_skewsymmetric(:,:)
   MATRIX_TYPE_COMPLEX, allocatable :: z_complex(:,:)
Andreas Marek's avatar
Andreas Marek committed
145
   ! eigenvalues
146
   EV_TYPE, allocatable             :: ev_skewsymmetric(:), ev_complex(:)
Andreas Marek's avatar
Andreas Marek committed
147

148
149
   TEST_INT_TYPE                    :: status, i, j
   integer(kind=c_int)              :: error_elpa
Andreas Marek's avatar
Andreas Marek committed
150

151
152
   type(output_t)                   :: write_to_file
   class(elpa_t), pointer           :: e_complex, e_skewsymmetric
153
           
Andreas Marek's avatar
Andreas Marek committed
154
155
156
157
158
159
160
161
162
163
164
165
166
   call read_input_parameters(na, nev, nblk, write_to_file)
   call setup_mpi(myid, nprocs)
#ifdef HAVE_REDIRECT
#ifdef WITH_MPI
   call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
   call redirect_stdout(myid)
#endif
#endif

   if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
     print *, "ELPA API version not supported"
     stop 1
   endif
Carolin Penke's avatar
Carolin Penke committed
167
! 
Andreas Marek's avatar
Andreas Marek committed
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
   layout = 'C'
   do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
      if(mod(nprocs,np_cols) == 0 ) exit
   enddo
   np_rows = nprocs/np_cols
   assert(nprocs == np_rows * np_cols)

   if (myid == 0) then
     print '((a,i0))', 'Matrix size: ', na
     print '((a,i0))', 'Num eigenvectors: ', nev
     print '((a,i0))', 'Blocksize: ', nblk
#ifdef WITH_MPI
     print '((a,i0))', 'Num MPI proc: ', nprocs
     print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
     print '(a)',      'Process layout: ' // layout
#endif
     print *,''
   endif

187
188
189
   call set_up_blacsgrid(int(mpi_comm_world,kind=BLAS_KIND), np_rows, &
                             np_cols, layout, &
                             my_blacs_ctxt, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
190
191
192
193

   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)

194
195
   allocate(a_skewsymmetric (na_rows,na_cols))
   allocate(as_skewsymmetric(na_rows,na_cols))
Carolin Penke's avatar
Carolin Penke committed
196
   allocate(z_skewsymmetric (na_rows,2*na_cols))
197
   allocate(ev_skewsymmetric(na))
Andreas Marek's avatar
Andreas Marek committed
198

199
200
201
   a_skewsymmetric(:,:) = 0.0
   z_skewsymmetric(:,:) = 0.0
   ev_skewsymmetric(:) = 0.0
Andreas Marek's avatar
Andreas Marek committed
202

Carolin Penke's avatar
Carolin Penke committed
203
204
   call prepare_matrix_random(na, myid, sc_desc, a_skewsymmetric, &
   z_skewsymmetric(:,1:na_cols), as_skewsymmetric, is_skewsymmetric=1)
Carolin Penke's avatar
Carolin Penke committed
205
   
206
   !call MPI_BARRIER(MPI_COMM_WORLD, mpierr)  
207
   as_skewsymmetric(:,:) = a_skewsymmetric(:,:)
Carolin Penke's avatar
Carolin Penke committed
208
   
Andreas Marek's avatar
Andreas Marek committed
209

210
211
212
213
214
   ! prepare the complex matrix for the "brute force" case
   allocate(a_complex (na_rows,na_cols))
   allocate(as_complex(na_rows,na_cols))
   allocate(z_complex (na_rows,na_cols))
   allocate(ev_complex(na))
Andreas Marek's avatar
Andreas Marek committed
215

Carolin Penke's avatar
Carolin Penke committed
216
217
218
   a_complex(1:na_rows,1:na_cols) = 0.0
   z_complex(1:na_rows,1:na_cols) = 0.0
   as_complex(1:na_rows,1:na_cols) = 0.0
Carolin Penke's avatar
Carolin Penke committed
219
   
Andreas Marek's avatar
Andreas Marek committed
220

Carolin Penke's avatar
Carolin Penke committed
221
      do j=1, na_cols
Andreas Marek's avatar
Andreas Marek committed
222
223
224
225
226
227
228
229
        do i=1,na_rows
#ifdef TEST_DOUBLE
          a_complex(i,j) = dcmplx(0.0, a_skewsymmetric(i,j))
#endif
#ifdef TEST_SINGLE
          a_complex(i,j) = cmplx(0.0, a_skewsymmetric(i,j))
#endif
        enddo
Carolin Penke's avatar
Carolin Penke committed
230
231
      enddo
   
Andreas Marek's avatar
Andreas Marek committed
232
233


Carolin Penke's avatar
Carolin Penke committed
234
235
   z_complex(1:na_rows,1:na_cols)  = a_complex(1:na_rows,1:na_cols)
   as_complex(1:na_rows,1:na_cols) = a_complex(1:na_rows,1:na_cols)
Andreas Marek's avatar
Andreas Marek committed
236

237
   ! first set up and solve the brute force problem
238
   e_complex => elpa_allocate(error_elpa)
239
   call set_basic_params(e_complex, na, nev, na_rows, na_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
240

241
   call e_complex%set("timings",1, error_elpa)
Andreas Marek's avatar
Andreas Marek committed
242

243
   call e_complex%set("debug",1,error_elpa)
244
245
246
247
248
249
250
#if TEST_NVIDIA_GPU == 1 || (TEST_NVIDIA_GPU == 0) && (TEST_AMD_GPU == 0)   
   call e_complex%set("nvidia-gpu", TEST_GPU,error_elpa)
#endif
#if TEST_AMD_GPU == 1
   call e_complex%set("amd-gpu", TEST_GPU,error_elpa)
#endif

251
   call e_complex%set("omp_threads", 8, error_elpa)
Andreas Marek's avatar
Andreas Marek committed
252

253
   assert_elpa_ok(e_complex%setup())
254
   call e_complex%set("solver", elpa_solver_2stage, error_elpa)
Andreas Marek's avatar
Andreas Marek committed
255

256
   call e_complex%timer_start("eigenvectors: brute force as complex matrix")
257
   call e_complex%eigenvectors(a_complex, ev_complex, z_complex, error_elpa)
258
   call e_complex%timer_stop("eigenvectors: brute force as complex matrix")
Andreas Marek's avatar
Andreas Marek committed
259

260
261
   if (myid .eq. 0) then
     print *, ""
262
     call e_complex%print_times("eigenvectors: brute force as complex matrix")
263
   endif 
Carolin Penke's avatar
Carolin Penke committed
264
#ifdef WITH_MPI
265
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
Carolin Penke's avatar
Carolin Penke committed
266
267
#endif     
!      as_complex(:,:) = z_complex(:,:)
268
269
270
271
272
#ifdef TEST_SINGLE
     status = check_correctness_evp_numeric_residuals_complex_single(na, nev, as_complex, z_complex, ev_complex, sc_desc, &
                                                    nblk, myid, np_rows,np_cols, my_prow, my_pcol)
#else
     status = check_correctness_evp_numeric_residuals_complex_double(na, nev, as_complex, z_complex, ev_complex, sc_desc, &
273
                                                    nblk, myid, np_rows,np_cols, my_prow, my_pcol)
274
#endif
Andreas Marek's avatar
Andreas Marek committed
275
276
    status = 0
    call check_status(status, myid)
Andreas Marek's avatar
Andreas Marek committed
277
278

#ifdef WITH_MPI
279
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
Andreas Marek's avatar
Andreas Marek committed
280
#endif
281
   ! now run the skewsymmetric case
282
   e_skewsymmetric => elpa_allocate(error_elpa)
283
   call set_basic_params(e_skewsymmetric, na, nev, na_rows, na_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
284

285
   call e_skewsymmetric%set("timings",1, error_elpa)
286

287
   call e_skewsymmetric%set("debug",1,error_elpa)
288
289
290
291
292
293
294

#if TEST_NVIDIA_GPU == 1 || (TEST_NVIDIA_GPU == 0) && (TEST_AMD_GPU == 0)
   call e_skewsymmetric%set("nvidia-gpu", TEST_GPU,error_elpa)
#endif
#if TEST_AMD_GPU == 1
   call e_skewsymmetric%set("amd-gpu", TEST_GPU,error_elpa)
#endif
295
   call e_skewsymmetric%set("omp_threads",8, error_elpa)
296
297

   assert_elpa_ok(e_skewsymmetric%setup())
Carolin Penke's avatar
Carolin Penke committed
298
   
299
   call e_skewsymmetric%set("solver", elpa_solver_2stage, error_elpa)
300
301

   call e_skewsymmetric%timer_start("eigenvectors: skewsymmetric ")
Andreas Marek's avatar
Andreas Marek committed
302
   call e_skewsymmetric%skew_eigenvectors(a_skewsymmetric, ev_skewsymmetric, z_skewsymmetric, error_elpa)
303
   call e_skewsymmetric%timer_stop("eigenvectors: skewsymmetric ")
Andreas Marek's avatar
Andreas Marek committed
304
305
306

   if (myid .eq. 0) then
     print *, ""
Carolin Penke's avatar
Carolin Penke committed
307
     call e_skewsymmetric%print_times("eigenvectors: skewsymmetric")
Carolin Penke's avatar
Carolin Penke committed
308
   endif
Carolin Penke's avatar
Carolin Penke committed
309
   
Carolin Penke's avatar
Carolin Penke committed
310
311
   ! check eigenvalues
   do i=1, na
Carolin Penke's avatar
Carolin Penke committed
312
     if (myid == 0) then
Andreas Marek's avatar
Andreas Marek committed
313
#ifdef TEST_DOUBLE
314
       if (abs(ev_complex(i)-ev_skewsymmetric(i))/abs(ev_complex(i)) .gt. 1e-10) then
Andreas Marek's avatar
Andreas Marek committed
315
316
317
318
#endif
#ifdef TEST_SINGLE
       if (abs(ev_complex(i)-ev_skewsymmetric(i))/abs(ev_complex(i)) .gt. 1e-4) then
#endif
Carolin Penke's avatar
Carolin Penke committed
319
         print *,"ev: i=",i,ev_complex(i),ev_skewsymmetric(i)
Carolin Penke's avatar
Carolin Penke committed
320
321
         status = 1
     endif
Carolin Penke's avatar
Carolin Penke committed
322
323
     endif
   enddo
Andreas Marek's avatar
Andreas Marek committed
324
325


326
!    call check_status(status, myid)
Carolin Penke's avatar
Carolin Penke committed
327
   
Carolin Penke's avatar
Carolin Penke committed
328
329
330
   z_complex(:,:) = 0
   do j=1, na_cols
     do i=1,na_rows
Andreas Marek's avatar
Andreas Marek committed
331
#ifdef TEST_DOUBLE
332
       z_complex(i,j) = dcmplx(z_skewsymmetric(i,j), z_skewsymmetric(i,na_cols+j))
Andreas Marek's avatar
Andreas Marek committed
333
334
335
336
#endif
#ifdef TEST_SINGLE
       z_complex(i,j) = cmplx(z_skewsymmetric(i,j), z_skewsymmetric(i,na_cols+j))
#endif
Carolin Penke's avatar
Carolin Penke committed
337
338
     enddo
   enddo
339
340
341
#ifdef WITH_MPI
   call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
#endif
Carolin Penke's avatar
Carolin Penke committed
342

343
344
345
346
347
348
349
#ifdef TEST_SINGLE
   status = check_correctness_evp_numeric_residuals_ss_real_single(na, nev, as_skewsymmetric, z_complex, ev_skewsymmetric, &
                              sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
#else
   status = check_correctness_evp_numeric_residuals_ss_real_double(na, nev, as_skewsymmetric, z_complex, ev_skewsymmetric, &
                              sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)
#endif
Carolin Penke's avatar
Carolin Penke committed
350
   
351
#ifdef WITH_MPI
352
    call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
353
#endif
354
355
   call elpa_deallocate(e_complex,error_elpa)
   call elpa_deallocate(e_skewsymmetric,error_elpa)
356

Andreas Marek's avatar
Andreas Marek committed
357

358
359
360
   !to do 
   ! - check whether brute-force check_correctness_evp_numeric_residuals worsk (complex ev)
   ! - invent a test for skewsymmetric residuals
Andreas Marek's avatar
Andreas Marek committed
361

362
363
364
365
   deallocate(a_complex)
   deallocate(as_complex)
   deallocate(z_complex)
   deallocate(ev_complex)
Andreas Marek's avatar
Andreas Marek committed
366

367
368
369
370
   deallocate(a_skewsymmetric)
   deallocate(as_skewsymmetric)
   deallocate(z_skewsymmetric)
   deallocate(ev_skewsymmetric)
371
   call elpa_uninit(error_elpa)
Andreas Marek's avatar
Andreas Marek committed
372

373
374


Andreas Marek's avatar
Andreas Marek committed
375
376
377
378
379
380
381
382
383
384
385
#ifdef WITH_MPI
   call blacs_gridexit(my_blacs_ctxt)
   call mpi_finalize(mpierr)
#endif

   call exit(status)

contains
   subroutine set_basic_params(elpa, na, nev, na_rows, na_cols, my_prow, my_pcol)
     implicit none
     class(elpa_t), pointer      :: elpa
386
387
388
389
390
391
392
393
394
395
396
397
     TEST_INT_TYPE, intent(in)         :: na, nev, na_rows, na_cols, my_prow, my_pcol

     call elpa%set("na", int(na,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("nev", int(nev,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("local_nrows", int(na_rows,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("local_ncols", int(na_cols,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("nblk", int(nblk,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
Andreas Marek's avatar
Andreas Marek committed
398
399

#ifdef WITH_MPI
400
401
402
403
404
405
     call elpa%set("mpi_comm_parent", int(MPI_COMM_WORLD,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("process_row", int(my_prow,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
     call elpa%set("process_col", int(my_pcol,kind=c_int), error_elpa)
     assert_elpa_ok(error_elpa)
Andreas Marek's avatar
Andreas Marek committed
406
407
#endif
   end subroutine
408
409
   subroutine check_status(status, myid)
     implicit none
410
411
     TEST_INT_TYPE, intent(in) :: status, myid
     TEST_INT_MPI_TYPE         :: mpierr
412
413
414
415
416
417
418
419
     if (status /= 0) then
       if (myid == 0) print *, "Result incorrect!"
#ifdef WITH_MPI
       call mpi_finalize(mpierr)
#endif
       call exit(status)
     endif
   end subroutine
Andreas Marek's avatar
Andreas Marek committed
420
end program