Choral
R1d_mod.F90
Go to the documentation of this file.
1 
9 
10 module r1d_mod
11 
12  use real_type, only: rp, equal
13  use basic_tools
14  !$ use OMP_LIB
15 
16  implicit none
17  private
18 
19  ! %----------------------------------------%
20  ! | |
21  ! | PUBLIC DATA |
22  ! | |
23  ! %----------------------------------------%
24  ! GENERIC
25  public :: copy !! TESTED
26  public :: scale !! TESTED
27  public :: mult !! TESTED
28  public :: axpy !! TESTED
29  public :: xpay !! TESTED
30  public :: axpby !! TESTED
31  public :: scalprod !! TESTED
32  public :: norm_2 !! TESTED
33 
34  ! NON GENERIC
35  public :: reorder !! TESTED
36 
37  ! %----------------------------------------%
38  ! | |
39  ! | GENERIc SUBROUTINES |
40  ! | |
41  ! %----------------------------------------%
42 
44  interface copy
45  module procedure r1d_copy !! TESTED
46  end interface copy
47 
49  interface scale
50  module procedure ax_1 !! TESTED
51  module procedure ax_2 !! TESTED
52  end interface scale
53 
55  interface axpy
56  module procedure r1d_axpy !! TESTED
57  end interface axpy
58 
60  interface xpay
61  module procedure xpay_1 !! TESTED
62  module procedure xpay_2 !! TESTED
63  end interface xpay
64 
66  interface axpby
67  module procedure axpby_1 !! TESTED
68  module procedure axpby_2 !! TESTED
69  end interface axpby
70 
72  interface reorder
73  module procedure r1d_reorder
74  end interface reorder
75 
77  interface scalprod
78  module procedure r1d_scalprod !! TESTED
79  end interface scalprod
80 
82  interface norm_2
83  module procedure r1d_norm_2 !! TESTED
84  end interface norm_2
85 
87  interface mult
88  module procedure xy_2 !! TESTED
89  end interface mult
90 
91 contains
92 
95  subroutine r1d_axpy(a, x, y)
96  real(RP), dimension(:), intent(inout) :: x
97  real(RP) , intent(in) :: a
98  real(RP), dimension(:), intent(in) :: y
99 
100  integer :: ii, n2
101 
102 #if DBG>0
103  if ( size(y,1) < size(x,1)) call quit(&
104  &"R1d_mod: R1d_axpy: 1" )
105 #endif
106 
107  n2 = size(x,1)
108  !$OMP PARALLEL
109  !$OMP DO
110  do ii=1, n2
111  x(ii) = a*x(ii) + y(ii)
112  end do
113  !$OMP END DO
114  !$OMP END PARALLEL
115 
116  end subroutine r1d_axpy
117 
120  subroutine xpay_1(x, a, y)
121  real(RP), dimension(:), intent(inout) :: x
122  real(RP) , intent(in) :: a
123  real(RP), dimension(:), intent(in) :: y
124 
125  integer :: ii, n2
126 
127 #if DBG>0
128  if ( size(y,1) < size(x,1)) call quit(&
129  &"R1d_mod: xpay_1: 1" )
130 #endif
131 
132  n2 = size(x,1)
133  if (equal(a,1._rp)) then
134  !$OMP PARALLEL
135  !$OMP DO
136  do ii=1, n2
137  x(ii) = x(ii) + y(ii)
138  end do
139  !$OMP END DO
140  !$OMP END PARALLEL
141  else
142  !$OMP PARALLEL
143  !$OMP DO
144  do ii=1, n2
145  x(ii) = x(ii) + a*y(ii)
146  end do
147  !$OMP END DO
148  !$OMP END PARALLEL
149  end if
150 
151  end subroutine xpay_1
152 
153 
156  subroutine xpay_2(z, x, a, y)
157  real(RP), dimension(:), intent(out) :: z
158  real(RP) , intent(in) :: a
159  real(RP), dimension(:), intent(in) :: x, y
160 
161  integer :: ii, n2
162 
163 #if DBG>0
164  if ( size(x,1) < size(z,1)) call quit(&
165  &"R1d_mod: xpay_2: 1" )
166  if ( size(y,1) < size(z,1)) call quit(&
167  &"R1d_mod: xpay_2: 2" )
168 #endif
169 
170  n2 = size(z,1)
171  if (equal(a,1._rp)) then
172  !$OMP PARALLEL
173  !$OMP DO
174  do ii=1, n2
175  z(ii) = x(ii) + y(ii)
176  end do
177  !$OMP END DO
178  !$OMP END PARALLEL
179 
180  else
181  !$OMP PARALLEL
182  !$OMP DO
183  do ii=1, n2
184  z(ii) = x(ii) + a*y(ii)
185  end do
186  !$OMP END DO
187  !$OMP END PARALLEL
188 
189  end if
190 
191  end subroutine xpay_2
192 
193 
196  subroutine axpby_1(a, x, b, y)
197  real(RP), dimension(:), intent(inout) :: x
198  real(RP) , intent(in) :: a, b
199  real(RP), dimension(:), intent(in) :: y
200 
201  integer :: ii, n2
202 
203 #if DBG>0
204  if ( size(y,1) < size(x,1)) call quit(&
205  &"R1d_mod: axpby_1: 1" )
206 #endif
207 
208  n2 = size(x,1)
209  !$OMP PARALLEL
210  !$OMP DO
211  do ii=1, n2
212  x(ii) = a*x(ii) + b*y(ii)
213  end do
214  !$OMP END DO
215  !$OMP END PARALLEL
216 
217  end subroutine axpby_1
218 
221  subroutine axpby_2(z, a, x, b, y)
222  real(RP), dimension(:), intent(out) :: z
223  real(RP) , intent(in) :: a, b
224  real(RP), dimension(:), intent(in) :: x, y
225 
226  integer :: ii, n2
227 
228 #if DBG>0
229  if ( size(x,1) < size(z,1)) call quit(&
230  &"R1d_mod: axpby_2: 1" )
231  if ( size(y,1) < size(z,1)) call quit(&
232  &"R1d_mod: axpby_2: 2" )
233 #endif
234 
235  n2 = size(z,1)
236  !$OMP PARALLEL
237  !$OMP DO
238  do ii=1, n2
239  z(ii) = a*x(ii) + b*y(ii)
240  end do
241  !$OMP END DO
242  !$OMP END PARALLEL
243 
244  end subroutine axpby_2
245 
248  subroutine ax_1(x, a)
249  real(RP), dimension(:), intent(inout) :: x
250  real(RP) , intent(in) :: a
251 
252  integer :: ii, n2
253 
254  n2 = size(x,1)
255  !$OMP PARALLEL
256  !$OMP DO
257  do ii=1, n2
258  x(ii) = a*x(ii)
259  end do
260  !$OMP END DO
261  !$OMP END PARALLEL
262 
263  end subroutine ax_1
264 
265 
268  subroutine ax_2(y, a, x)
269  real(RP), dimension(:), intent(out) :: y
270  real(RP) , intent(in) :: a
271  real(RP), dimension(:), intent(in) :: x
272 
273  integer :: ii, n2
274 
275 #if DBG>0
276  if ( size(x,1) < size(y,1)) call quit(&
277  &"R1d_mod: ax_2: 1" )
278 #endif
279 
280  n2 = size(y,1)
281  !$OMP PARALLEL
282  !$OMP DO
283  do ii=1, n2
284  y(ii) = a*x(ii)
285  end do
286  !$OMP END DO
287  !$OMP END PARALLEL
288 
289  end subroutine ax_2
290 
291 
294  subroutine r1d_copy(y, x)
295  real(RP), dimension(:), intent(out) :: y
296  real(RP), dimension(:), intent(in) :: x
297 
298  integer :: ii, n2
299 
300 #if DBG>0
301  if ( size(x,1) < size(y,1)) call quit(&
302  &"R1d_mod: R1d_copy: 1" )
303 #endif
304 
305  n2 = size(y,1)
306  !$OMP PARALLEL
307  !$OMP DO
308  do ii=1, n2
309  y(ii) = x(ii)
310  end do
311  !$OMP END DO
312  !$OMP END PARALLEL
313 
314  end subroutine r1d_copy
315 
318  function r1d_scalprod(x, y) result(res)
319  real(RP), dimension(:), intent(in) :: x, y
320  real(RP) :: res
321  integer :: ii, n2
322 
323 #if DBG>0
324  if ( size(x,1) /= size(y,1)) call quit(&
325  &"R1d_mod: R1d_scalProd: 1" )
326 #endif
327 
328  res = 0._rp
329  n2 = size(x,1)
330  !$OMP PARALLEL
331  !$OMP DO reduction(+:res)
332  do ii=1, n2
333  res = res + x(ii) * y(ii)
334  end do
335  !$OMP END DO
336  !$OMP END PARALLEL
337 
338  end function r1d_scalprod
339 
342  function r1d_norm_2(x) result(res)
343  real(RP), dimension(:), intent(in) :: x
344  real(RP) :: res
345 
346  res = scalprod(x,x)
347  res = sqrt(res)
348 
349  end function r1d_norm_2
350 
351 
354  subroutine r1d_reorder(y, sigma, x)
355  real(RP), dimension(:), intent(out) :: y
356  real(RP), dimension(:), intent(in) :: x
357  integer , dimension(:), intent(in) :: sigma
358 
359  integer :: n2, ii
360 
361 #if DBG>0
362  if ( size(x,1) /= size(y,1)) call quit(&
363  &"R1d_mod: : R1d_reorder: 1" )
364  if ( size(x,1) /= size(sigma,1)) call quit(&
365  &"R1d_mod: : R1d_reorder: 2" )
366 #endif
367 
368  n2 = size(x, 1)
369 
370  !$OMP PARALLEL
371  !$OMP DO
372  do ii=1, n2
373  y(ii) = x(sigma(ii))
374  end do
375  !$OMP END DO
376  !$OMP END PARALLEL
377 
378  end subroutine r1d_reorder
379 
380 
383  subroutine xy_2(z, x, y)
384  real(RP), dimension(:), intent(out) :: z
385  real(RP), dimension(:), intent(in) :: x, y
386 
387  integer :: ii, n2
388 
389 #if DBG>0
390  if ( size(x,1) < size(z,1)) call quit(&
391  &"R1d_mod: xy_2: 1" )
392  if ( size(y,1) < size(z,1)) call quit(&
393  &"R1d_mod: xy_2: 2" )
394 #endif
395 
396  n2 = size(z,1)
397  !$OMP PARALLEL
398  !$OMP DO
399  do ii=1, n2
400  z(ii) = x(ii) * y(ii)
401  end do
402  !$OMP END DO
403  !$OMP END PARALLEL
404 
405  end subroutine xy_2
406 
407 
408 end module r1d_mod
subroutine xpay_2(z, x, a, y)
z(:) := x(:) + a*y(:)
Definition: R1d_mod.F90:157
subroutine r1d_axpy(a, x, y)
x(:) := a*x(:) + y(:)
Definition: R1d_mod.F90:96
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine ax_1(x, a)
x(:) := a*x(:)
Definition: R1d_mod.F90:249
x = a*x // OR // y = a*x
Definition: R1d_mod.F90:49
subroutine ax_2(y, a, x)
y(:) := a*x(:)
Definition: R1d_mod.F90:269
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
Definition: real_type.F90:90
subroutine, public quit(message)
Stop program execution, display an error messahe.
L2 vector norm.
Definition: R1d_mod.F90:82
z = x*y
Definition: R1d_mod.F90:87
Test real equality.
Definition: real_type.F90:146
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
subroutine xpay_1(x, a, y)
x(:) := x(:) + a*y(:)
Definition: R1d_mod.F90:121
y = x
Definition: R1d_mod.F90:44
x = a*x + y
Definition: R1d_mod.F90:55
scalar product
Definition: R1d_mod.F90:77
subroutine r1d_copy(y, x)
y(:) := x(:)
Definition: R1d_mod.F90:295
real(rp) function r1d_scalprod(x, y)
scalar product
Definition: R1d_mod.F90:319
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
subroutine xy_2(z, x, y)
z(:) := x(:)*y(:)
Definition: R1d_mod.F90:384
real(rp) function r1d_norm_2(x)
Euclidian norm.
Definition: R1d_mod.F90:343
z = a*x + b*y
Definition: R1d_mod.F90:66
subroutine r1d_reorder(y, sigma, x)
y(ii) := x(sigma(ii))
Definition: R1d_mod.F90:355
subroutine axpby_1(a, x, b, y)
x(:) := a*x(:) + b*y(:)
Definition: R1d_mod.F90:197
subroutine axpby_2(z, a, x, b, y)
z(:) := a*x(:) + b*y(:)
Definition: R1d_mod.F90:222
array entry permutation
Definition: R1d_mod.F90:72