Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
math.f90
Go to the documentation of this file.
1! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
2!
3! The UChicago Argonne, LLC as Operator of Argonne National
4! Laboratory holds copyright in the Software. The copyright holder
5! reserves all rights except those expressly granted to licensees,
6! and U.S. Government license rights.
7!
8! Redistribution and use in source and binary forms, with or without
9! modification, are permitted provided that the following conditions
10! are met:
11!
12! 1. Redistributions of source code must retain the above copyright
13! notice, this list of conditions and the disclaimer below.
14!
15! 2. Redistributions in binary form must reproduce the above copyright
16! notice, this list of conditions and the disclaimer (as noted below)
17! in the documentation and/or other materials provided with the
18! distribution.
19!
20! 3. Neither the name of ANL nor the names of its contributors
21! may be used to endorse or promote products derived from this software
22! without specific prior written permission.
23!
24! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
27! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
28! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
29! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
31! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36!
37! Additional BSD Notice
38! ---------------------
39! 1. This notice is required to be provided under our contract with
40! the U.S. Department of Energy (DOE). This work was produced at
41! Argonne National Laboratory under Contract
42! No. DE-AC02-06CH11357 with the DOE.
43!
44! 2. Neither the United States Government nor UCHICAGO ARGONNE,
45! LLC nor any of their employees, makes any warranty,
46! express or implied, or assumes any liability or responsibility for the
47! accuracy, completeness, or usefulness of any information, apparatus,
48! product, or process disclosed, or represents that its use would not
49! infringe privately-owned rights.
50!
51! 3. Also, reference herein to any specific commercial products, process,
52! or services by trade name, trademark, manufacturer or otherwise does
53! not necessarily constitute or imply its endorsement, recommendation,
54! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
55! The views and opinions of authors expressed
56! herein do not necessarily state or reflect those of the United States
57! Government or UCHICAGO ARGONNE, LLC, and shall
58! not be used for advertising or product endorsement purposes.
59!
60module math
61 use num_types, only : rp, dp, sp, qp, i4, xp
63 use mpi_f08, only: mpi_min, mpi_max, mpi_sum, mpi_in_place, mpi_integer, &
64 mpi_allreduce
65 implicit none
66 private
67
69 real(kind=rp), public, parameter :: neko_eps = epsilon(1.0_rp)
70
72 real(kind=rp), public, parameter :: neko_m_ln2 = log(2.0_rp)
73
75 real(kind=rp), public, parameter :: pi = 4._rp*atan(1._rp)
76
77 interface abscmp
78 module procedure sabscmp, dabscmp, qabscmp
79 end interface abscmp
80
81 interface sort
82 module procedure sortrp, sorti4
83 end interface sort
84
85 interface swap
86 module procedure swapdp, swapi4
87 end interface swap
88
89 interface reord
90 module procedure reorddp, reordi4
91 end interface reord
92
93 interface flipv
94 module procedure flipvdp, flipvi4
95 end interface flipv
96
97 interface relcmp
98 module procedure srelcmp, drelcmp, qrelcmp
99 end interface relcmp
100
101 interface pwmax
102 module procedure pwmax_vec2, pwmax_vec3, pwmax_scal2, pwmax_scal3
103 end interface pwmax
104
105 interface pwmin
106 module procedure pwmin_vec2, pwmin_vec3, pwmin_sca2, pwmin_sca3
107 end interface pwmin
108
109 public :: abscmp, rzero, izero, row_zero, rone, copy, cmult, cadd, cfill, &
116
117contains
118
120 pure function sabscmp(x, y)
121 real(kind=sp), intent(in) :: x
122 real(kind=sp), intent(in) :: y
123 logical :: sabscmp
124
125 sabscmp = abs(x - y) .lt. neko_eps
126
127 end function sabscmp
128
130 pure function dabscmp(x, y)
131 real(kind=dp), intent(in) :: x
132 real(kind=dp), intent(in) :: y
133 logical :: dabscmp
134
135 dabscmp = abs(x - y) .lt. neko_eps
136
137 end function dabscmp
138
140 pure function qabscmp(x, y)
141 real(kind=qp), intent(in) :: x
142 real(kind=qp), intent(in) :: y
143 logical :: qabscmp
144
145 qabscmp = abs(x - y) .lt. neko_eps
146
147 end function qabscmp
148
150 pure function srelcmp(x, y, eps)
151 real(kind=sp), intent(in) :: x
152 real(kind=sp), intent(in) :: y
153 real(kind=sp), intent(in), optional :: eps
154 logical :: srelcmp
155 if (present(eps)) then
156 srelcmp = abs(x - y) .le. eps*abs(y)
157 else
158 srelcmp = abs(x - y) .le. neko_eps*abs(y)
159 end if
160
161 end function srelcmp
162
164 pure function drelcmp(x, y, eps)
165 real(kind=dp), intent(in) :: x
166 real(kind=dp), intent(in) :: y
167 real(kind=dp), intent(in), optional :: eps
168 logical :: drelcmp
169 if (present(eps)) then
170 drelcmp = abs(x - y) .le. eps*abs(y)
171 else
172 drelcmp = abs(x - y) .le. neko_eps*abs(y)
173 end if
174
175 end function drelcmp
176
177
179 pure function qrelcmp(x, y, eps)
180 real(kind=qp), intent(in) :: x
181 real(kind=qp), intent(in) :: y
182 real(kind=qp), intent(in), optional :: eps
183 logical :: qrelcmp
184 if (present(eps)) then
185 qrelcmp = abs(x - y)/abs(y) .lt. eps
186 else
187 qrelcmp = abs(x - y)/abs(y) .lt. neko_eps
188 end if
189
190 end function qrelcmp
191
193 subroutine rzero(a, n)
194 integer, intent(in) :: n
195 real(kind=rp), dimension(n), intent(inout) :: a
196 integer :: i
197
198 do i = 1, n
199 a(i) = 0.0_rp
200 end do
201 end subroutine rzero
202
204 subroutine izero(a, n)
205 integer, intent(in) :: n
206 integer, dimension(n), intent(inout) :: a
207 integer :: i
208
209 do i = 1, n
210 a(i) = 0
211 end do
212 end subroutine izero
213
215 subroutine row_zero(a, m, n, e)
216 integer, intent(in) :: m, n, e
217 real(kind=rp), intent(inout) :: a(m,n)
218 integer :: j
219
220 do j = 1,n
221 a(e,j) = 0.0_rp
222 end do
223 end subroutine row_zero
224
226 subroutine rone(a, n)
227 integer, intent(in) :: n
228 real(kind=rp), dimension(n), intent(inout) :: a
229 integer :: i
230
231 do i = 1, n
232 a(i) = 1.0_rp
233 end do
234 end subroutine rone
235
237 subroutine copy(a, b, n)
238 integer, intent(in) :: n
239 real(kind=rp), dimension(n), intent(in) :: b
240 real(kind=rp), dimension(n), intent(inout) :: a
241 integer :: i
242
243 do i = 1, n
244 a(i) = b(i)
245 end do
246
247 end subroutine copy
248
256 subroutine masked_copy(a, b, mask, n, m)
257 integer, intent(in) :: n, m
258 real(kind=rp), dimension(n), intent(in) :: b
259 real(kind=rp), dimension(n), intent(inout) :: a
260 integer, dimension(0:m) :: mask
261 integer :: i, j
262
263 do i = 1, m
264 j = mask(i)
265 a(j) = b(j)
266 end do
267
268 end subroutine masked_copy
269
278 subroutine masked_red_copy(a, b, mask, n, m)
279 integer, intent(in) :: n, m
280 real(kind=rp), dimension(n), intent(in) :: b
281 real(kind=rp), dimension(m), intent(inout) :: a
282 integer, dimension(0:m) :: mask
283 integer :: i, j
284
285 do i = 1, m
286 j = mask(i)
287 a(i) = b(j)
288 end do
289
290 end subroutine masked_red_copy
291
292
295 subroutine cfill_mask(a, c, size, mask, mask_size)
296 integer, intent(in) :: size, mask_size
297 real(kind=rp), dimension(size), intent(inout) :: a
298 real(kind=rp), intent(in) :: c
299 integer, dimension(mask_size), intent(in) :: mask
300 integer :: i
301
302 do i = 1, mask_size
303 a(mask(i)) = c
304 end do
305
306 end subroutine cfill_mask
307
309 subroutine cmult(a, c, n)
310 integer, intent(in) :: n
311 real(kind=rp), dimension(n), intent(inout) :: a
312 real(kind=rp), intent(in) :: c
313 integer :: i
314
315 do i = 1, n
316 a(i) = c * a(i)
317 end do
318 end subroutine cmult
319
321 subroutine cadd(a, s, n)
322 integer, intent(in) :: n
323 real(kind=rp), dimension(n), intent(inout) :: a
324 real(kind=rp), intent(in) :: s
325 integer :: i
326
327 do i = 1, n
328 a(i) = a(i) + s
329 end do
330 end subroutine cadd
331
333 subroutine cadd2(a, b, s, n)
334 integer, intent(in) :: n
335 real(kind=rp), dimension(n), intent(inout) :: a
336 real(kind=rp), dimension(n), intent(in) :: b
337 real(kind=rp), intent(in) :: s
338 integer :: i
339
340 do i = 1, n
341 a(i) = b(i) + s
342 end do
343 end subroutine cadd2
344
346 subroutine cfill(a, c, n)
347 integer, intent(in) :: n
348 real(kind=rp), dimension(n), intent(inout) :: a
349 real(kind=rp), intent(in) :: c
350 integer :: i
351
352 do i = 1, n
353 a(i) = c
354 end do
355 end subroutine cfill
356
358 function glsum(a, n)
359 integer, intent(in) :: n
360 real(kind=rp), dimension(n) :: a
361 real(kind=rp) :: glsum
362 real(kind=xp) :: tmp
363 integer :: i, ierr
364 tmp = 0.0_rp
365 do i = 1, n
366 tmp = tmp + a(i)
367 end do
368 call mpi_allreduce(mpi_in_place, tmp, 1, &
369 mpi_extra_precision, mpi_sum, neko_comm, ierr)
370 glsum = tmp
371
372 end function glsum
373
375 function glmax(a, n)
376 integer, intent(in) :: n
377 real(kind=rp), dimension(n) :: a
378 real(kind=rp) :: tmp, glmax
379 integer :: i, ierr
380
381 tmp = -huge(0.0_rp)
382 do i = 1, n
383 tmp = max(tmp,a(i))
384 end do
385 call mpi_allreduce(tmp, glmax, 1, &
386 mpi_real_precision, mpi_max, neko_comm, ierr)
387 end function glmax
388
390 function glimax(a, n)
391 integer, intent(in) :: n
392 integer, dimension(n) :: a
393 integer :: tmp, glimax
394 integer :: i, ierr
395
396 tmp = -huge(0)
397 do i = 1, n
398 tmp = max(tmp,a(i))
399 end do
400 call mpi_allreduce(tmp, glimax, 1, &
401 mpi_integer, mpi_max, neko_comm, ierr)
402 end function glimax
403
405 function glmin(a, n)
406 integer, intent(in) :: n
407 real(kind=rp), dimension(n) :: a
408 real(kind=rp) :: tmp, glmin
409 integer :: i, ierr
410
411 tmp = huge(0.0_rp)
412 do i = 1, n
413 tmp = min(tmp,a(i))
414 end do
415 call mpi_allreduce(tmp, glmin, 1, &
416 mpi_real_precision, mpi_min, neko_comm, ierr)
417 end function glmin
418
420 function glimin(a, n)
421 integer, intent(in) :: n
422 integer, dimension(n) :: a
423 integer :: tmp, glimin
424 integer :: i, ierr
425
426 tmp = huge(0)
427 do i = 1, n
428 tmp = min(tmp,a(i))
429 end do
430 call mpi_allreduce(tmp, glimin, 1, &
431 mpi_integer, mpi_min, neko_comm, ierr)
432 end function glimin
433
434
435
436
438 subroutine chsign(a, n)
439 integer, intent(in) :: n
440 real(kind=rp), dimension(n), intent(inout) :: a
441 integer :: i
442
443 do i = 1, n
444 a(i) = -a(i)
445 end do
446
447 end subroutine chsign
448
450 function vlmax(vec,n) result(tmax)
451 integer :: n, i
452 real(kind=rp), intent(in) :: vec(n)
453 real(kind=rp) :: tmax
454 tmax = real(-99d20, rp)
455 do i = 1, n
456 tmax = max(tmax, vec(i))
457 end do
458 end function vlmax
459
461 function vlmin(vec,n) result(tmin)
462 integer, intent(in) :: n
463 real(kind=rp), intent(in) :: vec(n)
464 real(kind=rp) :: tmin
465 integer :: i
466 tmin = real(99.0e20, rp)
467 do i = 1, n
468 tmin = min(tmin, vec(i))
469 end do
470 end function vlmin
471
473 subroutine invcol1(a, n)
474 integer, intent(in) :: n
475 real(kind=rp), dimension(n), intent(inout) :: a
476 integer :: i
477
478 do i = 1, n
479 a(i) = 1.0_xp / real(a(i),xp)
480 end do
481
482 end subroutine invcol1
483
485 subroutine invcol3(a, b, c, n)
486 integer, intent(in) :: n
487 real(kind=rp), dimension(n), intent(inout) :: a
488 real(kind=rp), dimension(n), intent(in) :: b,c
489 integer :: i
490
491 do i = 1, n
492 a(i) = real(b(i),xp) / c(i)
493 end do
494
495 end subroutine invcol3
496
498 subroutine invers2(a, b, n)
499 integer, intent(in) :: n
500 real(kind=rp), dimension(n), intent(inout) :: a
501 real(kind=rp), dimension(n), intent(in) :: b
502 integer :: i
503
504 do i = 1, n
505 a(i) = 1.0_xp / real(b(i),xp)
506 end do
507
508 end subroutine invers2
509
512 subroutine vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
513 integer, intent(in) :: n
514 real(kind=rp), dimension(n), intent(in) :: v1, v2, v3
515 real(kind=rp), dimension(n), intent(in) :: w1, w2, w3
516 real(kind=rp), dimension(n), intent(out) :: u1, u2, u3
517 integer :: i
518
519 do i = 1, n
520 u1(i) = v2(i)*w3(i) - v3(i)*w2(i)
521 u2(i) = v3(i)*w1(i) - v1(i)*w3(i)
522 u3(i) = v1(i)*w2(i) - v2(i)*w1(i)
523 end do
524
525 end subroutine vcross
526
529 subroutine vdot2(dot, u1, u2, v1, v2, n)
530 integer, intent(in) :: n
531 real(kind=rp), dimension(n), intent(in) :: u1, u2
532 real(kind=rp), dimension(n), intent(in) :: v1, v2
533 real(kind=rp), dimension(n), intent(out) :: dot
534 integer :: i
535 do i = 1, n
536 dot(i) = u1(i)*v1(i) + u2(i)*v2(i)
537 end do
538
539 end subroutine vdot2
540
543 subroutine vdot3(dot, u1, u2, u3, v1, v2, v3, n)
544 integer, intent(in) :: n
545 real(kind=rp), dimension(n), intent(in) :: u1, u2, u3
546 real(kind=rp), dimension(n), intent(in) :: v1, v2, v3
547 real(kind=rp), dimension(n), intent(out) :: dot
548 integer :: i
549
550 do i = 1, n
551 dot(i) = u1(i)*v1(i) + u2(i)*v2(i) + u3(i)*v3(i)
552 end do
553
554 end subroutine vdot3
555
557 function vlsc3(u, v, w, n) result(s)
558 integer, intent(in) :: n
559 real(kind=rp), dimension(n), intent(in) :: u, v, w
560 real(kind=rp) :: s
561 integer :: i
562
563 s = 0.0_rp
564 do i = 1, n
565 s = s + u(i)*v(i)*w(i)
566 end do
567
568 end function vlsc3
569
571 function vlsc2(u, v, n) result(s)
572 integer, intent(in) :: n
573 real(kind=rp), dimension(n), intent(in) :: u, v
574 real(kind=rp) :: s
575 integer :: i
576
577 s = 0.0_rp
578 do i = 1, n
579 s = s + u(i)*v(i)
580 end do
581
582 end function vlsc2
583
585 subroutine add2(a, b, n)
586 integer, intent(in) :: n
587 real(kind=rp), dimension(n), intent(inout) :: a
588 real(kind=rp), dimension(n), intent(in) :: b
589 integer :: i
590
591 do i = 1, n
592 a(i) = a(i) + b(i)
593 end do
594
595 end subroutine add2
596
598 subroutine add3(a, b, c, n)
599 integer, intent(in) :: n
600 real(kind=rp), dimension(n), intent(inout) :: a
601 real(kind=rp), dimension(n), intent(in) :: b
602 real(kind=rp), dimension(n), intent(in) :: c
603 integer :: i
604
605 do i = 1, n
606 a(i) = b(i) + c(i)
607 end do
608
609 end subroutine add3
610
612 subroutine add4(a, b, c, d, n)
613 integer, intent(in) :: n
614 real(kind=rp), dimension(n), intent(out) :: a
615 real(kind=rp), dimension(n), intent(in) :: d
616 real(kind=rp), dimension(n), intent(in) :: c
617 real(kind=rp), dimension(n), intent(in) :: b
618 integer :: i
619
620 do i = 1, n
621 a(i) = b(i) + c(i) + d(i)
622 end do
623
624 end subroutine add4
625
627 subroutine sub2(a, b, n)
628 integer, intent(in) :: n
629 real(kind=rp), dimension(n), intent(inout) :: a
630 real(kind=rp), dimension(n), intent(inout) :: b
631 integer :: i
632
633 do i = 1, n
634 a(i) = a(i) - b(i)
635 end do
636
637 end subroutine sub2
638
640 subroutine sub3(a, b, c, n)
641 integer, intent(in) :: n
642 real(kind=rp), dimension(n), intent(inout) :: a
643 real(kind=rp), dimension(n), intent(in) :: b
644 real(kind=rp), dimension(n), intent(in) :: c
645 integer :: i
646
647 do i = 1, n
648 a(i) = b(i) - c(i)
649 end do
650
651 end subroutine sub3
652
653
656 subroutine add2s1(a, b, c1, n)
657 integer, intent(in) :: n
658 real(kind=rp), dimension(n), intent(inout) :: a
659 real(kind=rp), dimension(n), intent(inout) :: b
660 real(kind=rp), intent(in) :: c1
661 integer :: i
662
663 do i = 1, n
664 a(i) = c1 * a(i) + b(i)
665 end do
666
667 end subroutine add2s1
668
671 subroutine add2s2(a, b, c1, n)
672 integer, intent(in) :: n
673 real(kind=rp), dimension(n), intent(inout) :: a
674 real(kind=rp), dimension(n), intent(inout) :: b
675 real(kind=rp), intent(in) :: c1
676 integer :: i
677
678 do i = 1, n
679 a(i) = a(i) + c1 * b(i)
680 end do
681
682 end subroutine add2s2
683
685 subroutine addsqr2s2(a, b, c1, n)
686 integer, intent(in) :: n
687 real(kind=rp), dimension(n), intent(inout) :: a
688 real(kind=rp), dimension(n), intent(in) :: b
689 real(kind=rp), intent(in) :: c1
690 integer :: i
691
692 do i = 1,n
693 a(i) = a(i) + c1 * ( b(i) * b(i) )
694 end do
695
696 end subroutine addsqr2s2
697
699 subroutine cmult2(a, b, c, n)
700 integer, intent(in) :: n
701 real(kind=rp), dimension(n), intent(inout) :: a
702 real(kind=rp), dimension(n), intent(in) :: b
703 real(kind=rp), intent(in) :: c
704 integer :: i
705
706 do i = 1, n
707 a(i) = c * b(i)
708 end do
709
710 end subroutine cmult2
711
713 subroutine invcol2(a, b, n)
714 integer, intent(in) :: n
715 real(kind=rp), dimension(n), intent(inout) :: a
716 real(kind=rp), dimension(n), intent(in) :: b
717 integer :: i
718
719 do i = 1, n
720 a(i) = real(a(i),xp) /b(i)
721 end do
722
723 end subroutine invcol2
724
725
727 subroutine col2(a, b, n)
728 integer, intent(in) :: n
729 real(kind=rp), dimension(n), intent(inout) :: a
730 real(kind=rp), dimension(n), intent(in) :: b
731 integer :: i
732
733 do i = 1, n
734 a(i) = a(i) * b(i)
735 end do
736
737 end subroutine col2
738
740 subroutine col3(a, b, c, n)
741 integer, intent(in) :: n
742 real(kind=rp), dimension(n), intent(inout) :: a
743 real(kind=rp), dimension(n), intent(in) :: b
744 real(kind=rp), dimension(n), intent(in) :: c
745 integer :: i
746
747 do i = 1, n
748 a(i) = b(i) * c(i)
749 end do
750
751 end subroutine col3
752
754 subroutine subcol3(a, b, c, n)
755 integer, intent(in) :: n
756 real(kind=rp), dimension(n), intent(inout) :: a
757 real(kind=rp), dimension(n), intent(in) :: b
758 real(kind=rp), dimension(n), intent(in) :: c
759 integer :: i
760
761 do i = 1,n
762 a(i) = a(i) - b(i) * c(i)
763 end do
764
765 end subroutine subcol3
766
768 subroutine add3s2(a, b, c, c1, c2 ,n)
769 integer, intent(in) :: n
770 real(kind=rp), dimension(n), intent(inout) :: a
771 real(kind=rp), dimension(n), intent(in) :: b
772 real(kind=rp), dimension(n), intent(in) :: c
773 real(kind=rp), intent(in) :: c1, c2
774 integer :: i
775
776 do i = 1,n
777 a(i) = c1 * b(i) + c2 * c(i)
778 end do
779
780 end subroutine add3s2
781
782
784 subroutine subcol4(a, b, c, d, n)
785 integer, intent(in) :: n
786 real(kind=rp), dimension(n), intent(inout) :: a
787 real(kind=rp), dimension(n), intent(in) :: b
788 real(kind=rp), dimension(n), intent(in) :: c
789 real(kind=rp), dimension(n), intent(in) :: d
790 integer :: i
791
792 do i = 1,n
793 a(i) = a(i) - b(i) * c(i) * d(i)
794 end do
795
796 end subroutine subcol4
797
799 subroutine addcol3(a, b, c, n)
800 integer, intent(in) :: n
801 real(kind=rp), dimension(n), intent(inout) :: a
802 real(kind=rp), dimension(n), intent(in) :: b
803 real(kind=rp), dimension(n), intent(in) :: c
804 integer :: i
805
806 do i = 1,n
807 a(i) = a(i) + b(i) * c(i)
808 end do
809
810 end subroutine addcol3
811
813 subroutine addcol4(a, b, c, d, n)
814 integer, intent(in) :: n
815 real(kind=rp), dimension(n), intent(inout) :: a
816 real(kind=rp), dimension(n), intent(in) :: b
817 real(kind=rp), dimension(n), intent(in) :: c
818 real(kind=rp), dimension(n), intent(in) :: d
819 integer :: i
820
821 do i = 1,n
822 a(i) = a(i) + b(i) * c(i) * d(i)
823 end do
824
825 end subroutine addcol4
826
828 subroutine ascol5(a, b, c, d, e, n)
829 integer, intent(in) :: n
830 real(kind=rp), dimension(n), intent(inout) :: a
831 real(kind=rp), dimension(n), intent(in) :: b
832 real(kind=rp), dimension(n), intent(in) :: c
833 real(kind=rp), dimension(n), intent(in) :: d
834 real(kind=rp), dimension(n), intent(in) :: e
835 integer :: i
836
837 do i = 1,n
838 a(i) = b(i)*c(i)-d(i)*e(i)
839 end do
840
841 end subroutine ascol5
842
844 subroutine p_update(a, b, c, c1, c2, n)
845 integer, intent(in) :: n
846 real(kind=rp), dimension(n), intent(inout) :: a
847 real(kind=rp), dimension(n), intent(in) :: b
848 real(kind=rp), dimension(n), intent(in) :: c
849 real(kind=rp), intent(in) :: c1, c2
850 integer :: i
851
852 do i = 1,n
853 a(i) = b(i) + c1*(a(i)-c2*c(i))
854 end do
855
856 end subroutine p_update
857
859 subroutine x_update(a, b, c, c1, c2, n)
860 integer, intent(in) :: n
861 real(kind=rp), dimension(n), intent(inout) :: a
862 real(kind=rp), dimension(n), intent(in) :: b
863 real(kind=rp), dimension(n), intent(in) :: c
864 real(kind=rp), intent(in) :: c1, c2
865 integer :: i
866
867 do i = 1,n
868 a(i) = a(i) + c1*b(i)+c2*c(i)
869 end do
870
871 end subroutine x_update
872
874 function glsc2(a, b, n)
875 integer, intent(in) :: n
876 real(kind=rp), dimension(n), intent(in) :: a
877 real(kind=rp), dimension(n), intent(in) :: b
878 real(kind=rp) :: glsc2
879 real(kind=xp) :: tmp
880 integer :: i, ierr
881
882 tmp = 0.0_xp
883 do i = 1, n
884 tmp = tmp + a(i) * b(i)
885 end do
886
887 call mpi_allreduce(mpi_in_place, tmp, 1, &
888 mpi_extra_precision, mpi_sum, neko_comm, ierr)
889 glsc2 = tmp
890 end function glsc2
891
893 function glsc3(a, b, c, n)
894 integer, intent(in) :: n
895 real(kind=rp), dimension(n), intent(in) :: a
896 real(kind=rp), dimension(n), intent(in) :: b
897 real(kind=rp), dimension(n), intent(in) :: c
898 real(kind=rp) :: glsc3
899 real(kind=xp) :: tmp
900 integer :: i, ierr
901
902 tmp = 0.0_xp
903 do i = 1, n
904 tmp = tmp + a(i) * b(i) * c(i)
905 end do
906
907 call mpi_allreduce(mpi_in_place, tmp, 1, &
908 mpi_extra_precision, mpi_sum, neko_comm, ierr)
909 glsc3 = tmp
910
911 end function glsc3
912 function glsc4(a, b, c, d, n)
913 integer, intent(in) :: n
914 real(kind=rp), dimension(n), intent(in) :: a
915 real(kind=rp), dimension(n), intent(in) :: b
916 real(kind=rp), dimension(n), intent(in) :: c
917 real(kind=rp), dimension(n), intent(in) :: d
918 real(kind=rp) :: glsc4
919 real(kind=xp) :: tmp
920 integer :: i, ierr
921
922 tmp = 0.0_xp
923 do i = 1, n
924 tmp = tmp + a(i) * b(i) * c(i) * d(i)
925 end do
926
927 call mpi_allreduce(mpi_in_place, tmp, 1, &
928 mpi_extra_precision, mpi_sum, neko_comm, ierr)
929 glsc4 = tmp
930
931 end function glsc4
932
938 subroutine sortrp(a, ind, n)
939 integer, intent(in) :: n
940 real(kind=rp), intent(inout) :: a(n)
941 integer, intent(out) :: ind(n)
942 real(kind=rp) :: aa
943 integer :: j, ir, i, ii, l
944
945 do j = 1, n
946 ind(j) = j
947 end do
948
949 if (n .le. 1) return
950
951
952 l = n/2+1
953 ir = n
954 do while (.true.)
955 if (l .gt. 1) then
956 l = l-1
957 aa = a(l)
958 ii = ind(l)
959 else
960 aa = a(ir)
961 ii = ind(ir)
962 a(ir) = a(1)
963 ind(ir) = ind(1)
964 ir = ir - 1
965 if (ir .eq. 1) then
966 a(1) = aa
967 ind(1) = ii
968 return
969 end if
970 end if
971 i = l
972 j = l+l
973 do while (j .le. ir)
974 if (j .lt. ir) then
975 if ( a(j) .lt. a(j+1) ) j = j + 1
976 end if
977 if (aa .lt. a(j)) then
978 a(i) = a(j)
979 ind(i) = ind(j)
980 i = j
981 j = j+j
982 else
983 j = ir+1
984 end if
985 end do
986 a(i) = aa
987 ind(i) = ii
988 end do
989 end subroutine sortrp
990
996 subroutine sorti4(a, ind, n)
997 integer, intent(in) :: n
998 integer(i4), intent(inout) :: a(n)
999 integer, intent(out) :: ind(n)
1000 integer(i4) :: aa
1001 integer :: j, ir, i, ii, l
1002
1003 do j = 1, n
1004 ind(j) = j
1005 end do
1006
1007 if (n .le. 1) return
1008
1009 l = n/2+1
1010 ir = n
1011 do while (.true.)
1012 if (l .gt. 1) then
1013 l = l - 1
1014 aa = a(l)
1015 ii = ind(l)
1016 else
1017 aa = a(ir)
1018 ii = ind(ir)
1019 a(ir) = a( 1)
1020 ind(ir) = ind( 1)
1021 ir = ir - 1
1022 if (ir .eq. 1) then
1023 a(1) = aa
1024 ind(1) = ii
1025 return
1026 end if
1027 end if
1028 i = l
1029 j = l + l
1030 do while (j .le. ir)
1031 if (j .lt. ir) then
1032 if ( a(j) .lt. a(j + 1) ) j = j + 1
1033 end if
1034 if (aa .lt. a(j)) then
1035 a(i) = a(j)
1036 ind(i) = ind(j)
1037 i = j
1038 j = j + j
1039 else
1040 j = ir + 1
1041 end if
1042 end do
1043 a(i) = aa
1044 ind(i) = ii
1045 end do
1046 end subroutine sorti4
1047
1052 subroutine swapdp(b, ind, n)
1053 integer, intent(in) :: n
1054 real(kind=rp), intent(inout) :: b(n)
1055 integer, intent(in) :: ind(n)
1056 real(kind=rp) :: temp(n)
1057 integer :: i, jj
1058
1059 do i = 1, n
1060 temp(i) = b(i)
1061 end do
1062 do i = 1, n
1063 jj = ind(i)
1064 b(i) = temp(jj)
1065 end do
1066 end subroutine swapdp
1067
1072 subroutine swapi4(b, ind, n)
1073 integer, intent(in) :: n
1074 integer(i4), intent(inout) :: b(n)
1075 integer, intent(in) :: ind(n)
1076 integer(i4) :: temp(n)
1077 integer :: i, jj
1078
1079 do i = 1, n
1080 temp(i) = b(i)
1081 end do
1082 do i = 1, n
1083 jj = ind(i)
1084 b(i) = temp(jj)
1085 end do
1086 end subroutine swapi4
1087
1092 subroutine reorddp(b, ind, n)
1093 integer, intent(in) :: n
1094 real(kind=rp), intent(inout) :: b(n)
1095 integer, intent(in) :: ind(n)
1096 real(kind=rp) :: temp(n)
1097 integer :: i, jj
1098
1099 do i = 1, n
1100 temp(i) = b(i)
1101 end do
1102 do i = 1, n
1103 jj = ind(i)
1104 b(jj) = temp(i)
1105 end do
1106 end subroutine reorddp
1107
1112 subroutine reordi4(b, ind, n)
1113 integer, intent(in) :: n
1114 integer(i4), intent(inout) :: b(n)
1115 integer, intent(in) :: ind(n)
1116 integer(i4) :: temp(n)
1117 integer :: i, jj
1118
1119 do i = 1, n
1120 temp(i) = b(i)
1121 end do
1122 do i = 1, n
1123 jj = ind(i)
1124 b(jj) = temp(i)
1125 end do
1126 end subroutine reordi4
1127
1132 subroutine flipvdp(b, ind, n)
1133 integer, intent(in) :: n
1134 real(kind=rp), intent(inout) :: b(n)
1135 integer, intent(inout) :: ind(n)
1136 real(kind=rp) :: temp(n)
1137 integer :: tempind(n)
1138 integer :: i, jj
1139
1140 do i = 1, n
1141 jj = n+1-i
1142 temp(jj) = b(i)
1143 tempind(jj) = ind(i)
1144 end do
1145 do i = 1,n
1146 b(i) = temp(i)
1147 ind(i) = tempind(i)
1148 end do
1149 end subroutine flipvdp
1150
1155 subroutine flipvi4(b, ind, n)
1156 integer, intent(in) :: n
1157 integer(i4), intent(inout) :: b(n)
1158 integer, intent(inout) :: ind(n)
1159 integer(i4) :: temp(n)
1160 integer :: tempind(n)
1161 integer :: i, jj
1162
1163 do i = 1, n
1164 jj = n+1-i
1165 temp(jj) = b(i)
1166 tempind(jj) = ind(i)
1167 end do
1168 do i = 1,n
1169 b(i) = temp(i)
1170 ind(i) = tempind(i)
1171 end do
1172 end subroutine flipvi4
1173
1177 subroutine absval(a, n)
1178 integer, intent(in) :: n
1179 real(kind=rp), dimension(n), intent(inout) :: a
1180 integer :: i
1181 do i = 1, n
1182 a(i) = abs(a(i))
1183 end do
1184 end subroutine absval
1185
1186 ! ========================================================================== !
1187 ! Point-wise operations
1188
1190 subroutine pwmax_vec2(a, b, n)
1191 integer, intent(in) :: n
1192 real(kind=rp), dimension(n), intent(inout) :: a
1193 real(kind=rp), dimension(n), intent(in) :: b
1194 integer :: i
1195
1196 do i = 1, n
1197 a(i) = max(a(i), b(i))
1198 end do
1199 end subroutine pwmax_vec2
1200
1202 subroutine pwmax_vec3(a, b, c, n)
1203 integer, intent(in) :: n
1204 real(kind=rp), dimension(n), intent(inout) :: a
1205 real(kind=rp), dimension(n), intent(in) :: b, c
1206 integer :: i
1207
1208 do i = 1, n
1209 a(i) = max(b(i), c(i))
1210 end do
1211 end subroutine pwmax_vec3
1212
1214 subroutine pwmax_scal2(a, b, n)
1215 integer, intent(in) :: n
1216 real(kind=rp), dimension(n), intent(inout) :: a
1217 real(kind=rp), intent(in) :: b
1218 integer :: i
1219
1220 do i = 1, n
1221 a(i) = max(a(i), b)
1222 end do
1223 end subroutine pwmax_scal2
1224
1226 subroutine pwmax_scal3(a, b, c, n)
1227 integer, intent(in) :: n
1228 real(kind=rp), dimension(n), intent(inout) :: a
1229 real(kind=rp), dimension(n), intent(in) :: b
1230 real(kind=rp), intent(in) :: c
1231 integer :: i
1232
1233 do i = 1, n
1234 a(i) = max(b(i), c)
1235 end do
1236 end subroutine pwmax_scal3
1237
1239 subroutine pwmin_vec2(a, b, n)
1240 integer, intent(in) :: n
1241 real(kind=rp), dimension(n), intent(inout) :: a
1242 real(kind=rp), dimension(n), intent(in) :: b
1243 integer :: i
1244
1245 do i = 1, n
1246 a(i) = min(a(i), b(i))
1247 end do
1248 end subroutine pwmin_vec2
1249
1251 subroutine pwmin_vec3(a, b, c, n)
1252 integer, intent(in) :: n
1253 real(kind=rp), dimension(n), intent(inout) :: a
1254 real(kind=rp), dimension(n), intent(in) :: b, c
1255 integer :: i
1256
1257 do i = 1, n
1258 a(i) = min(b(i), c(i))
1259 end do
1260 end subroutine pwmin_vec3
1261
1263 subroutine pwmin_sca2(a, b, n)
1264 integer, intent(in) :: n
1265 real(kind=rp), dimension(n), intent(inout) :: a
1266 real(kind=rp), intent(in) :: b
1267 integer :: i
1268
1269 do i = 1, n
1270 a(i) = min(a(i), b)
1271 end do
1272 end subroutine pwmin_sca2
1273
1275 subroutine pwmin_sca3(a, b, c, n)
1276 integer, intent(in) :: n
1277 real(kind=rp), dimension(n), intent(inout) :: a
1278 real(kind=rp), dimension(n), intent(in) :: b
1279 real(kind=rp), intent(in) :: c
1280 integer :: i
1281
1282 do i = 1, n
1283 a(i) = min(b(i), c)
1284 end do
1285 end subroutine pwmin_sca3
1286
1287end module math
double real
Definition comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition comm.F90:16
type(mpi_datatype) mpi_real_precision
MPI type for working precision of REAL types.
Definition comm.F90:23
type(mpi_datatype) mpi_extra_precision
Definition comm.F90:24
Definition math.f90:60
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition math.f90:310
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
Definition math.f90:700
subroutine, public row_zero(a, m, n, e)
Sets row e to 0 in matrix a.
Definition math.f90:216
subroutine, public invcol2(a, b, n)
Vector division .
Definition math.f90:714
real(kind=rp) function, public vlsc2(u, v, n)
Compute multiplication sum .
Definition math.f90:572
subroutine pwmin_sca3(a, b, c, n)
Point-wise minimum of scalar and vector .
Definition math.f90:1276
real(kind=rp), parameter, public pi
Definition math.f90:75
subroutine pwmax_scal2(a, b, n)
Point-wise maximum of scalar and vector .
Definition math.f90:1215
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition math.f90:894
pure logical function qabscmp(x, y)
Return double precision absolute comparison .
Definition math.f90:141
subroutine, public ascol5(a, b, c, d, e, n)
Returns .
Definition math.f90:829
subroutine, public invers2(a, b, n)
Compute inverted vector .
Definition math.f90:499
subroutine, public cadd2(a, b, s, n)
Add a scalar to vector .
Definition math.f90:334
subroutine pwmin_sca2(a, b, n)
Point-wise minimum of scalar and vector .
Definition math.f90:1264
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:322
subroutine reorddp(b, ind, n)
reorder double precision array - inverse of swap
Definition math.f90:1093
subroutine, public addsqr2s2(a, b, c1, n)
Returns .
Definition math.f90:686
real(kind=rp) function, public glsc4(a, b, c, d, n)
Definition math.f90:913
subroutine, public masked_red_copy(a, b, mask, n, m)
Copy a masked vector to reduced contigous vector .
Definition math.f90:279
subroutine pwmax_vec3(a, b, c, n)
Point-wise maximum of two vectors .
Definition math.f90:1203
subroutine swapdp(b, ind, n)
sort double precision array acording to ind vector
Definition math.f90:1053
subroutine flipvi4(b, ind, n)
Flip single integer vector b and ind.
Definition math.f90:1156
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
Definition math.f90:657
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition math.f90:875
subroutine, public subcol3(a, b, c, n)
Returns .
Definition math.f90:755
subroutine, public rone(a, n)
Set all elements to one.
Definition math.f90:227
subroutine flipvdp(b, ind, n)
Flip double precision vector b and ind.
Definition math.f90:1133
subroutine, public x_update(a, b, c, c1, c2, n)
Returns .
Definition math.f90:860
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:599
subroutine swapi4(b, ind, n)
sort single integer array acording to ind vector
Definition math.f90:1073
integer function, public glimin(a, n)
Min of an integer vector of length n.
Definition math.f90:421
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
Definition math.f90:359
subroutine, public sub3(a, b, c, n)
Vector subtraction .
Definition math.f90:641
subroutine, public addcol4(a, b, c, d, n)
Returns .
Definition math.f90:814
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:586
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition math.f90:347
subroutine, public absval(a, n)
Take the absolute value of an array.
Definition math.f90:1178
subroutine, public invcol3(a, b, c, n)
Invert a vector .
Definition math.f90:486
pure logical function sabscmp(x, y)
Return single precision absolute comparison .
Definition math.f90:121
subroutine, public add3s2(a, b, c, c1, c2, n)
Returns .
Definition math.f90:769
subroutine, public subcol4(a, b, c, d, n)
Returns .
Definition math.f90:785
subroutine sorti4(a, ind, n)
Heap Sort for single integer arrays.
Definition math.f90:997
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:800
pure logical function dabscmp(x, y)
Return double precision absolute comparison .
Definition math.f90:131
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:474
subroutine, public masked_copy(a, b, mask, n, m)
Copy a masked vector .
Definition math.f90:257
real(kind=rp), parameter, public neko_m_ln2
Definition math.f90:72
subroutine, public chsign(a, n)
Change sign of vector .
Definition math.f90:439
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:728
subroutine, public izero(a, n)
Zero an integer vector.
Definition math.f90:205
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition math.f90:376
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
subroutine pwmax_scal3(a, b, c, n)
Point-wise maximum of scalar and vector .
Definition math.f90:1227
subroutine, public add4(a, b, c, d, n)
Vector addition .
Definition math.f90:613
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:741
subroutine pwmin_vec2(a, b, n)
Point-wise minimum of two vectors .
Definition math.f90:1240
subroutine pwmin_vec3(a, b, c, n)
Point-wise minimum of two vectors .
Definition math.f90:1252
subroutine pwmax_vec2(a, b, n)
Point-wise maximum of two vectors .
Definition math.f90:1191
pure logical function drelcmp(x, y, eps)
Return double precision relative comparison .
Definition math.f90:165
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:69
subroutine, public vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Compute a dot product (3-d version) assuming vector components etc.
Definition math.f90:544
pure logical function qrelcmp(x, y, eps)
Return quad precision relative comparison .
Definition math.f90:180
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:194
subroutine, public vdot2(dot, u1, u2, v1, v2, n)
Compute a dot product (2-d version) assuming vector components etc.
Definition math.f90:530
real(kind=rp) function, public vlmin(vec, n)
minimun value of a vector of length n
Definition math.f90:462
real(kind=rp) function, public vlmax(vec, n)
maximum value of a vector of length n
Definition math.f90:451
integer function, public glimax(a, n)
Max of an integer vector of length n.
Definition math.f90:391
subroutine sortrp(a, ind, n)
Heap Sort for double precision arrays.
Definition math.f90:939
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
subroutine, public cfill_mask(a, c, size, mask, mask_size)
Fill a constant to a masked vector. .
Definition math.f90:296
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition math.f90:672
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition math.f90:406
subroutine, public vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Compute a cross product assuming vector components etc.
Definition math.f90:513
pure logical function srelcmp(x, y, eps)
Return single precision relative comparison .
Definition math.f90:151
real(kind=rp) function, public vlsc3(u, v, w, n)
Compute multiplication sum .
Definition math.f90:558
subroutine reordi4(b, ind, n)
reorder single integer array - inverse of swap
Definition math.f90:1113
subroutine, public p_update(a, b, c, c1, c2, n)
Returns .
Definition math.f90:845
integer, parameter, public qp
Definition num_types.f90:10
integer, parameter, public i4
Definition num_types.f90:6
integer, parameter, public xp
Definition num_types.f90:14
integer, parameter, public dp
Definition num_types.f90:9
integer, parameter, public sp
Definition num_types.f90:8
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
#define max(a, b)
Definition tensor.cu:40