test_skewsymmetric.F90 12.8 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
48
49
50
51
52
53
54
55
56
57
58
59
60
!    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]
! 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
61
62
#  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
63
64
65
66
67
68
#  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
69
70
#  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
71
72
73
74
75
76
77
78
79
80
81
82
83
#  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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#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
Andreas Marek's avatar
Andreas Marek committed
98
99
100
101
102
#include "assert.h"

program test
   use elpa

103
   !use test_util
Andreas Marek's avatar
Andreas Marek committed
104
105
106
107
108
   use test_setup_mpi
   use test_prepare_matrix
   use test_read_input_parameters
   use test_blacs_infrastructure
   use test_check_correctness
109
   use precision_for_tests
Andreas Marek's avatar
Andreas Marek committed
110
111
112
113
114
115
116
117
   use iso_fortran_env

#ifdef HAVE_REDIRECT
   use test_redirect
#endif
   implicit none

   ! matrix dimensions
118
   TEST_INT_TYPE                          :: na, nev, nblk
Andreas Marek's avatar
Andreas Marek committed
119
120

   ! mpi
121
122
123
124
125
   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
126
127

   ! blacs
128
   character(len=1)                 :: layout
129
   TEST_INT_TYPE                          :: my_blacs_ctxt, sc_desc(9), info, nprow, npcol
Andreas Marek's avatar
Andreas Marek committed
130
131

   ! The Matrix
132
133
   MATRIX_TYPE, allocatable         :: a_skewsymmetric(:,:), as_skewsymmetric(:,:)
   MATRIX_TYPE_COMPLEX, allocatable :: a_complex(:,:), as_complex(:,:)
Andreas Marek's avatar
Andreas Marek committed
134
   ! eigenvectors
135
136
   MATRIX_TYPE, allocatable         :: z_skewsymmetric(:,:)
   MATRIX_TYPE_COMPLEX, allocatable :: z_complex(:,:)
Andreas Marek's avatar
Andreas Marek committed
137
   ! eigenvalues
138
   EV_TYPE, allocatable             :: ev_skewsymmetric(:), ev_complex(:)
Andreas Marek's avatar
Andreas Marek committed
139

140
141
   TEST_INT_TYPE                    :: status, i, j
   integer(kind=c_int)              :: error_elpa
Andreas Marek's avatar
Andreas Marek committed
142

143
144
   type(output_t)                   :: write_to_file
   class(elpa_t), pointer           :: e_complex, e_skewsymmetric
145
           
Andreas Marek's avatar
Andreas Marek committed
146
147
148
149
150
151
152
153
154
155
156
157
158
   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
159
! 
Andreas Marek's avatar
Andreas Marek committed
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
   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

   call set_up_blacsgrid(mpi_comm_world, np_rows, np_cols, layout, &
                         my_blacs_ctxt, my_prow, my_pcol)

   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)

185
186
   allocate(a_skewsymmetric (na_rows,na_cols))
   allocate(as_skewsymmetric(na_rows,na_cols))
Carolin Penke's avatar
Carolin Penke committed
187
   allocate(z_skewsymmetric (na_rows,2*na_cols))
188
   allocate(ev_skewsymmetric(na))
Andreas Marek's avatar
Andreas Marek committed
189

190
191
192
   a_skewsymmetric(:,:) = 0.0
   z_skewsymmetric(:,:) = 0.0
   ev_skewsymmetric(:) = 0.0
Andreas Marek's avatar
Andreas Marek committed
193

Carolin Penke's avatar
Carolin Penke committed
194
195
   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
196
   
197
   !call MPI_BARRIER(MPI_COMM_WORLD, mpierr)  
198
   as_skewsymmetric(:,:) = a_skewsymmetric(:,:)
Carolin Penke's avatar
Carolin Penke committed
199
   
Andreas Marek's avatar
Andreas Marek committed
200

201
202
203
204
205
   ! 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
206

Carolin Penke's avatar
Carolin Penke committed
207
208
209
   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
210
   
Andreas Marek's avatar
Andreas Marek committed
211

Carolin Penke's avatar
Carolin Penke committed
212
      do j=1, na_cols
Andreas Marek's avatar
Andreas Marek committed
213
214
215
216
217
218
219
220
        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
221
222
      enddo
   
Andreas Marek's avatar
Andreas Marek committed
223
224


Carolin Penke's avatar
Carolin Penke committed
225
226
   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
227

228
229
230
   ! first set up and solve the brute force problem
   e_complex => elpa_allocate()
   call set_basic_params(e_complex, na, nev, na_rows, na_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
231

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

234
   call e_complex%set("debug",1)
Carolin Penke's avatar
Carolin Penke committed
235
   call e_complex%set("gpu", 0)
236
   call e_complex%set("omp_threads", 8, error_elpa)
Andreas Marek's avatar
Andreas Marek committed
237

238
   assert_elpa_ok(e_complex%setup())
239
   call e_complex%set("solver", elpa_solver_2stage, error_elpa)
Andreas Marek's avatar
Andreas Marek committed
240

241
   call e_complex%timer_start("eigenvectors: brute force as complex matrix")
242
   call e_complex%eigenvectors(a_complex, ev_complex, z_complex, error_elpa)
243
   call e_complex%timer_stop("eigenvectors: brute force as complex matrix")
Andreas Marek's avatar
Andreas Marek committed
244

245
246
   if (myid .eq. 0) then
     print *, ""
247
     call e_complex%print_times("eigenvectors: brute force as complex matrix")
248
   endif 
Carolin Penke's avatar
Carolin Penke committed
249
#ifdef WITH_MPI
250
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
Carolin Penke's avatar
Carolin Penke committed
251
252
#endif     
!      as_complex(:,:) = z_complex(:,:)
253
254
255
256
257
#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, &
258
                                                    nblk, myid, np_rows,np_cols, my_prow, my_pcol)
259
#endif
Andreas Marek's avatar
Andreas Marek committed
260
261
    status = 0
    call check_status(status, myid)
Andreas Marek's avatar
Andreas Marek committed
262
263

#ifdef WITH_MPI
264
     call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
Andreas Marek's avatar
Andreas Marek committed
265
#endif
266
267
268
   ! now run the skewsymmetric case
   e_skewsymmetric => elpa_allocate()
   call set_basic_params(e_skewsymmetric, na, nev, na_rows, na_cols, my_prow, my_pcol)
Andreas Marek's avatar
Andreas Marek committed
269

270
   call e_skewsymmetric%set("timings",1, error_elpa)
271
272

   call e_skewsymmetric%set("debug",1)
Carolin Penke's avatar
Carolin Penke committed
273
   call e_skewsymmetric%set("gpu", 0)
274
   call e_skewsymmetric%set("omp_threads",8, error_elpa)
275
276
277

   call e_skewsymmetric%set("is_skewsymmetric",1)
   assert_elpa_ok(e_skewsymmetric%setup())
Carolin Penke's avatar
Carolin Penke committed
278
   
279
   call e_skewsymmetric%set("solver", elpa_solver_2stage, error_elpa)
280

281
   call e_skewsymmetric%get("is_skewsymmetric", i,error_elpa)
Carolin Penke's avatar
Carolin Penke committed
282
   
283
   call e_skewsymmetric%timer_start("eigenvectors: skewsymmetric ")
284
   call e_skewsymmetric%eigenvectors(a_skewsymmetric, ev_skewsymmetric, z_skewsymmetric, error_elpa)
285
   call e_skewsymmetric%timer_stop("eigenvectors: skewsymmetric ")
Andreas Marek's avatar
Andreas Marek committed
286
287
288

   if (myid .eq. 0) then
     print *, ""
Carolin Penke's avatar
Carolin Penke committed
289
     call e_skewsymmetric%print_times("eigenvectors: skewsymmetric")
Carolin Penke's avatar
Carolin Penke committed
290
   endif
Carolin Penke's avatar
Carolin Penke committed
291
   
Carolin Penke's avatar
Carolin Penke committed
292
293
   ! check eigenvalues
   do i=1, na
Carolin Penke's avatar
Carolin Penke committed
294
     if (myid == 0) then
Andreas Marek's avatar
Andreas Marek committed
295
#ifdef TEST_DOUBLE
296
       if (abs(ev_complex(i)-ev_skewsymmetric(i))/abs(ev_complex(i)) .gt. 1e-10) then
Andreas Marek's avatar
Andreas Marek committed
297
298
299
300
#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
301
         print *,"ev: i=",i,ev_complex(i),ev_skewsymmetric(i)
Carolin Penke's avatar
Carolin Penke committed
302
303
         status = 1
     endif
Carolin Penke's avatar
Carolin Penke committed
304
305
     endif
   enddo
Andreas Marek's avatar
Andreas Marek committed
306
307


308
!    call check_status(status, myid)
Carolin Penke's avatar
Carolin Penke committed
309
   
Carolin Penke's avatar
Carolin Penke committed
310
311
312
   z_complex(:,:) = 0
   do j=1, na_cols
     do i=1,na_rows
Andreas Marek's avatar
Andreas Marek committed
313
#ifdef TEST_DOUBLE
314
       z_complex(i,j) = dcmplx(z_skewsymmetric(i,j), z_skewsymmetric(i,na_cols+j))
Andreas Marek's avatar
Andreas Marek committed
315
316
317
318
#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
319
320
     enddo
   enddo
321
322
323
#ifdef WITH_MPI
   call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
#endif
Carolin Penke's avatar
Carolin Penke committed
324
325
326
   status = check_correctness_evp_numeric_residuals_ss(na, nev, as_skewsymmetric, z_complex, ev_skewsymmetric, &
                              sc_desc, nblk, myid, np_rows,np_cols, my_prow, my_pcol)

Carolin Penke's avatar
Carolin Penke committed
327
   
328
#ifdef WITH_MPI
329
    call MPI_BARRIER(MPI_COMM_WORLD, mpierr)
330
331
332
333
#endif
   call elpa_deallocate(e_complex)
   call elpa_deallocate(e_skewsymmetric)

Andreas Marek's avatar
Andreas Marek committed
334

335
336
337
   !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
338

339
340
341
342
   deallocate(a_complex)
   deallocate(as_complex)
   deallocate(z_complex)
   deallocate(ev_complex)
Andreas Marek's avatar
Andreas Marek committed
343

344
345
346
347
   deallocate(a_skewsymmetric)
   deallocate(as_skewsymmetric)
   deallocate(z_skewsymmetric)
   deallocate(ev_skewsymmetric)
Andreas Marek's avatar
Andreas Marek committed
348
349
   call elpa_uninit()

350
351


Andreas Marek's avatar
Andreas Marek committed
352
353
354
355
356
357
358
359
360
361
362
#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
363
364
365
366
367
368
369
370
371
372
373
374
     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
375
376

#ifdef WITH_MPI
377
378
379
380
381
382
     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
383
384
#endif
   end subroutine
385
386
   subroutine check_status(status, myid)
     implicit none
387
388
     TEST_INT_TYPE, intent(in) :: status, myid
     TEST_INT_MPI_TYPE         :: mpierr
389
390
391
392
393
394
395
396
     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
397
end program