Line data Source code
1 : ! (c) Copyright Pavel Kus, 2017, MPCDF
2 : !
3 : ! This file is part of ELPA.
4 : !
5 : ! The ELPA library was originally created by the ELPA consortium,
6 : ! consisting of the following organizations:
7 : !
8 : ! - Max Planck Computing and Data Facility (MPCDF), formerly known as
9 : ! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
10 : ! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
11 : ! Informatik,
12 : ! - Technische Universität München, Lehrstuhl für Informatik mit
13 : ! Schwerpunkt Wissenschaftliches Rechnen ,
14 : ! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
15 : ! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
16 : ! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
17 : ! and
18 : ! - IBM Deutschland GmbH
19 : !
20 : !
21 : ! More information can be found here:
22 : ! http://elpa.mpcdf.mpg.de/
23 : !
24 : ! ELPA is free software: you can redistribute it and/or modify
25 : ! it under the terms of the version 3 of the license of the
26 : ! GNU Lesser General Public License as published by the Free
27 : ! Software Foundation.
28 : !
29 : ! ELPA is distributed in the hope that it will be useful,
30 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
31 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 : ! GNU Lesser General Public License for more details.
33 : !
34 : ! You should have received a copy of the GNU Lesser General Public License
35 : ! along with ELPA. If not, see <http://www.gnu.org/licenses/>
36 : !
37 : ! ELPA reflects a substantial effort on the part of the original
38 : ! ELPA consortium, and we ask you to respect the spirit of the
39 : ! license that we chose: i.e., please contribute any changes you
40 : ! may have back to the original ELPA library distribution, and keep
41 : ! any derivatives of ELPA under the same license that we chose for
42 : ! the original distribution, the GNU Lesser General Public License.
43 :
44 : subroutine prepare_matrix_analytic_&
45 : &MATH_DATATYPE&
46 : &_&
47 2112 : &PRECISION&
48 2112 : &(na, a, nblk, myid, np_rows, np_cols, my_prow, my_pcol)
49 : implicit none
50 : integer(kind=ik), intent(in) :: na, nblk, myid, np_rows, np_cols, my_prow, my_pcol
51 : MATH_DATATYPE(kind=REAL_DATATYPE), intent(inout) :: a(:,:)
52 :
53 : integer(kind=ik) :: globI, globJ, locI, locJ, levels(num_primes)
54 :
55 : type(timer_t) :: timer
56 :
57 2112 : call timer%enable()
58 2112 : call timer%start("prepare_matrix_analytic")
59 :
60 : ! for debug only, do it systematicaly somehow ... unit tests
61 : call check_module_sanity_&
62 : &MATH_DATATYPE&
63 : &_&
64 : &PRECISION&
65 2112 : &(myid)
66 :
67 2112 : if(.not. decompose(na, levels)) then
68 0 : if(myid == 0) then
69 0 : print *, "Analytic test can be run only with matrix sizes of the form 2^n * 3^m * 5^o"
70 0 : stop 1
71 : end if
72 : end if
73 :
74 2112 : call timer%start("loop")
75 318912 : do globI = 1, na
76 47836800 : do globJ = 1, na
77 47520000 : if(map_global_array_index_to_local_index(globI, globJ, locI, locJ, &
78 : nblk, np_rows, np_cols, my_prow, my_pcol)) then
79 30240000 : call timer%start("evaluation")
80 : a(locI, locJ) = analytic_matrix_&
81 : &MATH_DATATYPE&
82 : &_&
83 : &PRECISION&
84 30240000 : &(na, globI, globJ)
85 30240000 : call timer%stop("evaluation")
86 : end if
87 : end do
88 : end do
89 2112 : call timer%stop("loop")
90 :
91 2112 : call timer%stop("prepare_matrix_analytic")
92 2112 : if(myid == 0) then
93 1344 : call timer%print("prepare_matrix_analytic")
94 : end if
95 2112 : call timer%free()
96 4224 : end subroutine
97 :
98 16608 : function check_correctness_analytic_&
99 : &MATH_DATATYPE&
100 : &_&
101 : &PRECISION&
102 16608 : &(na, nev, ev, z, nblk, myid, np_rows, np_cols, my_prow, my_pcol, check_all_evals, check_eigenvectors) result(status)
103 : implicit none
104 : #include "../../src/general/precision_kinds.F90"
105 : integer(kind=ik), intent(in) :: na, nev, nblk, myid, np_rows, &
106 : np_cols, my_prow, my_pcol
107 : integer(kind=ik) :: status, mpierr
108 : MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:)
109 : real(kind=rk), intent(inout) :: ev(:)
110 : logical, intent(in) :: check_all_evals, check_eigenvectors
111 :
112 : integer(kind=ik) :: globI, globJ, locI, locJ, &
113 : levels(num_primes)
114 : real(kind=rk) :: diff, max_z_diff, max_ev_diff, &
115 : glob_max_z_diff, max_curr_z_diff
116 : #ifdef DOUBLE_PRECISION
117 : real(kind=rk), parameter :: tol_eigenvalues = 5e-14_rk8
118 : real(kind=rk), parameter :: tol_eigenvectors = 6e-11_rk8
119 : #endif
120 : #ifdef SINGLE_PRECISION
121 : ! tolerance needs to be very high due to qr tests
122 : ! it should be distinguished somehow!
123 : real(kind=rk), parameter :: tol_eigenvalues = 7e-6_rk4
124 : real(kind=rk), parameter :: tol_eigenvectors = 4e-3_rk4
125 : #endif
126 : real(kind=rk) :: computed_ev, expected_ev
127 : MATH_DATATYPE(kind=rck) :: computed_z, expected_z
128 :
129 : MATH_DATATYPE(kind=rck) :: max_value_for_normalization, &
130 : computed_z_on_max_position, &
131 : normalization_quotient
132 22656 : MATH_DATATYPE(kind=rck) :: max_values_array(np_rows * np_cols), &
133 : corresponding_exact_value
134 : integer(kind=ik) :: max_value_idx, rank_with_max, &
135 : rank_with_max_reduced, &
136 : num_checked_evals
137 22656 : integer(kind=ik) :: max_idx_array(np_rows * np_cols), &
138 : rank
139 :
140 : type(timer_t) :: timer
141 :
142 16608 : call timer%enable()
143 16608 : call timer%start("check_correctness_analytic")
144 :
145 :
146 16608 : if(.not. decompose(na, levels)) then
147 0 : print *, "can not decomopse matrix size"
148 0 : stop 1
149 : end if
150 :
151 16608 : if(check_all_evals) then
152 16224 : num_checked_evals = na
153 : else
154 384 : num_checked_evals = nev
155 : endif
156 : !call print_matrix(myid, na, z, "z")
157 16608 : max_z_diff = 0.0_rk
158 16608 : max_ev_diff = 0.0_rk
159 16608 : call timer%start("loop_eigenvalues")
160 2507808 : do globJ = 1, num_checked_evals
161 2491200 : computed_ev = ev(globJ)
162 2491200 : call timer%start("evaluation")
163 : expected_ev = analytic_eigenvalues_real_&
164 : &PRECISION&
165 2491200 : &(na, globJ)
166 2491200 : call timer%stop("evaluation")
167 2491200 : diff = abs(computed_ev - expected_ev)
168 2491200 : max_ev_diff = max(diff, max_ev_diff)
169 : end do
170 16608 : call timer%stop("loop_eigenvalues")
171 :
172 16608 : call timer%start("loop_eigenvectors")
173 2507808 : do globJ = 1, nev
174 2491200 : max_curr_z_diff = 0.0_rk
175 :
176 : ! eigenvectors are unique up to multiplication by scalar (complex in complex case)
177 : ! to be able to compare them with analytic, we have to normalize them somehow
178 : ! we will find a value in computed eigenvector with highest absolut value and enforce
179 : ! such multiple of computed eigenvector, that the value on corresponding position is the same
180 : ! as an corresponding value in the analytical eigenvector
181 :
182 : ! find the maximal value in the local part of given eigenvector (with index globJ)
183 2491200 : max_value_for_normalization = 0.0_rk
184 2491200 : max_value_idx = -1
185 376171200 : do globI = 1, na
186 373680000 : if(map_global_array_index_to_local_index(globI, globJ, locI, locJ, &
187 : nblk, np_rows, np_cols, my_prow, my_pcol)) then
188 246240000 : computed_z = z(locI, locJ)
189 246240000 : if(abs(computed_z) > abs(max_value_for_normalization)) then
190 17599984 : max_value_for_normalization = computed_z
191 17599984 : max_value_idx = globI
192 : end if
193 : end if
194 : end do
195 :
196 : ! find the global maximum and its position. From technical reasons (looking for a
197 : ! maximum of complex number), it is not so easy to do it nicely. Therefore we
198 : ! communicate local maxima to mpi rank 0 and resolve there. If we wanted to do
199 : ! it without this, it would be tricky.. question of uniquness - two complex numbers
200 : ! with the same absolut values, but completely different...
201 : #ifdef WITH_MPI
202 : call MPI_Gather(max_value_for_normalization, 1, MPI_MATH_DATATYPE_PRECISION, &
203 1699200 : max_values_array, 1, MPI_MATH_DATATYPE_PRECISION, 0, MPI_COMM_WORLD, mpierr)
204 1699200 : call MPI_Gather(max_value_idx, 1, MPI_INT, max_idx_array, 1, MPI_INT, 0, MPI_COMM_WORLD, mpierr)
205 1699200 : max_value_for_normalization = 0.0_rk
206 1699200 : max_value_idx = -1
207 5097600 : do rank = 1, np_cols * np_rows
208 3398400 : if(abs(max_values_array(rank)) > abs(max_value_for_normalization)) then
209 2426448 : max_value_for_normalization = max_values_array(rank)
210 2426448 : max_value_idx = max_idx_array(rank)
211 : end if
212 : end do
213 1699200 : call MPI_Bcast(max_value_for_normalization, 1, MPI_MATH_DATATYPE_PRECISION, 0, MPI_COMM_WORLD, mpierr)
214 1699200 : call MPI_Bcast(max_value_idx, 1, MPI_INT, 0, MPI_COMM_WORLD, mpierr)
215 : #endif
216 : ! we decided what the maximum computed value is. Calculate expected value on the same
217 2491200 : if(abs(max_value_for_normalization) < 0.0001_rk) then
218 0 : if(myid == 0) print *, 'Maximal value in eigenvector too small :', max_value_for_normalization
219 0 : status =1
220 0 : return
221 : end if
222 2491200 : call timer%start("evaluation_helper")
223 : corresponding_exact_value = analytic_eigenvectors_&
224 : &MATH_DATATYPE&
225 : &_&
226 : &PRECISION&
227 2491200 : &(na, max_value_idx, globJ)
228 2491200 : call timer%stop("evaluation_helper")
229 2491200 : normalization_quotient = corresponding_exact_value / max_value_for_normalization
230 : ! write(*,*) "normalization q", normalization_quotient
231 :
232 : ! compare computed and expected eigenvector values, but take into account normalization quotient
233 376171200 : do globI = 1, na
234 373680000 : if(map_global_array_index_to_local_index(globI, globJ, locI, locJ, &
235 : nblk, np_rows, np_cols, my_prow, my_pcol)) then
236 246240000 : computed_z = z(locI, locJ)
237 246240000 : call timer%start("evaluation")
238 : expected_z = analytic_eigenvectors_&
239 : &MATH_DATATYPE&
240 : &_&
241 : &PRECISION&
242 246240000 : &(na, globI, globJ)
243 246240000 : call timer%stop("evaluation")
244 246240000 : max_curr_z_diff = max(abs(normalization_quotient * computed_z - expected_z), max_curr_z_diff)
245 : end if
246 : end do
247 : ! we have max difference of one of the eigenvectors, update global
248 2491200 : max_z_diff = max(max_z_diff, max_curr_z_diff)
249 : end do !globJ
250 16608 : call timer%stop("loop_eigenvectors")
251 :
252 : #ifdef WITH_MPI
253 11328 : call mpi_allreduce(max_z_diff, glob_max_z_diff, 1, MPI_REAL_PRECISION, MPI_MAX, MPI_COMM_WORLD, mpierr)
254 : #else
255 5280 : glob_max_z_diff = max_z_diff
256 : #endif
257 16608 : if(myid == 0) print *, 'Maximum error in eigenvalues :', max_ev_diff
258 16608 : if (check_eigenvectors) then
259 8304 : if(myid == 0) print *, 'Maximum error in eigenvectors :', glob_max_z_diff
260 : endif
261 :
262 16608 : status = 0
263 16608 : if (nev .gt. 2) then
264 16608 : if (max_ev_diff .gt. tol_eigenvalues .or. max_ev_diff .eq. 0.0_rk) status = 1
265 16608 : if (check_eigenvectors) then
266 8304 : if (glob_max_z_diff .gt. tol_eigenvectors .or. glob_max_z_diff .eq. 0.0_rk) status = 1
267 : endif
268 : else
269 0 : if (max_ev_diff .gt. tol_eigenvalues) status = 1
270 0 : if (check_eigenvectors) then
271 0 : if (glob_max_z_diff .gt. tol_eigenvectors) status = 1
272 : endif
273 : endif
274 :
275 16608 : call timer%stop("check_correctness_analytic")
276 16608 : if(myid == 0) then
277 10944 : call timer%print("check_correctness_analytic")
278 : end if
279 16608 : call timer%free()
280 33216 : end function
281 :
282 :
283 79447488 : function analytic_matrix_&
284 : &MATH_DATATYPE&
285 : &_&
286 : &PRECISION&
287 : &(na, i, j) result(element)
288 : implicit none
289 : integer(kind=ik), intent(in) :: na, i, j
290 : MATH_DATATYPE(kind=REAL_DATATYPE) :: element
291 :
292 : element = analytic_&
293 : &MATH_DATATYPE&
294 : &_&
295 : &PRECISION&
296 79447488 : &(na, i, j, ANALYTIC_MATRIX)
297 :
298 79447488 : end function
299 :
300 297938688 : function analytic_eigenvectors_&
301 : &MATH_DATATYPE&
302 : &_&
303 : &PRECISION&
304 : &(na, i, j) result(element)
305 : implicit none
306 : integer(kind=ik), intent(in) :: na, i, j
307 : MATH_DATATYPE(kind=REAL_DATATYPE) :: element
308 :
309 : element = analytic_&
310 : &MATH_DATATYPE&
311 : &_&
312 : &PRECISION&
313 297938688 : &(na, i, j, ANALYTIC_EIGENVECTORS)
314 :
315 297938688 : end function
316 :
317 2491200 : function analytic_eigenvalues_&
318 : &MATH_DATATYPE&
319 : &_&
320 : &PRECISION&
321 : &(na, i) result(element)
322 : implicit none
323 : integer(kind=ik), intent(in) :: na, i
324 : real(kind=REAL_DATATYPE) :: element
325 :
326 : element = analytic_real_&
327 : &PRECISION&
328 2491200 : &(na, i, i, ANALYTIC_EIGENVALUES)
329 :
330 2491200 : end function
331 :
332 429084864 : function analytic_&
333 : &MATH_DATATYPE&
334 : &_&
335 : &PRECISION&
336 : &(na, i, j, what) result(element)
337 : implicit none
338 : #include "../../src/general/precision_kinds.F90"
339 : integer(kind=ik), intent(in) :: na, i, j, what
340 : MATH_DATATYPE(kind=rck) :: element, mat2x2(2,2), mat(5,5)
341 : real(kind=rk) :: a, am, amp
342 : integer(kind=ik) :: levels(num_primes)
343 : integer(kind=ik) :: ii, jj, m, prime_id, prime, total_level, level
344 :
345 : real(kind=rk), parameter :: s = 0.5_rk
346 : real(kind=rk), parameter :: c = 0.86602540378443864679_rk
347 : real(kind=rk), parameter :: sq2 = 1.4142135623730950488_rk
348 :
349 : real(kind=rk), parameter :: largest_ev = 2.0_rk
350 :
351 429084864 : assert(i <= na)
352 429084864 : assert(j <= na)
353 429084864 : assert(i >= 0)
354 429084864 : assert(j >= 0)
355 429084864 : assert(decompose(na, levels))
356 : ! go to zero-based indexing
357 429084864 : ii = i - 1
358 429084864 : jj = j - 1
359 429084864 : if (na .gt. 2) then
360 429059520 : a = exp(log(largest_ev)/(na-1))
361 : else
362 25344 : a = exp(log(largest_ev)/(1))
363 : endif
364 :
365 429084864 : element = 1.0_rck
366 : #ifdef COMPLEXCASE
367 198025632 : element = (1.0_rk, 0.0_rk)
368 : #endif
369 429084864 : total_level = 0
370 429084864 : am = a
371 1716339456 : do prime_id = 1,num_primes
372 1287254592 : prime = primes(prime_id)
373 2993228352 : do level = 1, levels(prime_id)
374 1705973760 : amp = am**(prime-1)
375 1705973760 : total_level = total_level + 1
376 1705973760 : if(what == ANALYTIC_MATRIX) then
377 : #ifdef REALCASE
378 : mat2x2 = reshape((/ c*c + amp * s*s, (amp - 1.0_rk) * s*c, &
379 : (amp - 1.0_rk) * s*c, s*s + amp * c*c /), &
380 157167360 : (/2, 2/), order=(/2,1/))
381 : #endif
382 : #ifdef COMPLEXCASE
383 : mat2x2 = reshape((/ 0.5_rck * (amp + 1.0_rck) * (1.0_rk, 0.0_rk), sq2/4.0_rk * (amp - 1.0_rk) * (1.0_rk, 1.0_rk), &
384 : sq2/4.0_rk * (amp - 1.0_rk) * (1.0_rk, -1.0_rk), 0.5_rck * (amp + 1.0_rck) * (1.0_rk, 0.0_rk) /), &
385 157167360 : (/2, 2/), order=(/2,1/))
386 : #endif
387 1391639040 : else if(what == ANALYTIC_EIGENVECTORS) then
388 : #ifdef REALCASE
389 : mat2x2 = reshape((/ c, s, &
390 : -s, c /), &
391 655234560 : (/2, 2/), order=(/2,1/))
392 : #endif
393 : #ifdef COMPLEXCASE
394 : mat2x2 = reshape((/ -sq2/2.0_rck * (1.0_rk, 0.0_rk), -sq2/2.0_rck * (1.0_rk, 0.0_rk), &
395 : 0.5_rk * (1.0_rk, -1.0_rk), 0.5_rk * (-1.0_rk, 1.0_rk) /), &
396 533064960 : (/2, 2/), order=(/2,1/))
397 : #endif
398 203339520 : else if(what == ANALYTIC_EIGENVALUES) then
399 : mat2x2 = reshape((/ 1.0_rck, 0.0_rck, &
400 : 0.0_rck, amp /), &
401 203339520 : (/2, 2/), order=(/2,1/))
402 : else
403 0 : assert(.false.)
404 : end if
405 :
406 1705973760 : mat = 0.0_rck
407 1705973760 : if(prime == 2) then
408 424909440 : mat(1:2, 1:2) = mat2x2
409 1281064320 : else if(prime == 3) then
410 424307520 : mat((/1,3/),(/1,3/)) = mat2x2
411 424307520 : if(what == ANALYTIC_EIGENVECTORS) then
412 296346240 : mat(2,2) = 1.0_rck
413 : else
414 127961280 : mat(2,2) = am
415 : end if
416 856756800 : else if(prime == 5) then
417 856756800 : mat((/1,5/),(/1,5/)) = mat2x2
418 856756800 : if(what == ANALYTIC_EIGENVECTORS) then
419 595406400 : mat(2,2) = 1.0_rck
420 595406400 : mat(3,3) = 1.0_rck
421 595406400 : mat(4,4) = 1.0_rck
422 : else
423 261350400 : mat(2,2) = am
424 261350400 : mat(3,3) = am**2
425 261350400 : mat(4,4) = am**3
426 : end if
427 : else
428 0 : assert(.false.)
429 : end if
430 :
431 : ! write(*,*) "calc value, elem: ", element, ", mat: ", mod(ii,2), mod(jj,2), mat(mod(ii,2), mod(jj,2)), "am ", am
432 : ! write(*,*) " matrix mat", mat
433 1705973760 : element = element * mat(mod(ii,prime) + 1, mod(jj,prime) + 1)
434 1705973760 : ii = ii / prime
435 1705973760 : jj = jj / prime
436 :
437 1705973760 : am = am**prime
438 : end do
439 : end do
440 : !write(*,*) "returning value ", element
441 429084864 : end function
442 :
443 :
444 : subroutine print_matrix_&
445 : &MATH_DATATYPE&
446 : &_&
447 0 : &PRECISION&
448 0 : &(myid, na, mat, mat_name)
449 : implicit none
450 : #include "../../src/general/precision_kinds.F90"
451 : integer(kind=ik), intent(in) :: myid, na
452 : character(len=*), intent(in) :: mat_name
453 : MATH_DATATYPE(kind=rck) :: mat(na, na)
454 : integer(kind=ik) :: i,j
455 : character(len=20) :: na_str
456 :
457 0 : if(myid .ne. 0) &
458 0 : return
459 0 : write(*,*) "Matrix: "//trim(mat_name)
460 0 : write(na_str, *) na
461 0 : do i = 1, na
462 : #ifdef REALCASE
463 0 : write(*, '('//trim(na_str)//'f8.3)') mat(i, :)
464 : #endif
465 : #ifdef COMPLEXCASE
466 0 : write(*,'('//trim(na_str)//'(A,f8.3,A,f8.3,A))') ('(', real(mat(i,j)), ',', aimag(mat(i,j)), ')', j=1,na)
467 : #endif
468 : end do
469 0 : write(*,*)
470 0 : end subroutine
471 :
472 :
473 : subroutine check_matrices_&
474 : &MATH_DATATYPE&
475 : &_&
476 14784 : &PRECISION&
477 : &(myid, na)
478 : implicit none
479 : #include "../../src/general/precision_kinds.F90"
480 : integer(kind=ik), intent(in) :: myid, na
481 29568 : MATH_DATATYPE(kind=rck) :: A(na, na), S(na, na), L(na, na), res(na, na)
482 : integer(kind=ik) :: i, j, decomposition(num_primes)
483 :
484 14784 : assert(decompose(na, decomposition))
485 :
486 439296 : do i = 1, na
487 49632000 : do j = 1, na
488 : A(i,j) = analytic_matrix_&
489 : &MATH_DATATYPE&
490 : &_&
491 : &PRECISION&
492 49207488 : &(na, i, j)
493 : S(i,j) = analytic_eigenvectors_&
494 : &MATH_DATATYPE&
495 : &_&
496 : &PRECISION&
497 49207488 : &(na, i, j)
498 : L(i,j) = analytic_&
499 : &MATH_DATATYPE&
500 : &_&
501 : &PRECISION&
502 49207488 : &(na, i, j, ANALYTIC_EIGENVALUES)
503 : end do
504 : end do
505 :
506 14784 : res = matmul(A,S) - matmul(S,L)
507 : #ifdef DOUBLE_PRECISION
508 9856 : assert(maxval(abs(res)) < 1e-8)
509 : #elif SINGLE_PRECISION
510 4928 : assert(maxval(abs(res)) < 1e-4)
511 : #else
512 : assert(.false.)
513 : #endif
514 : if(.false.) then
515 : !if(na == 2 .or. na == 5) then
516 : call print_matrix(myid, na, A, "A")
517 : call print_matrix(myid, na, S, "S")
518 : call print_matrix(myid, na, L, "L")
519 :
520 : call print_matrix(myid, na, matmul(A,S), "AS")
521 : call print_matrix(myid, na, matmul(S,L), "SL")
522 :
523 : call print_matrix(myid, na, res , "res")
524 : end if
525 :
526 14784 : end subroutine
527 :
528 : subroutine check_module_sanity_&
529 : &MATH_DATATYPE&
530 : &_&
531 2112 : &PRECISION&
532 : &(myid)
533 : implicit none
534 : integer(kind=ik), intent(in) :: myid
535 : integer(kind=ik) :: decomposition(num_primes), i
536 : integer(kind=ik), parameter :: check_sizes(7) = (/2, 3, 5, 6, 10, 25, 150/)
537 2112 : if(myid == 0) print *, "Checking test_analytic module sanity.... "
538 2112 : assert(decompose(1500, decomposition))
539 2112 : assert(all(decomposition == (/2,1,3/)))
540 2112 : assert(decompose(6,decomposition))
541 2112 : assert(all(decomposition == (/1,1,0/)))
542 :
543 16896 : do i =1, size(check_sizes)
544 : call check_matrices_&
545 : &MATH_DATATYPE&
546 : &_&
547 : &PRECISION&
548 14784 : &(myid, check_sizes(i))
549 : end do
550 :
551 2112 : if(myid == 0) print *, "Checking test_analytic module sanity.... DONE"
552 :
553 2112 : end subroutine
|
