Line data Source code
1 : ! This file is part of ELPA.
2 : !
3 : ! The ELPA library was originally created by the ELPA consortium,
4 : ! consisting of the following organizations:
5 : !
6 : ! - Max Planck Computing and Data Facility (MPCDF), formerly known as
7 : ! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
8 : ! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
9 : ! Informatik,
10 : ! - Technische Universität München, Lehrstuhl für Informatik mit
11 : ! Schwerpunkt Wissenschaftliches Rechnen ,
12 : ! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
13 : ! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
14 : ! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
15 : ! and
16 : ! - IBM Deutschland GmbH
17 : !
18 : !
19 : ! More information can be found here:
20 : ! http://elpa.mpcdf.mpg.de/
21 : !
22 : ! ELPA is free software: you can redistribute it and/or modify
23 : ! it under the terms of the version 3 of the license of the
24 : ! GNU Lesser General Public License as published by the Free
25 : ! Software Foundation.
26 : !
27 : ! ELPA is distributed in the hope that it will be useful,
28 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : ! GNU Lesser General Public License for more details.
31 : !
32 : ! You should have received a copy of the GNU Lesser General Public License
33 : ! along with ELPA. If not, see <http://www.gnu.org/licenses/>
34 : !
35 : ! ELPA reflects a substantial effort on the part of the original
36 : ! ELPA consortium, and we ask you to respect the spirit of the
37 : ! license that we chose: i.e., please contribute any changes you
38 : ! may have back to the original ELPA library distribution, and keep
39 : ! any derivatives of ELPA under the same license that we chose for
40 : ! the original distribution, the GNU Lesser General Public License.
41 : !
42 : !
43 : #include "config-f90.h"
44 : !>
45 : !> Fortran test programm to demonstrates the use of
46 : !> the elpa_solve_tridi library function
47 : !>
48 192 : program test_solve_tridi
49 :
50 : !
51 : ! Copyright of the original code rests with the authors inside the ELPA
52 : ! consortium. The copyright of any additional modifications shall rest
53 : ! with their original authors, but shall adhere to the licensing terms
54 : ! distributed along with the original code in the file "COPYING".
55 : !
56 : !-------------------------------------------------------------------------------
57 192 : use elpa1
58 : use test_util
59 :
60 : use test_read_input_parameters
61 : use test_check_correctness
62 : use test_setup_mpi
63 : use test_blacs_infrastructure
64 : use test_prepare_matrix
65 :
66 : #ifdef HAVE_REDIRECT
67 : use test_redirect
68 : #endif
69 : use test_output_type
70 :
71 : implicit none
72 :
73 : !-------------------------------------------------------------------------------
74 : ! Please set system size parameters below!
75 : ! na: System size
76 : ! nev: Number of eigenvectors to be calculated
77 : ! nblk: Blocking factor in block cyclic distribution
78 : !-------------------------------------------------------------------------------
79 : integer(kind=ik) :: nblk
80 : integer(kind=ik) :: na, nev
81 :
82 : integer(kind=ik) :: np_rows, np_cols, na_rows, na_cols
83 :
84 : integer(kind=ik) :: myid, nprocs, my_prow, my_pcol, mpi_comm_rows, mpi_comm_cols
85 : integer(kind=ik) :: i, mpierr, my_blacs_ctxt, sc_desc(9), info, nprow, npcol
86 :
87 : integer, external :: numroc
88 :
89 384 : real(kind=rk8), allocatable :: a(:,:), d(:), e(:), ev_analytic(:), ev(:)
90 : real(kind=rk8) :: diagonalELement, subdiagonalElement
91 :
92 192 : real(kind=rk8), allocatable :: as(:,:)
93 : real(kind=rk8) :: tmp
94 : integer(kind=ik) :: loctmp ,rowLocal, colLocal
95 :
96 :
97 : real(kind=rk8) :: maxerr
98 :
99 : logical :: wantDebug
100 :
101 : real(kind=rk8), parameter :: pi = 3.141592653589793238462643383279_rk8
102 :
103 : integer(kind=ik) :: STATUS
104 : #ifdef WITH_OPENMP
105 : integer(kind=ik) :: omp_get_max_threads, required_mpi_thread_level, &
106 : provided_mpi_thread_level
107 : #endif
108 : type(output_t) :: write_to_file
109 : logical :: success
110 : character(len=8) :: task_suffix
111 : integer(kind=ik) :: j
112 : !-------------------------------------------------------------------------------
113 :
114 192 : success = .true.
115 :
116 192 : call read_input_parameters(na, nev, nblk, write_to_file)
117 :
118 : !-------------------------------------------------------------------------------
119 : ! MPI Initialization
120 192 : call setup_mpi(myid, nprocs)
121 :
122 192 : STATUS = 0
123 :
124 : !#define DATATYPE REAL
125 : !#define ELPA1
126 : !#include "../../elpa_print_headers.F90"
127 :
128 192 : do np_cols = NINT(SQRT(REAL(nprocs))),2,-1
129 0 : if(mod(nprocs,np_cols) == 0 ) exit
130 : enddo
131 :
132 : ! at the end of the above loop, nprocs is always divisible by np_cols
133 :
134 192 : np_rows = nprocs/np_cols
135 :
136 192 : if(myid==0) then
137 128 : print *
138 128 : print '(a)','Test program for elpa_solve_tridi with a Toeplitz matrix'
139 128 : print *
140 128 : print '(3(a,i0))','Matrix size=',na,', Number of eigenvectors=',nev,', Block size=',nblk
141 128 : print '(3(a,i0))','Number of processor rows=',np_rows,', cols=',np_cols,', total=',nprocs
142 128 : print *
143 : endif
144 :
145 : !-------------------------------------------------------------------------------
146 : ! Set up BLACS context and MPI communicators
147 : !
148 : ! The BLACS context is only necessary for using Scalapack.
149 : !
150 : ! For ELPA, the MPI communicators along rows/cols are sufficient,
151 : ! and the grid setup may be done in an arbitrary way as long as it is
152 : ! consistent (i.e. 0<=my_prow<np_rows, 0<=my_pcol<np_cols and every
153 : ! process has a unique (my_prow,my_pcol) pair).
154 :
155 : call set_up_blacsgrid(mpi_comm_world, np_rows, np_cols, 'C', &
156 192 : my_blacs_ctxt, my_prow, my_pcol)
157 :
158 192 : if (myid==0) then
159 128 : print '(a)','| Past BLACS_Gridinfo.'
160 : end if
161 :
162 : ! All ELPA routines need MPI communicators for communicating within
163 : ! rows or columns of processes, these are set in elpa_get_communicators.
164 :
165 : mpierr = elpa_get_communicators(mpi_comm_world, my_prow, my_pcol, &
166 192 : mpi_comm_rows, mpi_comm_cols)
167 :
168 192 : if (myid==0) then
169 128 : print '(a)','| Past split communicator setup for rows and columns.'
170 : end if
171 :
172 : call set_up_blacs_descriptor(na ,nblk, my_prow, my_pcol, np_rows, np_cols, &
173 192 : na_rows, na_cols, sc_desc, my_blacs_ctxt, info)
174 :
175 192 : if (myid==0) then
176 128 : print '(a)','| Past scalapack descriptor setup.'
177 : end if
178 :
179 : !-------------------------------------------------------------------------------
180 : ! Allocate matrices and set up a test toeplitz matrix for solve_tridi
181 :
182 192 : allocate(a (na_rows,na_cols))
183 192 : allocate(as(na_rows,na_cols))
184 :
185 192 : allocate(d (na))
186 192 : allocate(e (na))
187 192 : allocate(ev_analytic(na))
188 192 : allocate(ev(na))
189 :
190 192 : a(:,:) = 0.0_rk8
191 :
192 :
193 : ! changeable numbers here would be nice
194 192 : diagonalElement = 0.45_rk8
195 192 : subdiagonalElement = 0.78_rk8
196 :
197 192 : d(:) = diagonalElement
198 192 : e(:) = subdiagonalElement
199 :
200 :
201 : ! set up the diagonal and subdiagonals (for general solver test)
202 28992 : do i=1, na ! for diagonal elements
203 28800 : if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
204 19200 : a(rowLocal,colLocal) = diagonalElement
205 : endif
206 : enddo
207 :
208 28800 : do i=1, na-1
209 28608 : if (map_global_array_index_to_local_index(i, i+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
210 19072 : a(rowLocal,colLocal) = subdiagonalElement
211 : endif
212 : enddo
213 :
214 28800 : do i=2, na
215 28608 : if (map_global_array_index_to_local_index(i, i-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
216 19072 : a(rowLocal,colLocal) = subdiagonalElement
217 : endif
218 : enddo
219 :
220 192 : as = a
221 :
222 :
223 :
224 : !-------------------------------------------------------------------------------
225 : ! Calculate eigenvalues/eigenvectors
226 :
227 192 : if (myid==0) then
228 128 : print '(a)','| Entering elpa_solve_tridi ... '
229 128 : print *
230 : end if
231 :
232 : #ifdef WITH_MPI
233 128 : call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
234 : #endif
235 :
236 : success = elpa_solve_tridi_double(na, nev, d, e, a, na_rows, nblk, na_cols, mpi_comm_rows, &
237 192 : mpi_comm_cols, wantDebug)
238 :
239 192 : if (.not.(success)) then
240 0 : write(error_unit,*) "elpa_solve_tridi produced an error! Aborting..."
241 : #ifdef WITH_MPI
242 0 : call MPI_ABORT(mpi_comm_world, 1, mpierr)
243 : #else
244 0 : call exit(1)
245 : #endif
246 : endif
247 :
248 192 : if (myid==0) then
249 128 : print '(a)','| elpa_solve_tridi complete.'
250 128 : print *
251 : end if
252 :
253 :
254 192 : ev = d
255 :
256 : !
257 : ! if(myid == 0) print *,'Time tridiag_real :',time_evp_fwd
258 : ! if(myid == 0) print *,'Time solve_tridi :',time_evp_solve
259 : ! if(myid == 0) print *,'Time trans_ev_real :',time_evp_back
260 : ! if(myid == 0) print *,'Total time (sum above):',time_evp_back+time_evp_solve+time_evp_fwd
261 : !
262 : ! if(write_to_file%eigenvectors) then
263 : ! write(unit = task_suffix, fmt = '(i8.8)') myid
264 : ! open(17,file="EVs_real_out_task_"//task_suffix(1:8)//".txt",form='formatted',status='new')
265 : ! write(17,*) "Part of eigenvectors: na_rows=",na_rows,"of na=",na," na_cols=",na_cols," of na=",na
266 : !
267 : ! do i=1,na_rows
268 : ! do j=1,na_cols
269 : ! write(17,*) "row=",i," col=",j," element of eigenvector=",z(i,j)
270 : ! enddo
271 : ! enddo
272 : ! close(17)
273 : ! endif
274 : !
275 : ! if(write_to_file%eigenvalues) then
276 : ! if (myid == 0) then
277 : ! open(17,file="Eigenvalues_real_out.txt",form='formatted',status='new')
278 : ! do i=1,na
279 : ! write(17,*) i,ev(i)
280 : ! enddo
281 : ! close(17)
282 : ! endif
283 : ! endif
284 : !
285 : !
286 : !-------------------------------------------------------------------------------
287 :
288 : ! analytic solution
289 28992 : do i=1, na
290 28800 : ev_analytic(i) = diagonalElement + 2.0_rk8 * subdiagonalElement *cos( pi*real(i,kind=rk8)/ real(na+1,kind=rk8) )
291 : enddo
292 :
293 : ! sort analytic solution:
294 :
295 : ! this hack is neihter elegant, nor optimized: for huge matrixes it might be expensive
296 : ! a proper sorting algorithmus might be implemented here
297 :
298 192 : tmp = minval(ev_analytic)
299 192 : loctmp = minloc(ev_analytic, 1)
300 :
301 192 : ev_analytic(loctmp) = ev_analytic(1)
302 192 : ev_analytic(1) = tmp
303 :
304 28800 : do i=2, na
305 28608 : tmp = ev_analytic(i)
306 2174208 : do j= i, na
307 2145600 : if (ev_analytic(j) .lt. tmp) then
308 1051392 : tmp = ev_analytic(j)
309 1051392 : loctmp = j
310 : endif
311 : enddo
312 28608 : ev_analytic(loctmp) = ev_analytic(i)
313 28608 : ev_analytic(i) = tmp
314 : enddo
315 :
316 : ! compute a simple error max of eigenvalues
317 192 : maxerr = 0.0_rk8
318 192 : maxerr = maxval( (d(:) - ev_analytic(:))/ev_analytic(:) , 1)
319 :
320 192 : if (maxerr .gt. 8.e-13) then
321 0 : if (myid .eq. 0) then
322 0 : print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
323 : endif
324 :
325 0 : status = 1
326 :
327 : endif
328 :
329 : ! Test correctness of result (using plain scalapack routines)
330 192 : status = check_correctness_evp_numeric_residuals(na, nev, as, a, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
331 :
332 192 : deallocate(a)
333 :
334 192 : deallocate(as)
335 192 : deallocate(d)
336 :
337 192 : deallocate(e)
338 192 : deallocate(ev)
339 192 : deallocate(ev_analytic)
340 :
341 : #ifdef WITH_MPI
342 128 : call blacs_gridexit(my_blacs_ctxt)
343 128 : call mpi_finalize(mpierr)
344 : #endif
345 :
346 192 : call EXIT(STATUS)
347 :
348 :
349 : end
350 :
351 : !-------------------------------------------------------------------------------
|
