Neko 0.9.1
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
62 use utils, only: neko_error
64 use mpi_f08, only: mpi_min, mpi_max, mpi_sum, mpi_in_place, mpi_integer, &
65 mpi_allreduce
66 implicit none
67 private
68
70 real(kind=rp), public, parameter :: neko_eps = epsilon(1.0_rp)
71
73 real(kind=rp), public, parameter :: neko_m_ln2 = log(2.0_rp)
74
76 real(kind=rp), public, parameter :: pi = 4._rp*atan(1._rp)
77
78 interface abscmp
79 module procedure sabscmp, dabscmp, qabscmp
80 end interface abscmp
81
82 interface sort
83 module procedure sortrp, sorti4
84 end interface sort
85
86 interface swap
87 module procedure swapdp, swapi4
88 end interface swap
89
90 interface reord
91 module procedure reorddp, reordi4
92 end interface reord
93
94 interface flipv
95 module procedure flipvdp, flipvi4
96 end interface flipv
97
98 interface relcmp
99 module procedure srelcmp, drelcmp, qrelcmp
100 end interface relcmp
101
102 interface pwmax
103 module procedure pwmax_vec2, pwmax_vec3, pwmax_scal2, pwmax_scal3
104 end interface pwmax
105
106 interface pwmin
107 module procedure pwmin_vec2, pwmin_vec3, pwmin_sca2, pwmin_sca3
108 end interface pwmin
109
110 public :: abscmp, rzero, izero, row_zero, rone, copy, cmult, cadd, cfill, &
117
118contains
119
121 pure function sabscmp(x, y)
122 real(kind=sp), intent(in) :: x
123 real(kind=sp), intent(in) :: y
124 logical :: sabscmp
125
126 sabscmp = abs(x - y) .lt. neko_eps
127
128 end function sabscmp
129
131 pure function dabscmp(x, y)
132 real(kind=dp), intent(in) :: x
133 real(kind=dp), intent(in) :: y
134 logical :: dabscmp
135
136 dabscmp = abs(x - y) .lt. neko_eps
137
138 end function dabscmp
139
141 pure function qabscmp(x, y)
142 real(kind=qp), intent(in) :: x
143 real(kind=qp), intent(in) :: y
144 logical :: qabscmp
145
146 qabscmp = abs(x - y) .lt. neko_eps
147
148 end function qabscmp
149
151 pure function srelcmp(x, y, eps)
152 real(kind=sp), intent(in) :: x
153 real(kind=sp), intent(in) :: y
154 real(kind=sp), intent(in), optional :: eps
155 logical :: srelcmp
156 if (present(eps)) then
157 srelcmp = abs(x - y) .le. eps*abs(y)
158 else
159 srelcmp = abs(x - y) .le. neko_eps*abs(y)
160 end if
161
162 end function srelcmp
163
165 pure function drelcmp(x, y, eps)
166 real(kind=dp), intent(in) :: x
167 real(kind=dp), intent(in) :: y
168 real(kind=dp), intent(in), optional :: eps
169 logical :: drelcmp
170 if (present(eps)) then
171 drelcmp = abs(x - y) .le. eps*abs(y)
172 else
173 drelcmp = abs(x - y) .le. neko_eps*abs(y)
174 end if
175
176 end function drelcmp
177
178
180 pure function qrelcmp(x, y, eps)
181 real(kind=qp), intent(in) :: x
182 real(kind=qp), intent(in) :: y
183 real(kind=qp), intent(in), optional :: eps
184 logical :: qrelcmp
185 if (present(eps)) then
186 qrelcmp = abs(x - y)/abs(y) .lt. eps
187 else
188 qrelcmp = abs(x - y)/abs(y) .lt. neko_eps
189 end if
190
191 end function qrelcmp
192
194 subroutine rzero(a, n)
195 integer, intent(in) :: n
196 real(kind=rp), dimension(n), intent(inout) :: a
197 integer :: i
198
199 do i = 1, n
200 a(i) = 0.0_rp
201 end do
202 end subroutine rzero
203
205 subroutine izero(a, n)
206 integer, intent(in) :: n
207 integer, dimension(n), intent(inout) :: a
208 integer :: i
209
210 do i = 1, n
211 a(i) = 0
212 end do
213 end subroutine izero
214
216 subroutine row_zero(a, m, n, e)
217 integer, intent(in) :: m, n, e
218 real(kind=rp), intent(inout) :: a(m,n)
219 integer :: j
220
221 do j = 1,n
222 a(e,j) = 0.0_rp
223 end do
224 end subroutine row_zero
225
227 subroutine rone(a, n)
228 integer, intent(in) :: n
229 real(kind=rp), dimension(n), intent(inout) :: a
230 integer :: i
231
232 do i = 1, n
233 a(i) = 1.0_rp
234 end do
235 end subroutine rone
236
238 subroutine copy(a, b, n)
239 integer, intent(in) :: n
240 real(kind=rp), dimension(n), intent(in) :: b
241 real(kind=rp), dimension(n), intent(inout) :: a
242 integer :: i
243
244 do i = 1, n
245 a(i) = b(i)
246 end do
247
248 end subroutine copy
249
257 subroutine masked_copy(a, b, mask, n, m)
258 integer, intent(in) :: n, m
259 real(kind=rp), dimension(n), intent(in) :: b
260 real(kind=rp), dimension(n), intent(inout) :: a
261 integer, dimension(0:m) :: mask
262 integer :: i, j
263
264 do i = 1, m
265 j = mask(i)
266 a(j) = b(j)
267 end do
268
269 end subroutine masked_copy
270
279 subroutine masked_red_copy(a, b, mask, n, m)
280 integer, intent(in) :: n, m
281 real(kind=rp), dimension(n), intent(in) :: b
282 real(kind=rp), dimension(m), intent(inout) :: a
283 integer, dimension(0:m) :: mask
284 integer :: i, j
285
286 do i = 1, m
287 j = mask(i)
288 a(i) = b(j)
289 end do
290
291 end subroutine masked_red_copy
292
293
296 subroutine cfill_mask(a, c, size, mask, mask_size)
297 integer, intent(in) :: size, mask_size
298 real(kind=rp), dimension(size), intent(inout) :: a
299 real(kind=rp), intent(in) :: c
300 integer, dimension(mask_size), intent(in) :: mask
301 integer :: i
302
303 do i = 1, mask_size
304 a(mask(i)) = c
305 end do
306
307 end subroutine cfill_mask
308
310 subroutine cmult(a, c, n)
311 integer, intent(in) :: n
312 real(kind=rp), dimension(n), intent(inout) :: a
313 real(kind=rp), intent(in) :: c
314 integer :: i
315
316 do i = 1, n
317 a(i) = c * a(i)
318 end do
319 end subroutine cmult
320
322 subroutine cadd(a, s, n)
323 integer, intent(in) :: n
324 real(kind=rp), dimension(n), intent(inout) :: a
325 real(kind=rp), intent(in) :: s
326 integer :: i
327
328 do i = 1, n
329 a(i) = a(i) + s
330 end do
331 end subroutine cadd
332
334 subroutine cadd2(a, b, s, n)
335 integer, intent(in) :: n
336 real(kind=rp), dimension(n), intent(inout) :: a
337 real(kind=rp), dimension(n), intent(in) :: b
338 real(kind=rp), intent(in) :: s
339 integer :: i
340
341 do i = 1, n
342 a(i) = b(i) + s
343 end do
344 end subroutine cadd2
345
347 subroutine cfill(a, c, n)
348 integer, intent(in) :: n
349 real(kind=rp), dimension(n), intent(inout) :: a
350 real(kind=rp), intent(in) :: c
351 integer :: i
352
353 do i = 1, n
354 a(i) = c
355 end do
356 end subroutine cfill
357
359 function glsum(a, n)
360 integer, intent(in) :: n
361 real(kind=rp), dimension(n) :: a
362 real(kind=rp) :: glsum
363 real(kind=xp) :: tmp
364 integer :: i, ierr
365 tmp = 0.0_rp
366 do i = 1, n
367 tmp = tmp + a(i)
368 end do
369 call mpi_allreduce(mpi_in_place, tmp, 1, &
370 mpi_extra_precision, mpi_sum, neko_comm, ierr)
371 glsum = tmp
372
373 end function glsum
374
376 function glmax(a, n)
377 integer, intent(in) :: n
378 real(kind=rp), dimension(n) :: a
379 real(kind=rp) :: tmp, glmax
380 integer :: i, ierr
381
382 tmp = -huge(0.0_rp)
383 do i = 1, n
384 tmp = max(tmp,a(i))
385 end do
386 call mpi_allreduce(tmp, glmax, 1, &
387 mpi_real_precision, mpi_max, neko_comm, ierr)
388 end function glmax
389
391 function glimax(a, n)
392 integer, intent(in) :: n
393 integer, dimension(n) :: a
394 integer :: tmp, glimax
395 integer :: i, ierr
396
397 tmp = -huge(0)
398 do i = 1, n
399 tmp = max(tmp,a(i))
400 end do
401 call mpi_allreduce(tmp, glimax, 1, &
402 mpi_integer, mpi_max, neko_comm, ierr)
403 end function glimax
404
406 function glmin(a, n)
407 integer, intent(in) :: n
408 real(kind=rp), dimension(n) :: a
409 real(kind=rp) :: tmp, glmin
410 integer :: i, ierr
411
412 tmp = huge(0.0_rp)
413 do i = 1, n
414 tmp = min(tmp,a(i))
415 end do
416 call mpi_allreduce(tmp, glmin, 1, &
417 mpi_real_precision, mpi_min, neko_comm, ierr)
418 end function glmin
419
421 function glimin(a, n)
422 integer, intent(in) :: n
423 integer, dimension(n) :: a
424 integer :: tmp, glimin
425 integer :: i, ierr
426
427 tmp = huge(0)
428 do i = 1, n
429 tmp = min(tmp,a(i))
430 end do
431 call mpi_allreduce(tmp, glimin, 1, &
432 mpi_integer, mpi_min, neko_comm, ierr)
433 end function glimin
434
435
436
437
439 subroutine chsign(a, n)
440 integer, intent(in) :: n
441 real(kind=rp), dimension(n), intent(inout) :: a
442 integer :: i
443
444 do i = 1, n
445 a(i) = -a(i)
446 end do
447
448 end subroutine chsign
449
451 function vlmax(vec,n) result(tmax)
452 integer :: n, i
453 real(kind=rp), intent(in) :: vec(n)
454 real(kind=rp) :: tmax
455 tmax = real(-99d20, rp)
456 do i = 1, n
457 tmax = max(tmax, vec(i))
458 end do
459 end function vlmax
460
462 function vlmin(vec,n) result(tmin)
463 integer, intent(in) :: n
464 real(kind=rp), intent(in) :: vec(n)
465 real(kind=rp) :: tmin
466 integer :: i
467 tmin = real(99.0e20, rp)
468 do i = 1, n
469 tmin = min(tmin, vec(i))
470 end do
471 end function vlmin
472
474 subroutine invcol1(a, n)
475 integer, intent(in) :: n
476 real(kind=rp), dimension(n), intent(inout) :: a
477 integer :: i
478
479 do i = 1, n
480 a(i) = 1.0_xp / real(a(i),xp)
481 end do
482
483 end subroutine invcol1
484
486 subroutine invcol3(a, b, c, n)
487 integer, intent(in) :: n
488 real(kind=rp), dimension(n), intent(inout) :: a
489 real(kind=rp), dimension(n), intent(in) :: b,c
490 integer :: i
491
492 do i = 1, n
493 a(i) = real(b(i),xp) / c(i)
494 end do
495
496 end subroutine invcol3
497
499 subroutine invers2(a, b, n)
500 integer, intent(in) :: n
501 real(kind=rp), dimension(n), intent(inout) :: a
502 real(kind=rp), dimension(n), intent(in) :: b
503 integer :: i
504
505 do i = 1, n
506 a(i) = 1.0_xp / real(b(i),xp)
507 end do
508
509 end subroutine invers2
510
513 subroutine vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
514 integer, intent(in) :: n
515 real(kind=rp), dimension(n), intent(in) :: v1, v2, v3
516 real(kind=rp), dimension(n), intent(in) :: w1, w2, w3
517 real(kind=rp), dimension(n), intent(out) :: u1, u2, u3
518 integer :: i
519
520 do i = 1, n
521 u1(i) = v2(i)*w3(i) - v3(i)*w2(i)
522 u2(i) = v3(i)*w1(i) - v1(i)*w3(i)
523 u3(i) = v1(i)*w2(i) - v2(i)*w1(i)
524 end do
525
526 end subroutine vcross
527
530 subroutine vdot2(dot, u1, u2, v1, v2, n)
531 integer, intent(in) :: n
532 real(kind=rp), dimension(n), intent(in) :: u1, u2
533 real(kind=rp), dimension(n), intent(in) :: v1, v2
534 real(kind=rp), dimension(n), intent(out) :: dot
535 integer :: i
536 do i = 1, n
537 dot(i) = u1(i)*v1(i) + u2(i)*v2(i)
538 end do
539
540 end subroutine vdot2
541
544 subroutine vdot3(dot, u1, u2, u3, v1, v2, v3, n)
545 integer, intent(in) :: n
546 real(kind=rp), dimension(n), intent(in) :: u1, u2, u3
547 real(kind=rp), dimension(n), intent(in) :: v1, v2, v3
548 real(kind=rp), dimension(n), intent(out) :: dot
549 integer :: i
550
551 do i = 1, n
552 dot(i) = u1(i)*v1(i) + u2(i)*v2(i) + u3(i)*v3(i)
553 end do
554
555 end subroutine vdot3
556
558 function vlsc3(u, v, w, n) result(s)
559 integer, intent(in) :: n
560 real(kind=rp), dimension(n), intent(in) :: u, v, w
561 real(kind=rp) :: s
562 integer :: i
563
564 s = 0.0_rp
565 do i = 1, n
566 s = s + u(i)*v(i)*w(i)
567 end do
568
569 end function vlsc3
570
572 function vlsc2(u, v, n) result(s)
573 integer, intent(in) :: n
574 real(kind=rp), dimension(n), intent(in) :: u, v
575 real(kind=rp) :: s
576 integer :: i
577
578 s = 0.0_rp
579 do i = 1, n
580 s = s + u(i)*v(i)
581 end do
582
583 end function vlsc2
584
586 subroutine add2(a, b, n)
587 integer, intent(in) :: n
588 real(kind=rp), dimension(n), intent(inout) :: a
589 real(kind=rp), dimension(n), intent(in) :: b
590 integer :: i
591
592 do i = 1, n
593 a(i) = a(i) + b(i)
594 end do
595
596 end subroutine add2
597
599 subroutine add3(a, b, c, n)
600 integer, intent(in) :: n
601 real(kind=rp), dimension(n), intent(inout) :: a
602 real(kind=rp), dimension(n), intent(in) :: b
603 real(kind=rp), dimension(n), intent(in) :: c
604 integer :: i
605
606 do i = 1, n
607 a(i) = b(i) + c(i)
608 end do
609
610 end subroutine add3
611
613 subroutine add4(a, b, c, d, n)
614 integer, intent(in) :: n
615 real(kind=rp), dimension(n), intent(out) :: a
616 real(kind=rp), dimension(n), intent(in) :: d
617 real(kind=rp), dimension(n), intent(in) :: c
618 real(kind=rp), dimension(n), intent(in) :: b
619 integer :: i
620
621 do i = 1, n
622 a(i) = b(i) + c(i) + d(i)
623 end do
624
625 end subroutine add4
626
628 subroutine sub2(a, b, n)
629 integer, intent(in) :: n
630 real(kind=rp), dimension(n), intent(inout) :: a
631 real(kind=rp), dimension(n), intent(inout) :: b
632 integer :: i
633
634 do i = 1, n
635 a(i) = a(i) - b(i)
636 end do
637
638 end subroutine sub2
639
641 subroutine sub3(a, b, c, n)
642 integer, intent(in) :: n
643 real(kind=rp), dimension(n), intent(inout) :: a
644 real(kind=rp), dimension(n), intent(in) :: b
645 real(kind=rp), dimension(n), intent(in) :: c
646 integer :: i
647
648 do i = 1, n
649 a(i) = b(i) - c(i)
650 end do
651
652 end subroutine sub3
653
654
657 subroutine add2s1(a, b, c1, n)
658 integer, intent(in) :: n
659 real(kind=rp), dimension(n), intent(inout) :: a
660 real(kind=rp), dimension(n), intent(inout) :: b
661 real(kind=rp), intent(in) :: c1
662 integer :: i
663
664 do i = 1, n
665 a(i) = c1 * a(i) + b(i)
666 end do
667
668 end subroutine add2s1
669
672 subroutine add2s2(a, b, c1, n)
673 integer, intent(in) :: n
674 real(kind=rp), dimension(n), intent(inout) :: a
675 real(kind=rp), dimension(n), intent(inout) :: b
676 real(kind=rp), intent(in) :: c1
677 integer :: i
678
679 do i = 1, n
680 a(i) = a(i) + c1 * b(i)
681 end do
682
683 end subroutine add2s2
684
686 subroutine addsqr2s2(a, b, c1, n)
687 integer, intent(in) :: n
688 real(kind=rp), dimension(n), intent(inout) :: a
689 real(kind=rp), dimension(n), intent(in) :: b
690 real(kind=rp), intent(in) :: c1
691 integer :: i
692
693 do i = 1,n
694 a(i) = a(i) + c1 * ( b(i) * b(i) )
695 end do
696
697 end subroutine addsqr2s2
698
700 subroutine cmult2(a, b, c, n)
701 integer, intent(in) :: n
702 real(kind=rp), dimension(n), intent(inout) :: a
703 real(kind=rp), dimension(n), intent(in) :: b
704 real(kind=rp), intent(in) :: c
705 integer :: i
706
707 do i = 1, n
708 a(i) = c * b(i)
709 end do
710
711 end subroutine cmult2
712
714 subroutine invcol2(a, b, n)
715 integer, intent(in) :: n
716 real(kind=rp), dimension(n), intent(inout) :: a
717 real(kind=rp), dimension(n), intent(in) :: b
718 integer :: i
719
720 do i = 1, n
721 a(i) = real(a(i),xp) /b(i)
722 end do
723
724 end subroutine invcol2
725
726
728 subroutine col2(a, b, n)
729 integer, intent(in) :: n
730 real(kind=rp), dimension(n), intent(inout) :: a
731 real(kind=rp), dimension(n), intent(in) :: b
732 integer :: i
733
734 do i = 1, n
735 a(i) = a(i) * b(i)
736 end do
737
738 end subroutine col2
739
741 subroutine col3(a, b, c, n)
742 integer, intent(in) :: n
743 real(kind=rp), dimension(n), intent(inout) :: a
744 real(kind=rp), dimension(n), intent(in) :: b
745 real(kind=rp), dimension(n), intent(in) :: c
746 integer :: i
747
748 do i = 1, n
749 a(i) = b(i) * c(i)
750 end do
751
752 end subroutine col3
753
755 subroutine subcol3(a, b, c, n)
756 integer, intent(in) :: n
757 real(kind=rp), dimension(n), intent(inout) :: a
758 real(kind=rp), dimension(n), intent(in) :: b
759 real(kind=rp), dimension(n), intent(in) :: c
760 integer :: i
761
762 do i = 1,n
763 a(i) = a(i) - b(i) * c(i)
764 end do
765
766 end subroutine subcol3
767
769 subroutine add3s2(a, b, c, c1, c2 ,n)
770 integer, intent(in) :: n
771 real(kind=rp), dimension(n), intent(inout) :: a
772 real(kind=rp), dimension(n), intent(in) :: b
773 real(kind=rp), dimension(n), intent(in) :: c
774 real(kind=rp), intent(in) :: c1, c2
775 integer :: i
776
777 do i = 1,n
778 a(i) = c1 * b(i) + c2 * c(i)
779 end do
780
781 end subroutine add3s2
782
783
785 subroutine subcol4(a, b, c, d, n)
786 integer, intent(in) :: n
787 real(kind=rp), dimension(n), intent(inout) :: a
788 real(kind=rp), dimension(n), intent(in) :: b
789 real(kind=rp), dimension(n), intent(in) :: c
790 real(kind=rp), dimension(n), intent(in) :: d
791 integer :: i
792
793 do i = 1,n
794 a(i) = a(i) - b(i) * c(i) * d(i)
795 end do
796
797 end subroutine subcol4
798
800 subroutine addcol3(a, b, c, n)
801 integer, intent(in) :: n
802 real(kind=rp), dimension(n), intent(inout) :: a
803 real(kind=rp), dimension(n), intent(in) :: b
804 real(kind=rp), dimension(n), intent(in) :: c
805 integer :: i
806
807 do i = 1,n
808 a(i) = a(i) + b(i) * c(i)
809 end do
810
811 end subroutine addcol3
812
814 subroutine addcol4(a, b, c, d, n)
815 integer, intent(in) :: n
816 real(kind=rp), dimension(n), intent(inout) :: a
817 real(kind=rp), dimension(n), intent(in) :: b
818 real(kind=rp), dimension(n), intent(in) :: c
819 real(kind=rp), dimension(n), intent(in) :: d
820 integer :: i
821
822 do i = 1,n
823 a(i) = a(i) + b(i) * c(i) * d(i)
824 end do
825
826 end subroutine addcol4
827
829 subroutine ascol5(a, b, c, d, e, n)
830 integer, intent(in) :: n
831 real(kind=rp), dimension(n), intent(inout) :: a
832 real(kind=rp), dimension(n), intent(in) :: b
833 real(kind=rp), dimension(n), intent(in) :: c
834 real(kind=rp), dimension(n), intent(in) :: d
835 real(kind=rp), dimension(n), intent(in) :: e
836 integer :: i
837
838 do i = 1,n
839 a(i) = b(i)*c(i)-d(i)*e(i)
840 end do
841
842 end subroutine ascol5
843
845 subroutine p_update(a, b, c, c1, c2, n)
846 integer, intent(in) :: n
847 real(kind=rp), dimension(n), intent(inout) :: a
848 real(kind=rp), dimension(n), intent(in) :: b
849 real(kind=rp), dimension(n), intent(in) :: c
850 real(kind=rp), intent(in) :: c1, c2
851 integer :: i
852
853 do i = 1,n
854 a(i) = b(i) + c1*(a(i)-c2*c(i))
855 end do
856
857 end subroutine p_update
858
860 subroutine x_update(a, b, c, c1, c2, n)
861 integer, intent(in) :: n
862 real(kind=rp), dimension(n), intent(inout) :: a
863 real(kind=rp), dimension(n), intent(in) :: b
864 real(kind=rp), dimension(n), intent(in) :: c
865 real(kind=rp), intent(in) :: c1, c2
866 integer :: i
867
868 do i = 1,n
869 a(i) = a(i) + c1*b(i)+c2*c(i)
870 end do
871
872 end subroutine x_update
873
875 function glsc2(a, b, n)
876 integer, intent(in) :: n
877 real(kind=rp), dimension(n), intent(in) :: a
878 real(kind=rp), dimension(n), intent(in) :: b
879 real(kind=rp) :: glsc2
880 real(kind=xp) :: tmp
881 integer :: i, ierr
882
883 tmp = 0.0_xp
884 do i = 1, n
885 tmp = tmp + a(i) * b(i)
886 end do
887
888 call mpi_allreduce(mpi_in_place, tmp, 1, &
889 mpi_extra_precision, mpi_sum, neko_comm, ierr)
890 glsc2 = tmp
891 end function glsc2
892
894 function glsc3(a, b, c, n)
895 integer, intent(in) :: n
896 real(kind=rp), dimension(n), intent(in) :: a
897 real(kind=rp), dimension(n), intent(in) :: b
898 real(kind=rp), dimension(n), intent(in) :: c
899 real(kind=rp) :: glsc3
900 real(kind=xp) :: tmp
901 integer :: i, ierr
902
903 tmp = 0.0_xp
904 do i = 1, n
905 tmp = tmp + a(i) * b(i) * c(i)
906 end do
907
908 call mpi_allreduce(mpi_in_place, tmp, 1, &
909 mpi_extra_precision, mpi_sum, neko_comm, ierr)
910 glsc3 = tmp
911
912 end function glsc3
913 function glsc4(a, b, c, d, n)
914 integer, intent(in) :: n
915 real(kind=rp), dimension(n), intent(in) :: a
916 real(kind=rp), dimension(n), intent(in) :: b
917 real(kind=rp), dimension(n), intent(in) :: c
918 real(kind=rp), dimension(n), intent(in) :: d
919 real(kind=rp) :: glsc4
920 real(kind=xp) :: tmp
921 integer :: i, ierr
922
923 tmp = 0.0_xp
924 do i = 1, n
925 tmp = tmp + a(i) * b(i) * c(i) * d(i)
926 end do
927
928 call mpi_allreduce(mpi_in_place, tmp, 1, &
929 mpi_extra_precision, mpi_sum, neko_comm, ierr)
930 glsc4 = tmp
931
932 end function glsc4
933
939 subroutine sortrp(a, ind, n)
940 integer, intent(in) :: n
941 real(kind=rp), intent(inout) :: a(n)
942 integer, intent(out) :: ind(n)
943 real(kind=rp) :: aa
944 integer :: j, ir, i, ii, l
945
946 do j = 1, n
947 ind(j) = j
948 end do
949
950 if (n .le. 1) return
951
952
953 l = n/2+1
954 ir = n
955 do while (.true.)
956 if (l .gt. 1) then
957 l = l-1
958 aa = a(l)
959 ii = ind(l)
960 else
961 aa = a(ir)
962 ii = ind(ir)
963 a(ir) = a(1)
964 ind(ir) = ind(1)
965 ir = ir - 1
966 if (ir .eq. 1) then
967 a(1) = aa
968 ind(1) = ii
969 return
970 end if
971 end if
972 i = l
973 j = l+l
974 do while (j .le. ir)
975 if (j .lt. ir) then
976 if ( a(j) .lt. a(j+1) ) j = j + 1
977 end if
978 if (aa .lt. a(j)) then
979 a(i) = a(j)
980 ind(i) = ind(j)
981 i = j
982 j = j+j
983 else
984 j = ir+1
985 end if
986 end do
987 a(i) = aa
988 ind(i) = ii
989 end do
990 end subroutine sortrp
991
997 subroutine sorti4(a, ind, n)
998 integer, intent(in) :: n
999 integer(i4), intent(inout) :: a(n)
1000 integer, intent(out) :: ind(n)
1001 integer(i4) :: aa
1002 integer :: j, ir, i, ii, l
1003
1004 do j = 1, n
1005 ind(j) = j
1006 end do
1007
1008 if (n .le. 1) return
1009
1010 l = n/2+1
1011 ir = n
1012 do while (.true.)
1013 if (l .gt. 1) then
1014 l = l - 1
1015 aa = a(l)
1016 ii = ind(l)
1017 else
1018 aa = a(ir)
1019 ii = ind(ir)
1020 a(ir) = a( 1)
1021 ind(ir) = ind( 1)
1022 ir = ir - 1
1023 if (ir .eq. 1) then
1024 a(1) = aa
1025 ind(1) = ii
1026 return
1027 end if
1028 end if
1029 i = l
1030 j = l + l
1031 do while (j .le. ir)
1032 if (j .lt. ir) then
1033 if ( a(j) .lt. a(j + 1) ) j = j + 1
1034 end if
1035 if (aa .lt. a(j)) then
1036 a(i) = a(j)
1037 ind(i) = ind(j)
1038 i = j
1039 j = j + j
1040 else
1041 j = ir + 1
1042 end if
1043 end do
1044 a(i) = aa
1045 ind(i) = ii
1046 end do
1047 end subroutine sorti4
1048
1053 subroutine swapdp(b, ind, n)
1054 integer, intent(in) :: n
1055 real(kind=rp), intent(inout) :: b(n)
1056 integer, intent(in) :: ind(n)
1057 real(kind=rp) :: temp(n)
1058 integer :: i, jj
1059
1060 do i = 1, n
1061 temp(i) = b(i)
1062 end do
1063 do i = 1, n
1064 jj = ind(i)
1065 b(i) = temp(jj)
1066 end do
1067 end subroutine swapdp
1068
1073 subroutine swapi4(b, ind, n)
1074 integer, intent(in) :: n
1075 integer(i4), intent(inout) :: b(n)
1076 integer, intent(in) :: ind(n)
1077 integer(i4) :: temp(n)
1078 integer :: i, jj
1079
1080 do i = 1, n
1081 temp(i) = b(i)
1082 end do
1083 do i = 1, n
1084 jj = ind(i)
1085 b(i) = temp(jj)
1086 end do
1087 end subroutine swapi4
1088
1093 subroutine reorddp(b, ind, n)
1094 integer, intent(in) :: n
1095 real(kind=rp), intent(inout) :: b(n)
1096 integer, intent(in) :: ind(n)
1097 real(kind=rp) :: temp(n)
1098 integer :: i, jj
1099
1100 do i = 1, n
1101 temp(i) = b(i)
1102 end do
1103 do i = 1, n
1104 jj = ind(i)
1105 b(jj) = temp(i)
1106 end do
1107 end subroutine reorddp
1108
1113 subroutine reordi4(b, ind, n)
1114 integer, intent(in) :: n
1115 integer(i4), intent(inout) :: b(n)
1116 integer, intent(in) :: ind(n)
1117 integer(i4) :: temp(n)
1118 integer :: i, jj
1119
1120 do i = 1, n
1121 temp(i) = b(i)
1122 end do
1123 do i = 1, n
1124 jj = ind(i)
1125 b(jj) = temp(i)
1126 end do
1127 end subroutine reordi4
1128
1133 subroutine flipvdp(b, ind, n)
1134 integer, intent(in) :: n
1135 real(kind=rp), intent(inout) :: b(n)
1136 integer, intent(inout) :: ind(n)
1137 real(kind=rp) :: temp(n)
1138 integer :: tempind(n)
1139 integer :: i, jj
1140
1141 do i = 1, n
1142 jj = n+1-i
1143 temp(jj) = b(i)
1144 tempind(jj) = ind(i)
1145 end do
1146 do i = 1,n
1147 b(i) = temp(i)
1148 ind(i) = tempind(i)
1149 end do
1150 end subroutine flipvdp
1151
1156 subroutine flipvi4(b, ind, n)
1157 integer, intent(in) :: n
1158 integer(i4), intent(inout) :: b(n)
1159 integer, intent(inout) :: ind(n)
1160 integer(i4) :: temp(n)
1161 integer :: tempind(n)
1162 integer :: i, jj
1163
1164 do i = 1, n
1165 jj = n+1-i
1166 temp(jj) = b(i)
1167 tempind(jj) = ind(i)
1168 end do
1169 do i = 1,n
1170 b(i) = temp(i)
1171 ind(i) = tempind(i)
1172 end do
1173 end subroutine flipvi4
1174
1178 subroutine absval(a, n)
1179 integer, intent(in) :: n
1180 real(kind=rp), dimension(n), intent(inout) :: a
1181 integer :: i
1182 do i = 1, n
1183 a(i) = abs(a(i))
1184 end do
1185 end subroutine absval
1186
1187 ! ========================================================================== !
1188 ! Point-wise operations
1189
1191 subroutine pwmax_vec2(a, b, n)
1192 integer, intent(in) :: n
1193 real(kind=rp), dimension(n), intent(inout) :: a
1194 real(kind=rp), dimension(n), intent(in) :: b
1195 integer :: i
1196
1197 do i = 1, n
1198 a(i) = max(a(i), b(i))
1199 end do
1200 end subroutine pwmax_vec2
1201
1203 subroutine pwmax_vec3(a, b, c, n)
1204 integer, intent(in) :: n
1205 real(kind=rp), dimension(n), intent(inout) :: a
1206 real(kind=rp), dimension(n), intent(in) :: b, c
1207 integer :: i
1208
1209 do i = 1, n
1210 a(i) = max(b(i), c(i))
1211 end do
1212 end subroutine pwmax_vec3
1213
1215 subroutine pwmax_scal2(a, b, n)
1216 integer, intent(in) :: n
1217 real(kind=rp), dimension(n), intent(inout) :: a
1218 real(kind=rp), intent(in) :: b
1219 integer :: i
1220
1221 do i = 1, n
1222 a(i) = max(a(i), b)
1223 end do
1224 end subroutine pwmax_scal2
1225
1227 subroutine pwmax_scal3(a, b, c, n)
1228 integer, intent(in) :: n
1229 real(kind=rp), dimension(n), intent(inout) :: a
1230 real(kind=rp), dimension(n), intent(in) :: b
1231 real(kind=rp), intent(in) :: c
1232 integer :: i
1233
1234 do i = 1, n
1235 a(i) = max(b(i), c)
1236 end do
1237 end subroutine pwmax_scal3
1238
1240 subroutine pwmin_vec2(a, b, n)
1241 integer, intent(in) :: n
1242 real(kind=rp), dimension(n), intent(inout) :: a
1243 real(kind=rp), dimension(n), intent(in) :: b
1244 integer :: i
1245
1246 do i = 1, n
1247 a(i) = min(a(i), b(i))
1248 end do
1249 end subroutine pwmin_vec2
1250
1252 subroutine pwmin_vec3(a, b, c, n)
1253 integer, intent(in) :: n
1254 real(kind=rp), dimension(n), intent(inout) :: a
1255 real(kind=rp), dimension(n), intent(in) :: b, c
1256 integer :: i
1257
1258 do i = 1, n
1259 a(i) = min(b(i), c(i))
1260 end do
1261 end subroutine pwmin_vec3
1262
1264 subroutine pwmin_sca2(a, b, n)
1265 integer, intent(in) :: n
1266 real(kind=rp), dimension(n), intent(inout) :: a
1267 real(kind=rp), intent(in) :: b
1268 integer :: i
1269
1270 do i = 1, n
1271 a(i) = min(a(i), b)
1272 end do
1273 end subroutine pwmin_sca2
1274
1276 subroutine pwmin_sca3(a, b, c, n)
1277 integer, intent(in) :: n
1278 real(kind=rp), dimension(n), intent(inout) :: a
1279 real(kind=rp), dimension(n), intent(in) :: b
1280 real(kind=rp), intent(in) :: c
1281 integer :: i
1282
1283 do i = 1, n
1284 a(i) = min(b(i), c)
1285 end do
1286 end subroutine pwmin_sca3
1287
1288end 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:311
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
Definition math.f90:701
subroutine, public row_zero(a, m, n, e)
Sets row e to 0 in matrix a.
Definition math.f90:217
subroutine, public invcol2(a, b, n)
Vector division .
Definition math.f90:715
real(kind=rp) function, public vlsc2(u, v, n)
Compute multiplication sum .
Definition math.f90:573
subroutine pwmin_sca3(a, b, c, n)
Point-wise minimum of scalar and vector .
Definition math.f90:1277
real(kind=rp), parameter, public pi
Definition math.f90:76
subroutine pwmax_scal2(a, b, n)
Point-wise maximum of scalar and vector .
Definition math.f90:1216
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition math.f90:895
pure logical function qabscmp(x, y)
Return double precision absolute comparison .
Definition math.f90:142
subroutine, public ascol5(a, b, c, d, e, n)
Returns .
Definition math.f90:830
subroutine, public invers2(a, b, n)
Compute inverted vector .
Definition math.f90:500
subroutine, public cadd2(a, b, s, n)
Add a scalar to vector .
Definition math.f90:335
subroutine pwmin_sca2(a, b, n)
Point-wise minimum of scalar and vector .
Definition math.f90:1265
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:323
subroutine reorddp(b, ind, n)
reorder double precision array - inverse of swap
Definition math.f90:1094
subroutine, public addsqr2s2(a, b, c1, n)
Returns .
Definition math.f90:687
real(kind=rp) function, public glsc4(a, b, c, d, n)
Definition math.f90:914
subroutine, public masked_red_copy(a, b, mask, n, m)
Copy a masked vector to reduced contigous vector .
Definition math.f90:280
subroutine pwmax_vec3(a, b, c, n)
Point-wise maximum of two vectors .
Definition math.f90:1204
subroutine swapdp(b, ind, n)
sort double precision array acording to ind vector
Definition math.f90:1054
subroutine flipvi4(b, ind, n)
Flip single integer vector b and ind.
Definition math.f90:1157
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
Definition math.f90:658
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition math.f90:876
subroutine, public subcol3(a, b, c, n)
Returns .
Definition math.f90:756
subroutine, public rone(a, n)
Set all elements to one.
Definition math.f90:228
subroutine flipvdp(b, ind, n)
Flip double precision vector b and ind.
Definition math.f90:1134
subroutine, public x_update(a, b, c, c1, c2, n)
Returns .
Definition math.f90:861
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:600
subroutine swapi4(b, ind, n)
sort single integer array acording to ind vector
Definition math.f90:1074
integer function, public glimin(a, n)
Min of an integer vector of length n.
Definition math.f90:422
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
Definition math.f90:360
subroutine, public sub3(a, b, c, n)
Vector subtraction .
Definition math.f90:642
subroutine, public addcol4(a, b, c, d, n)
Returns .
Definition math.f90:815
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:587
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition math.f90:348
subroutine, public absval(a, n)
Take the absolute value of an array.
Definition math.f90:1179
subroutine, public invcol3(a, b, c, n)
Invert a vector .
Definition math.f90:487
pure logical function sabscmp(x, y)
Return single precision absolute comparison .
Definition math.f90:122
subroutine, public add3s2(a, b, c, c1, c2, n)
Returns .
Definition math.f90:770
subroutine, public subcol4(a, b, c, d, n)
Returns .
Definition math.f90:786
subroutine sorti4(a, ind, n)
Heap Sort for single integer arrays.
Definition math.f90:998
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:801
pure logical function dabscmp(x, y)
Return double precision absolute comparison .
Definition math.f90:132
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:475
subroutine, public masked_copy(a, b, mask, n, m)
Copy a masked vector .
Definition math.f90:258
real(kind=rp), parameter, public neko_m_ln2
Definition math.f90:73
subroutine, public chsign(a, n)
Change sign of vector .
Definition math.f90:440
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:729
subroutine, public izero(a, n)
Zero an integer vector.
Definition math.f90:206
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition math.f90:377
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:239
subroutine pwmax_scal3(a, b, c, n)
Point-wise maximum of scalar and vector .
Definition math.f90:1228
subroutine, public add4(a, b, c, d, n)
Vector addition .
Definition math.f90:614
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:742
subroutine pwmin_vec2(a, b, n)
Point-wise minimum of two vectors .
Definition math.f90:1241
subroutine pwmin_vec3(a, b, c, n)
Point-wise minimum of two vectors .
Definition math.f90:1253
subroutine pwmax_vec2(a, b, n)
Point-wise maximum of two vectors .
Definition math.f90:1192
pure logical function drelcmp(x, y, eps)
Return double precision relative comparison .
Definition math.f90:166
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:70
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:545
pure logical function qrelcmp(x, y, eps)
Return quad precision relative comparison .
Definition math.f90:181
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:195
subroutine, public vdot2(dot, u1, u2, v1, v2, n)
Compute a dot product (2-d version) assuming vector components etc.
Definition math.f90:531
real(kind=rp) function, public vlmin(vec, n)
minimun value of a vector of length n
Definition math.f90:463
real(kind=rp) function, public vlmax(vec, n)
maximum value of a vector of length n
Definition math.f90:452
integer function, public glimax(a, n)
Max of an integer vector of length n.
Definition math.f90:392
subroutine sortrp(a, ind, n)
Heap Sort for double precision arrays.
Definition math.f90:940
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:629
subroutine, public cfill_mask(a, c, size, mask, mask_size)
Fill a constant to a masked vector. .
Definition math.f90:297
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition math.f90:673
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition math.f90:407
subroutine, public vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Compute a cross product assuming vector components etc.
Definition math.f90:514
pure logical function srelcmp(x, y, eps)
Return single precision relative comparison .
Definition math.f90:152
real(kind=rp) function, public vlsc3(u, v, w, n)
Compute multiplication sum .
Definition math.f90:559
subroutine reordi4(b, ind, n)
reorder single integer array - inverse of swap
Definition math.f90:1114
subroutine, public p_update(a, b, c, c1, c2, n)
Returns .
Definition math.f90:846
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
Utilities.
Definition utils.f90:35
#define max(a, b)
Definition tensor.cu:40