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 : ! Author: A. Marek, MPCDF
43 :
44 31536 : function check_correctness_evp_numeric_residuals_&
45 : &MATH_DATATYPE&
46 : &_&
47 : &PRECISION&
48 31536 : & (na, nev, as, z, ev, sc_desc, nblk, myid, np_rows, np_cols, my_prow, my_pcol, bs) result(status)
49 : implicit none
50 : #include "../../src/general/precision_kinds.F90"
51 : integer(kind=ik) :: status
52 : integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, np_cols, my_prow, my_pcol
53 : MATH_DATATYPE(kind=rck), intent(in) :: as(:,:), z(:,:)
54 : MATH_DATATYPE(kind=rck), intent(in), optional :: bs(:,:)
55 : real(kind=rk) :: ev(:)
56 63072 : MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
57 : MATH_DATATYPE(kind=rck) :: xc
58 :
59 : #ifndef WITH_MPI
60 : #if REALCASE == 1
61 : real(kind=rck) :: dnrm2, snrm2
62 : #endif
63 : #if COMPLEXCASE == 1
64 : complex(kind=rck) :: zdotc, cdotc
65 : #endif /* COMPLEXCASE */
66 : #endif
67 :
68 : integer(kind=ik) :: sc_desc(:)
69 :
70 : integer(kind=ik) :: i, rowLocal, colLocal
71 : real(kind=rck) :: err, errmax
72 :
73 : integer :: mpierr
74 :
75 : ! tolerance for the residual test for different math type/precision setups
76 : real(kind=rk), parameter :: tol_res_real_double = 5e-12_rk
77 : real(kind=rk), parameter :: tol_res_real_single = 3e-2_rk
78 : real(kind=rk), parameter :: tol_res_complex_double = 5e-12_rk
79 : real(kind=rk), parameter :: tol_res_complex_single = 3e-2_rk
80 : real(kind=rk) :: tol_res = tol_res_&
81 : &MATH_DATATYPE&
82 : &_&
83 : &PRECISION
84 : ! precision of generalized problem is lower
85 : real(kind=rk), parameter :: generalized_penalty = 10.0_rk
86 :
87 : ! tolerance for the orthogonality test for different math type/precision setups
88 : real(kind=rk), parameter :: tol_orth_real_double = 5e-11_rk
89 : real(kind=rk), parameter :: tol_orth_real_single = 9e-2_rk
90 : real(kind=rk), parameter :: tol_orth_complex_double = 5e-11_rk
91 : real(kind=rk), parameter :: tol_orth_complex_single = 9e-3_rk
92 : real(kind=rk), parameter :: tol_orth = tol_orth_&
93 : &MATH_DATATYPE&
94 : &_&
95 : &PRECISION
96 :
97 31536 : if(present(bs)) then
98 0 : tol_res = generalized_penalty * tol_res
99 : endif
100 31536 : status = 0
101 :
102 : ! 1. Residual (maximum of || A*Zi - Zi*EVi ||)
103 :
104 : ! tmp1 = Zi*EVi
105 31536 : tmp1(:,:) = z(:,:)
106 5232336 : do i=1,nev
107 5200800 : xc = ev(i)
108 : #ifdef WITH_MPI
109 : call p&
110 : &BLAS_CHAR&
111 3486400 : &scal(na, xc, tmp1, 1, i, sc_desc, 1)
112 : #else /* WITH_MPI */
113 : call BLAS_CHAR&
114 1714400 : &scal(na,xc,tmp1(:,i),1)
115 : #endif /* WITH_MPI */
116 : enddo
117 :
118 : ! for generalized EV problem, multiply by bs as well
119 : ! tmp2 = B * tmp1
120 31536 : if(present(bs)) then
121 : #ifdef WITH_MPI
122 : call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, bs, 1, 1, sc_desc, &
123 0 : tmp1, 1, 1, sc_desc, ZERO, tmp2, 1, 1, sc_desc)
124 : #else /* WITH_MPI */
125 0 : call PRECISION_GEMM('N','N',na,nev,na,ONE,bs,na,tmp1,na,ZERO,tmp2,na)
126 : #endif /* WITH_MPI */
127 : else
128 : ! normal eigenvalue problem .. no need to multiply
129 31536 : tmp2(:,:) = tmp1(:,:)
130 : end if
131 :
132 : ! tmp1 = A * Z
133 : ! as is original stored matrix, Z are the EVs
134 : #ifdef WITH_MPI
135 : call scal_PRECISION_GEMM('N', 'N', na, nev, na, ONE, as, 1, 1, sc_desc, &
136 21152 : z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
137 : #else /* WITH_MPI */
138 10384 : call PRECISION_GEMM('N','N',na,nev,na,ONE,as,na,z,na,ZERO,tmp1,na)
139 : #endif /* WITH_MPI */
140 :
141 : ! tmp1 = A*Zi - Zi*EVi
142 31536 : tmp1(:,:) = tmp1(:,:) - tmp2(:,:)
143 :
144 : ! Get maximum norm of columns of tmp1
145 31536 : errmax = 0.0_rk
146 :
147 5232336 : do i=1,nev
148 : #if REALCASE == 1
149 2863200 : err = 0.0_rk
150 : #ifdef WITH_MPI
151 1918400 : call scal_PRECISION_NRM2(na, err, tmp1, 1, i, sc_desc, 1)
152 : #else /* WITH_MPI */
153 944800 : err = PRECISION_NRM2(na,tmp1(1,i),1)
154 : #endif /* WITH_MPI */
155 2863200 : errmax = max(errmax, err)
156 : #endif /* REALCASE */
157 :
158 : #if COMPLEXCASE == 1
159 2337600 : xc = 0
160 : #ifdef WITH_MPI
161 1568000 : call scal_PRECISION_DOTC(na, xc, tmp1, 1, i, sc_desc, 1, tmp1, 1, i, sc_desc, 1)
162 : #else /* WITH_MPI */
163 769600 : xc = PRECISION_DOTC(na,tmp1,1,tmp1,1)
164 : #endif /* WITH_MPI */
165 2337600 : errmax = max(errmax, sqrt(real(xc,kind=REAL_DATATYPE)))
166 : #endif /* COMPLEXCASE */
167 : enddo
168 :
169 : ! Get maximum error norm over all processors
170 31536 : err = errmax
171 : #ifdef WITH_MPI
172 21152 : call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
173 : #else /* WITH_MPI */
174 10384 : errmax = err
175 : #endif /* WITH_MPI */
176 31536 : if (myid==0) print *,'Results of numerical residual checks:'
177 31536 : if (myid==0) print *,'Error Residual :',errmax
178 31536 : if (nev .ge. 2) then
179 31536 : if (errmax .gt. tol_res .or. errmax .eq. 0.0_rk) then
180 0 : status = 1
181 : endif
182 : else
183 0 : if (errmax .gt. tol_res) then
184 0 : status = 1
185 : endif
186 : endif
187 :
188 : ! 2. Eigenvector orthogonality
189 : !TODO for the generalized eigenvector problem, the orthogonality test has to be altered
190 : !TODO at the moment, it is skipped
191 31536 : if(.not. present(bs)) then
192 : ! tmp1 = Z**T * Z
193 31536 : tmp1 = 0
194 : #ifdef WITH_MPI
195 : call scal_PRECISION_GEMM(BLAS_TRANS_OR_CONJ, 'N', nev, nev, na, ONE, z, 1, 1, &
196 21152 : sc_desc, z, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
197 : #else /* WITH_MPI */
198 10384 : call PRECISION_GEMM(BLAS_TRANS_OR_CONJ,'N',nev,nev,na,ONE,z,na,z,na,ZERO,tmp1,na)
199 : #endif /* WITH_MPI */
200 : !TODO for the C interface, not all information is passed (zeros instead)
201 : !TODO than this part of the test cannot be done
202 : !TODO either we will not have this part of test at all, or it will be improved
203 31536 : if(nblk > 0) then
204 : ! First check, whether the elements on diagonal are 1 .. "normality" of the vectors
205 28176 : err = 0.0_rk
206 4254576 : do i=1, nev
207 4226400 : if (map_global_array_index_to_local_index(i, i, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
208 2808000 : err = max(err, abs(tmp1(rowLocal,colLocal) - 1.0_rk))
209 : endif
210 : end do
211 : #ifdef WITH_MPI
212 18912 : call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
213 : #else /* WITH_MPI */
214 9264 : errmax = err
215 : #endif /* WITH_MPI */
216 28176 : if (myid==0) print *,'Maximal error in eigenvector lengths:',errmax
217 : end if
218 :
219 : ! Second, find the maximal error in the whole Z**T * Z matrix (its diference from identity matrix)
220 : ! Initialize tmp2 to unit matrix
221 31536 : tmp2 = 0
222 : #ifdef WITH_MPI
223 21152 : call scal_PRECISION_LASET('A', nev, nev, ZERO, ONE, tmp2, 1, 1, sc_desc)
224 : #else /* WITH_MPI */
225 10384 : call PRECISION_LASET('A',nev,nev,ZERO,ONE,tmp2,na)
226 : #endif /* WITH_MPI */
227 :
228 : ! ! tmp1 = Z**T * Z - Unit Matrix
229 31536 : tmp1(:,:) = tmp1(:,:) - tmp2(:,:)
230 :
231 : ! Get maximum error (max abs value in tmp1)
232 31536 : err = maxval(abs(tmp1))
233 : #ifdef WITH_MPI
234 21152 : call mpi_allreduce(err, errmax, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
235 : #else /* WITH_MPI */
236 10384 : errmax = err
237 : #endif /* WITH_MPI */
238 31536 : if (myid==0) print *,'Error Orthogonality:',errmax
239 31536 : if (nev .ge. 2) then
240 31536 : if (errmax .gt. tol_orth .or. errmax .eq. 0.0_rk) then
241 0 : status = 1
242 : endif
243 : else
244 0 : if (errmax .gt. tol_orth) then
245 0 : status = 1
246 : endif
247 : endif
248 : endif ! skiping test of orthogonality for generalized eigenproblem
249 63072 : end function
250 :
251 :
252 : #if REALCASE == 1
253 :
254 : #ifdef DOUBLE_PRECISION_REAL
255 : !c> int check_correctness_evp_numeric_residuals_real_double_f(int na, int nev, int na_rows, int na_cols,
256 : !c> double *as, double *z, double *ev,
257 : !c> int sc_desc[9], int myid);
258 : #else
259 : !c> int check_correctness_evp_numeric_residuals_real_single_f(int na, int nev, int na_rows, int na_cols,
260 : !c> float *as, float *z, float *ev,
261 : !c> int sc_desc[9], int myid);
262 : #endif
263 :
264 : #endif /* REALCASE */
265 :
266 : #if COMPLEXCASE == 1
267 : #ifdef DOUBLE_PRECISION_COMPLEX
268 : !c> int check_correctness_evp_numeric_residuals_complex_double_f(int na, int nev, int na_rows, int na_cols,
269 : !c> complex double *as, complex double *z, double *ev,
270 : !c> int sc_desc[9], int myid);
271 : #else
272 : !c> int check_correctness_evp_numeric_residuals_complex_single_f(int na, int nev, int na_rows, int na_cols,
273 : !c> complex float *as, complex float *z, float *ev,
274 : !c> int sc_desc[9], int myid);
275 : #endif
276 : #endif /* COMPLEXCASE */
277 :
278 3360 : function check_correctness_evp_numeric_residuals_&
279 : &MATH_DATATYPE&
280 : &_&
281 : &PRECISION&
282 3360 : &_f (na, nev, na_rows, na_cols, as, z, ev, sc_desc, myid) result(status) &
283 : bind(C,name="check_correctness_evp_numeric_residuals_&
284 : &MATH_DATATYPE&
285 : &_&
286 : &PRECISION&
287 : &_f")
288 :
289 : use iso_c_binding
290 :
291 : implicit none
292 : #include "../../src/general/precision_kinds.F90"
293 :
294 : integer(kind=c_int) :: status
295 : integer(kind=c_int), value :: na, nev, myid, na_rows, na_cols
296 : MATH_DATATYPE(kind=rck) :: as(1:na_rows,1:na_cols), z(1:na_rows,1:na_cols)
297 : real(kind=rck) :: ev(1:na)
298 : integer(kind=c_int) :: sc_desc(1:9)
299 :
300 : ! TODO: I did not want to add all the variables to the C interface as well
301 : ! TODO: I think that we should find a better way to pass this information
302 : ! TODO: to all the functions anyway (get it from sc_desc, pass elpa_t, etc..)
303 : status = check_correctness_evp_numeric_residuals_&
304 : &MATH_DATATYPE&
305 : &_&
306 : &PRECISION&
307 3360 : & (na, nev, as, z, ev, sc_desc, 0, myid, 0, 0, 0, 0)
308 :
309 3360 : end function
310 :
311 1440 : function check_correctness_eigenvalues_toeplitz_&
312 : &MATH_DATATYPE&
313 : &_&
314 : &PRECISION&
315 1440 : & (na, diagonalElement, subdiagonalElement, ev, z, myid) result(status)
316 : use iso_c_binding
317 : implicit none
318 : #include "../../src/general/precision_kinds.F90"
319 :
320 : integer :: status, ii, j, myid
321 : integer, intent(in) :: na
322 : real(kind=rck) :: diagonalElement, subdiagonalElement
323 2880 : real(kind=rck) :: ev_analytic(na), ev(na)
324 : MATH_DATATYPE(kind=rck) :: z(:,:)
325 :
326 : #if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
327 : real(kind=rck), parameter :: pi = 3.141592653589793238462643383279_c_double
328 : #else
329 : real(kind=rck), parameter :: pi = 3.1415926535897932_c_float
330 : #endif
331 : real(kind=rck) :: tmp, maxerr
332 : integer :: loctmp
333 1440 : status = 0
334 :
335 : ! analytic solution
336 217440 : do ii=1, na
337 : ev_analytic(ii) = diagonalElement + 2.0_rk * &
338 : subdiagonalElement *cos( pi*real(ii,kind=rk)/ &
339 216000 : real(na+1,kind=rk) )
340 : enddo
341 :
342 : ! sort analytic solution:
343 :
344 : ! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
345 : ! a proper sorting algorithmus might be implemented here
346 :
347 1440 : tmp = minval(ev_analytic)
348 1440 : loctmp = minloc(ev_analytic, 1)
349 :
350 1440 : ev_analytic(loctmp) = ev_analytic(1)
351 1440 : ev_analytic(1) = tmp
352 216000 : do ii=2, na
353 214560 : tmp = ev_analytic(ii)
354 16306560 : do j= ii, na
355 16092000 : if (ev_analytic(j) .lt. tmp) then
356 7885440 : tmp = ev_analytic(j)
357 7885440 : loctmp = j
358 : endif
359 : enddo
360 214560 : ev_analytic(loctmp) = ev_analytic(ii)
361 214560 : ev_analytic(ii) = tmp
362 : enddo
363 :
364 : ! compute a simple error max of eigenvalues
365 1440 : maxerr = 0.0
366 1440 : maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)
367 :
368 : #if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
369 960 : if (maxerr .gt. 8.e-13_c_double .or. maxerr .eq. 0.0_c_double) then
370 : #else
371 480 : if (maxerr .gt. 8.e-4_c_float .or. maxerr .eq. 0.0_c_float) then
372 : #endif
373 0 : status = 1
374 0 : if (myid .eq. 0) then
375 0 : print *,"Result of Toeplitz matrix test: "
376 0 : print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
377 : endif
378 : endif
379 :
380 1440 : if (status .eq. 0) then
381 1440 : if (myid .eq. 0) then
382 960 : print *,"Result of Toeplitz matrix test: test passed"
383 960 : print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
384 : endif
385 : endif
386 2880 : end function
387 :
388 576 : function check_correctness_cholesky_&
389 : &MATH_DATATYPE&
390 : &_&
391 : &PRECISION&
392 576 : & (na, a, as, na_rows, sc_desc, myid) result(status)
393 : implicit none
394 : #include "../../src/general/precision_kinds.F90"
395 : integer(kind=ik) :: status
396 : integer(kind=ik), intent(in) :: na, myid, na_rows
397 :
398 : MATH_DATATYPE(kind=rck), intent(in) :: a(:,:), as(:,:)
399 1152 : MATH_DATATYPE(kind=rck), dimension(size(as,dim=1),size(as,dim=2)) :: tmp1, tmp2
400 : real(kind=rk) :: norm, normmax
401 :
402 : #ifdef WITH_MPI
403 : real(kind=rck) :: p&
404 : &BLAS_CHAR&
405 : &lange
406 : #else /* WITH_MPI */
407 : real(kind=rck) :: BLAS_CHAR&
408 : &lange
409 : #endif /* WITH_MPI */
410 :
411 : integer(kind=ik) :: sc_desc(:)
412 : real(kind=rck) :: err, errmax
413 : integer :: mpierr
414 :
415 576 : status = 0
416 576 : tmp1(:,:) = 0.0_rck
417 :
418 :
419 : #if REALCASE == 1
420 : ! tmp1 = a**T
421 : #ifdef WITH_MPI
422 : call p&
423 : &BLAS_CHAR&
424 192 : &tran(na, na, 1.0_rck, a, 1, 1, sc_desc, 0.0_rck, tmp1, 1, 1, sc_desc)
425 : #else /* WITH_MPI */
426 96 : tmp1 = transpose(a)
427 : #endif /* WITH_MPI */
428 : #endif /* REALCASE == 1 */
429 :
430 : #if COMPLEXCASE == 1
431 : ! tmp1 = a**H
432 : #ifdef WITH_MPI
433 : call p&
434 : &BLAS_CHAR&
435 192 : &tranc(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
436 : #else /* WITH_MPI */
437 96 : tmp1 = transpose(conjg(a))
438 : #endif /* WITH_MPI */
439 : #endif /* COMPLEXCASE == 1 */
440 :
441 : ! tmp2 = a**T * a
442 : #ifdef WITH_MPI
443 : call p&
444 : &BLAS_CHAR&
445 : &gemm("N","N", na, na, na, ONE, tmp1, 1, 1, sc_desc, a, 1, 1, &
446 384 : sc_desc, ZERO, tmp2, 1, 1, sc_desc)
447 : #else /* WITH_MPI */
448 : call BLAS_CHAR&
449 192 : &gemm("N","N", na, na, na, ONE, tmp1, na, a, na, ZERO, tmp2, na)
450 : #endif /* WITH_MPI */
451 :
452 : ! compare tmp2 with original matrix
453 576 : tmp2(:,:) = tmp2(:,:) - as(:,:)
454 :
455 : #ifdef WITH_MPI
456 : norm = p&
457 : &BLAS_CHAR&
458 384 : &lange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
459 : #else /* WITH_MPI */
460 : norm = BLAS_CHAR&
461 192 : &lange("M", na, na, tmp2, na_rows, tmp1)
462 : #endif /* WITH_MPI */
463 :
464 :
465 : #ifdef WITH_MPI
466 384 : call mpi_allreduce(norm,normmax,1,MPI_REAL_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr)
467 : #else /* WITH_MPI */
468 192 : normmax = norm
469 : #endif /* WITH_MPI */
470 :
471 576 : if (myid .eq. 0) then
472 384 : print *," Maximum error of result: ", normmax
473 : endif
474 :
475 : #if REALCASE == 1
476 : #ifdef DOUBLE_PRECISION_REAL
477 : ! if (normmax .gt. 5e-12_rk8 .or. normmax .eq. 0.0_rk8) then
478 192 : if (normmax .gt. 5e-12_rk8) then
479 0 : status = 1
480 : endif
481 : #else
482 : ! if (normmax .gt. 5e-4_rk4 .or. normmax .eq. 0.0_rk4) then
483 96 : if (normmax .gt. 5e-4_rk4 ) then
484 0 : status = 1
485 : endif
486 : #endif
487 : #endif
488 :
489 : #if COMPLEXCASE == 1
490 : #ifdef DOUBLE_PRECISION_COMPLEX
491 : ! if (normmax .gt. 5e-11_rk8 .or. normmax .eq. 0.0_rk8) then
492 192 : if (normmax .gt. 5e-11_rk8 ) then
493 0 : status = 1
494 : endif
495 : #else
496 : ! if (normmax .gt. 5e-3_rk4 .or. normmax .eq. 0.0_rk4) then
497 96 : if (normmax .gt. 5e-3_rk4) then
498 0 : status = 1
499 : endif
500 : #endif
501 : #endif
502 1152 : end function
503 :
504 768 : function check_correctness_hermitian_multiply_&
505 : &MATH_DATATYPE&
506 : &_&
507 : &PRECISION&
508 768 : & (na, a, b, c, na_rows, sc_desc, myid) result(status)
509 : implicit none
510 : #include "../../src/general/precision_kinds.F90"
511 : integer(kind=ik) :: status
512 : integer(kind=ik), intent(in) :: na, myid, na_rows
513 : MATH_DATATYPE(kind=rck), intent(in) :: a(:,:), b(:,:), c(:,:)
514 1536 : MATH_DATATYPE(kind=rck), dimension(size(a,dim=1),size(a,dim=2)) :: tmp1, tmp2
515 : real(kind=rck) :: norm, normmax
516 :
517 : #ifdef WITH_MPI
518 : real(kind=rck) :: p&
519 : &BLAS_CHAR&
520 : &lange
521 : #else /* WITH_MPI */
522 : real(kind=rck) :: BLAS_CHAR&
523 : &lange
524 : #endif /* WITH_MPI */
525 :
526 : integer(kind=ik) :: sc_desc(:)
527 : real(kind=rck) :: err, errmax
528 : integer :: mpierr
529 :
530 768 : status = 0
531 768 : tmp1(:,:) = ZERO
532 :
533 : #if REALCASE == 1
534 : ! tmp1 = a**T
535 : #ifdef WITH_MPI
536 : call p&
537 : &BLAS_CHAR&
538 320 : &tran(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
539 : #else /* WITH_MPI */
540 160 : tmp1 = transpose(a)
541 : #endif /* WITH_MPI */
542 :
543 : #endif /* REALCASE == 1 */
544 :
545 : #if COMPLEXCASE == 1
546 : ! tmp1 = a**H
547 : #ifdef WITH_MPI
548 : call p&
549 : &BLAS_CHAR&
550 192 : &tranc(na, na, ONE, a, 1, 1, sc_desc, ZERO, tmp1, 1, 1, sc_desc)
551 : #else /* WITH_MPI */
552 96 : tmp1 = transpose(conjg(a))
553 : #endif /* WITH_MPI */
554 : #endif /* COMPLEXCASE == 1 */
555 :
556 : ! tmp2 = tmp1 * b
557 : #ifdef WITH_MPI
558 : call p&
559 : &BLAS_CHAR&
560 : &gemm("N","N", na, na, na, ONE, tmp1, 1, 1, sc_desc, b, 1, 1, &
561 512 : sc_desc, ZERO, tmp2, 1, 1, sc_desc)
562 : #else
563 : call BLAS_CHAR&
564 256 : &gemm("N","N", na, na, na, ONE, tmp1, na, b, na, ZERO, tmp2, na)
565 : #endif
566 :
567 : ! compare tmp2 with c
568 768 : tmp2(:,:) = tmp2(:,:) - c(:,:)
569 :
570 : #ifdef WITH_MPI
571 : norm = p&
572 : &BLAS_CHAR&
573 512 : &lange("M",na, na, tmp2, 1, 1, sc_desc, tmp1)
574 : #else /* WITH_MPI */
575 : norm = BLAS_CHAR&
576 256 : &lange("M", na, na, tmp2, na_rows, tmp1)
577 : #endif /* WITH_MPI */
578 :
579 : #ifdef WITH_MPI
580 512 : call mpi_allreduce(norm,normmax,1,MPI_REAL_PRECISION,MPI_MAX,MPI_COMM_WORLD,mpierr)
581 : #else /* WITH_MPI */
582 256 : normmax = norm
583 : #endif /* WITH_MPI */
584 :
585 768 : if (myid .eq. 0) then
586 512 : print *," Maximum error of result: ", normmax
587 : endif
588 :
589 : #ifdef DOUBLE_PRECISION_REAL
590 384 : if (normmax .gt. 5e-11_rk8 ) then
591 0 : status = 1
592 : endif
593 : #else
594 384 : if (normmax .gt. 5e-3_rk4 ) then
595 0 : status = 1
596 : endif
597 : #endif
598 :
599 : #ifdef DOUBLE_PRECISION_COMPLEX
600 384 : if (normmax .gt. 5e-11_rk8 ) then
601 0 : status = 1
602 : endif
603 : #else
604 384 : if (normmax .gt. 5e-3_rk4 ) then
605 0 : status = 1
606 : endif
607 : #endif
608 1536 : end function
609 :
610 3360 : function check_correctness_eigenvalues_frank_&
611 : &MATH_DATATYPE&
612 : &_&
613 : &PRECISION&
614 3360 : & (na, ev, z, myid) result(status)
615 : use iso_c_binding
616 : implicit none
617 : #include "../../src/general/precision_kinds.F90"
618 :
619 : integer :: status, i, j, myid
620 : integer, intent(in) :: na
621 6720 : real(kind=rck) :: ev_analytic(na), ev(na)
622 : MATH_DATATYPE(kind=rck) :: z(:,:)
623 :
624 : #if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
625 : real(kind=rck), parameter :: pi = 3.141592653589793238462643383279_c_double
626 : #else
627 : real(kind=rck), parameter :: pi = 3.1415926535897932_c_float
628 : #endif
629 : real(kind=rck) :: tmp, maxerr
630 : integer :: loctmp
631 3360 : status = 0
632 :
633 : ! analytic solution
634 507360 : do i = 1, na
635 504000 : j = na - i
636 : #if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
637 : ev_analytic(i) = pi * (2.0_c_double * real(j,kind=c_double) + 1.0_c_double) / &
638 504000 : (2.0_c_double * real(na,kind=c_double) + 1.0_c_double)
639 504000 : ev_analytic(i) = 0.5_c_double / (1.0_c_double - cos(ev_analytic(i)))
640 : #else
641 : ev_analytic(i) = pi * (2.0_c_float * real(j,kind=c_float) + 1.0_c_float) / &
642 0 : (2.0_c_float * real(na,kind=c_float) + 1.0_c_float)
643 0 : ev_analytic(i) = 0.5_c_float / (1.0_c_float - cos(ev_analytic(i)))
644 : #endif
645 : enddo
646 :
647 : ! sort analytic solution:
648 :
649 : ! this hack is neither elegant, nor optimized: for huge matrixes it might be expensive
650 : ! a proper sorting algorithmus might be implemented here
651 :
652 3360 : tmp = minval(ev_analytic)
653 3360 : loctmp = minloc(ev_analytic, 1)
654 :
655 3360 : ev_analytic(loctmp) = ev_analytic(1)
656 3360 : ev_analytic(1) = tmp
657 504000 : do i=2, na
658 500640 : tmp = ev_analytic(i)
659 38048640 : do j= i, na
660 37548000 : if (ev_analytic(j) .lt. tmp) then
661 0 : tmp = ev_analytic(j)
662 0 : loctmp = j
663 : endif
664 : enddo
665 500640 : ev_analytic(loctmp) = ev_analytic(i)
666 500640 : ev_analytic(i) = tmp
667 : enddo
668 :
669 : ! compute a simple error max of eigenvalues
670 3360 : maxerr = 0.0
671 3360 : maxerr = maxval( (ev(:) - ev_analytic(:))/ev_analytic(:) , 1)
672 :
673 : #if defined(DOUBLE_PRECISION_REAL) || defined(DOUBLE_PRECISION_COMPLEX)
674 3360 : if (maxerr .gt. 8.e-13_c_double) then
675 : #else
676 0 : if (maxerr .gt. 8.e-4_c_float) then
677 : #endif
678 0 : status = 1
679 0 : if (myid .eq. 0) then
680 0 : print *,"Result of Frank matrix test: "
681 0 : print *,"Eigenvalues differ from analytic solution: maxerr = ",maxerr
682 : endif
683 : endif
684 6720 : end function
685 :
686 : ! vim: syntax=fortran
|
