Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
tensor_cpu.f90
Go to the documentation of this file.
1module tensor_cpu
2 use num_types, only : rp
3 use mxm_wrapper, only : mxm
4 implicit none
5 private
6
7 public :: tnsr2d_el_cpu, tnsr3d_el_cpu, tnsr3d_cpu, tnsr1_3d_cpu
8
9contains
10
22 subroutine tnsr2d_el_cpu(v, nv, u, nu, A, Bt)
23 integer, intent(in) :: nv, nu
24 real(kind=rp), intent(inout) :: v(nv*nv), u(nu*nu)
25 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv)
26 real(kind=rp) :: work(0:nu**2*nv)
27
28 call mxm(a, nv, u, nu, work, nu)
29 call mxm(work, nv, bt, nu, v, nv)
30
31 end subroutine tnsr2d_el_cpu
32
49 subroutine tnsr3d_el_cpu(v, nv, u, nu, A, Bt, Ct)
50 integer, intent(in) :: nv, nu
51 real(kind=rp), intent(inout) :: v(nv*nv*nv), u(nu*nu*nu)
52 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
53
54 if (nv .eq. nu) then
55 select case (nv)
56 case (14)
57 call tnsr3d_el_n14_cpu(v, u, a, bt, ct)
58 case (13)
59 call tnsr3d_el_n13_cpu(v, u, a, bt, ct)
60 case (12)
61 call tnsr3d_el_n12_cpu(v, u, a, bt, ct)
62 case (11)
63 call tnsr3d_el_n11_cpu(v, u, a, bt, ct)
64 case (10)
65 call tnsr3d_el_n10_cpu(v, u, a, bt, ct)
66 case (9)
67 call tnsr3d_el_n9_cpu(v, u, a, bt, ct)
68 case (8)
69 call tnsr3d_el_n8_cpu(v, u, a, bt, ct)
70 case (7)
71 call tnsr3d_el_n7_cpu(v, u, a, bt, ct)
72 case (6)
73 call tnsr3d_el_n6_cpu(v, u, a, bt, ct)
74 case (5)
75 call tnsr3d_el_n5_cpu(v, u, a, bt, ct)
76 case (4)
77 call tnsr3d_el_n4_cpu(v, u, a, bt, ct)
78 case (3)
79 call tnsr3d_el_n3_cpu(v, u, a, bt, ct)
80 case (2)
81 call tnsr3d_el_n2_cpu(v, u, a, bt, ct)
82 case default
83 call tnsr3d_el_n_cpu(v, u, a, bt, ct, nv)
84 end select
85 else if (nv .eq. 1) then
86 select case (nu)
87 case (4)
88 call tnsr3d_el_1_4_cpu(v, u, a, bt, ct)
89 case (6)
90 call tnsr3d_el_1_6_cpu(v, u, a, bt, ct)
91 case (8)
92 call tnsr3d_el_1_8_cpu(v, u, a, bt, ct)
93 case (10)
94 call tnsr3d_el_1_10_cpu(v, u, a, bt, ct)
95 case (12)
96 call tnsr3d_el_1_12_cpu(v, u, a, bt, ct)
97 case default
98 call tnsr3d_el_1_nu_cpu(v, u, nu, a, bt, ct)
99 end select
100 else
101 call tnsr3d_el_nvnu_cpu(v, nv, u, nu, a, bt, ct)
102 end if
103
104 end subroutine tnsr3d_el_cpu
105
119 !OCL SERIAL
120 subroutine tnsr3d_el_nvnu_cpu(v, nv, u, nu, A, Bt, Ct)
121 integer, intent(in) :: nv, nu
122 real(kind=rp), intent(inout) :: v(nv*nv*nv), u(nu*nu*nu)
123 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
124 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2)
125 real(kind=rp) :: tmp
126 integer :: i, j, k, l, nunu, nvnu, nvnv
127 integer :: ii, jj
128 nvnu = nv * nu
129 nunu = nu * nu
130 nvnv = nv * nv
131
132 do j = 1, nunu
133 do i = 1, nv
134 ii = i + nv * (j - 1)
135 tmp = 0.0_rp
136 do k = 1, nu
137 tmp = tmp + a(i,k) * u(k + nu * (j - 1))
138 end do
139 work(ii) = tmp
140 end do
141 end do
142
143 do i = 1, nu
144 do j = 1, nv
145 do l = 1, nv
146 ii = l + nv * (j - 1) + nvnv * (i - 1)
147 tmp = 0.0_rp
148 do k = 1, nu
149 jj = l + nv * (k - 1) + nvnu * (i - 1)
150 tmp = tmp + work(jj) * bt(k,j)
151 end do
152 work2(ii) = tmp
153 end do
154 end do
155 end do
156
157 do j = 1, nv
158 do i = 1, nvnv
159 jj = i + nvnv * (j - 1)
160 tmp = 0.0_rp
161 do k = 1, nu
162 ii = i + nvnv * (k - 1)
163 tmp = tmp + work2(ii) * ct(k, j)
164 end do
165 v(jj) = tmp
166 end do
167 end do
168
169 end subroutine tnsr3d_el_nvnu_cpu
170
181 !OCL SERIAL
182 subroutine tnsr3d_el_1_nu_cpu(v, u, nu, A, Bt, Ct)
183 integer, intent(in) :: nu
184 real(kind=rp), intent(inout) :: v(1)
185 real(kind=rp), intent(in) :: u(nu*nu*nu)
186 real(kind=rp), intent(in) :: a(1, nu), bt(nu, 1), ct(nu, 1)
187 real(kind=rp) :: work(nu**2), work2(nu)
188 real(kind=rp) :: tmp
189 integer :: i, j, k, nunu
190 integer :: jj
191 nunu = nu * nu
192
193 do j = 1, nunu
194 tmp = 0.0_rp
195 do k = 1, nu
196 tmp = tmp + a(1,k) * u(k + nu * (j - 1))
197 end do
198 work(j) = tmp
199 end do
200
201 do i = 1, nu
202 tmp = 0.0_rp
203 do k = 1, nu
204 jj = k + nu * (i - 1)
205 tmp = tmp + work(jj) * bt(k,1)
206 end do
207 work2(i) = tmp
208 end do
209
210 tmp = 0.0_rp
211 do k = 1, nu
212 tmp = tmp + work2(k) * ct(k, 1)
213 end do
214 v(1) = tmp
215
216 end subroutine tnsr3d_el_1_nu_cpu
217
221 !OCL SERIAL
222 subroutine tnsr3d_el_1_4_cpu(v, u, A, Bt, Ct)
223 integer, parameter :: n = 4
224 integer, parameter :: nn = n**2
225 real(kind=rp), intent(inout) :: v(1)
226 real(kind=rp), intent(in) :: u(n*n*n)
227 real(kind=rp), intent(in) :: a(1,n), bt(n,1), ct(n,1)
228 real(kind=rp) :: work(n**2), work2(n)
229 integer :: i, j
230
231 do j = 1, nn
232 work(j) = a(1,1) * u(1 + n * (j - 1)) &
233 + a(1,2) * u(2 + n * (j - 1)) &
234 + a(1,3) * u(3 + n * (j - 1)) &
235 + a(1,4) * u(4 + n * (j - 1))
236 end do
237
238 do i = 1, n
239 work2(i) = work(1 + n * (i - 1)) * bt(1,1) &
240 + work(2 + n * (i - 1)) * bt(2,1) &
241 + work(3 + n * (i - 1)) * bt(3,1) &
242 + work(4 + n * (i - 1)) * bt(4,1)
243 end do
244
245 v(1) = work2(1) * ct(1, 1) &
246 + work2(2) * ct(2, 1) &
247 + work2(3) * ct(3, 1) &
248 + work2(4) * ct(4, 1)
249
250 end subroutine tnsr3d_el_1_4_cpu
251
255 !OCL SERIAL
256 subroutine tnsr3d_el_1_6_cpu(v, u, A, Bt, Ct)
257 integer, parameter :: n = 6
258 integer, parameter :: nn = n**2
259 real(kind=rp), intent(inout) :: v(1)
260 real(kind=rp), intent(in) :: u(n*n*n)
261 real(kind=rp), intent(in) :: a(1,n), bt(n,1), ct(n,1)
262 real(kind=rp) :: work(n**2), work2(n)
263 integer :: i, j
264
265 do j = 1, nn
266 work(j) = a(1,1) * u(1 + n * (j - 1)) &
267 + a(1,2) * u(2 + n * (j - 1)) &
268 + a(1,3) * u(3 + n * (j - 1)) &
269 + a(1,4) * u(4 + n * (j - 1)) &
270 + a(1,5) * u(5 + n * (j - 1)) &
271 + a(1,6) * u(6 + n * (j - 1))
272 end do
273
274 do i = 1, n
275 work2(i) = work(1 + n * (i - 1)) * bt(1,1) &
276 + work(2 + n * (i - 1)) * bt(2,1) &
277 + work(3 + n * (i - 1)) * bt(3,1) &
278 + work(4 + n * (i - 1)) * bt(4,1) &
279 + work(5 + n * (i - 1)) * bt(5,1) &
280 + work(6 + n * (i - 1)) * bt(6,1)
281 end do
282
283 v(1) = work2(1) * ct(1, 1) &
284 + work2(2) * ct(2, 1) &
285 + work2(3) * ct(3, 1) &
286 + work2(4) * ct(4, 1) &
287 + work2(5) * ct(5, 1) &
288 + work2(6) * ct(6, 1)
289
290 end subroutine tnsr3d_el_1_6_cpu
291
295 !OCL SERIAL
296 subroutine tnsr3d_el_1_8_cpu(v, u, A, Bt, Ct)
297 integer, parameter :: n = 8
298 integer, parameter :: nn = n**2
299 real(kind=rp), intent(inout) :: v(1)
300 real(kind=rp), intent(in) :: u(n*n*n)
301 real(kind=rp), intent(in) :: a(1,n), bt(n,1), ct(n,1)
302 real(kind=rp) :: work(n**2), work2(n)
303 integer :: i, j
304
305 do j = 1, nn
306 work(j) = a(1,1) * u(1 + n * (j - 1)) &
307 + a(1,2) * u(2 + n * (j - 1)) &
308 + a(1,3) * u(3 + n * (j - 1)) &
309 + a(1,4) * u(4 + n * (j - 1)) &
310 + a(1,5) * u(5 + n * (j - 1)) &
311 + a(1,6) * u(6 + n * (j - 1)) &
312 + a(1,7) * u(7 + n * (j - 1)) &
313 + a(1,8) * u(8 + n * (j - 1))
314 end do
315
316 do i = 1, n
317 work2(i) = work(1 + n * (i - 1)) * bt(1,1) &
318 + work(2 + n * (i - 1)) * bt(2,1) &
319 + work(3 + n * (i - 1)) * bt(3,1) &
320 + work(4 + n * (i - 1)) * bt(4,1) &
321 + work(5 + n * (i - 1)) * bt(5,1) &
322 + work(6 + n * (i - 1)) * bt(6,1) &
323 + work(7 + n * (i - 1)) * bt(7,1) &
324 + work(8 + n * (i - 1)) * bt(8,1)
325 end do
326
327 v(1) = work2(1) * ct(1, 1) &
328 + work2(2) * ct(2, 1) &
329 + work2(3) * ct(3, 1) &
330 + work2(4) * ct(4, 1) &
331 + work2(5) * ct(5, 1) &
332 + work2(6) * ct(6, 1) &
333 + work2(7) * ct(7, 1) &
334 + work2(8) * ct(8, 1)
335
336
337 end subroutine tnsr3d_el_1_8_cpu
338
342 !OCL SERIAL
343 subroutine tnsr3d_el_1_10_cpu(v, u, A, Bt, Ct)
344 integer, parameter :: n = 10
345 integer, parameter :: nn = n**2
346 real(kind=rp), intent(inout) :: v(1)
347 real(kind=rp), intent(in) :: u(n*n*n)
348 real(kind=rp), intent(in) :: a(1,n), bt(n,1), ct(n,1)
349 real(kind=rp) :: work(n**2), work2(n)
350 integer :: i, j
351
352 do j = 1, nn
353 work(j) = a(1,1) * u(1 + n * (j - 1)) &
354 + a(1,2) * u(2 + n * (j - 1)) &
355 + a(1,3) * u(3 + n * (j - 1)) &
356 + a(1,4) * u(4 + n * (j - 1)) &
357 + a(1,5) * u(5 + n * (j - 1)) &
358 + a(1,6) * u(6 + n * (j - 1)) &
359 + a(1,7) * u(7 + n * (j - 1)) &
360 + a(1,8) * u(8 + n * (j - 1)) &
361 + a(1,9) * u(9 + n * (j - 1)) &
362 + a(1,10) * u(10 + n * (j - 1))
363 end do
364
365 do i = 1, n
366 work2(i) = work(1 + n * (i - 1)) * bt(1,1) &
367 + work(2 + n * (i - 1)) * bt(2,1) &
368 + work(3 + n * (i - 1)) * bt(3,1) &
369 + work(4 + n * (i - 1)) * bt(4,1) &
370 + work(5 + n * (i - 1)) * bt(5,1) &
371 + work(6 + n * (i - 1)) * bt(6,1) &
372 + work(7 + n * (i - 1)) * bt(7,1) &
373 + work(8 + n * (i - 1)) * bt(8,1) &
374 + work(9 + n * (i - 1)) * bt(9,1) &
375 + work(10 + n * (i - 1)) * bt(10,1)
376 end do
377
378 v(1) = work2(1) * ct(1, 1) &
379 + work2(2) * ct(2, 1) &
380 + work2(3) * ct(3, 1) &
381 + work2(4) * ct(4, 1) &
382 + work2(5) * ct(5, 1) &
383 + work2(6) * ct(6, 1) &
384 + work2(7) * ct(7, 1) &
385 + work2(8) * ct(8, 1) &
386 + work2(9) * ct(9, 1) &
387 + work2(10) * ct(10, 1)
388
389 end subroutine tnsr3d_el_1_10_cpu
390
394 !OCL SERIAL
395 subroutine tnsr3d_el_1_12_cpu(v, u, A, Bt, Ct)
396 integer, parameter :: n = 12
397 integer, parameter :: nn = n**2
398 real(kind=rp), intent(inout) :: v(1)
399 real(kind=rp), intent(in) :: u(n*n*n)
400 real(kind=rp), intent(in) :: a(1,n), bt(n,1), ct(n,1)
401 real(kind=rp) :: work(n**2), work2(n)
402 integer :: i, j
403
404 do j = 1, nn
405 work(j) = a(1,1) * u(1 + n * (j - 1)) &
406 + a(1,2) * u(2 + n * (j - 1)) &
407 + a(1,3) * u(3 + n * (j - 1)) &
408 + a(1,4) * u(4 + n * (j - 1)) &
409 + a(1,5) * u(5 + n * (j - 1)) &
410 + a(1,6) * u(6 + n * (j - 1)) &
411 + a(1,7) * u(7 + n * (j - 1)) &
412 + a(1,8) * u(8 + n * (j - 1)) &
413 + a(1,9) * u(9 + n * (j - 1)) &
414 + a(1,10) * u(10 + n * (j - 1)) &
415 + a(1,11) * u(11 + n * (j - 1)) &
416 + a(1,12) * u(12 + n * (j - 1))
417 end do
418
419 do i = 1, n
420 work2(i) = work(1 + n * (i - 1)) * bt(1,1) &
421 + work(2 + n * (i - 1)) * bt(2,1) &
422 + work(3 + n * (i - 1)) * bt(3,1) &
423 + work(4 + n * (i - 1)) * bt(4,1) &
424 + work(5 + n * (i - 1)) * bt(5,1) &
425 + work(6 + n * (i - 1)) * bt(6,1) &
426 + work(7 + n * (i - 1)) * bt(7,1) &
427 + work(8 + n * (i - 1)) * bt(8,1) &
428 + work(9 + n * (i - 1)) * bt(9,1) &
429 + work(10 + n * (i - 1)) * bt(10,1) &
430 + work(11 + n * (i - 1)) * bt(11,1) &
431 + work(12 + n * (i - 1)) * bt(12,1)
432 end do
433
434 v(1) = work2(1) * ct(1, 1) &
435 + work2(2) * ct(2, 1) &
436 + work2(3) * ct(3, 1) &
437 + work2(4) * ct(4, 1) &
438 + work2(5) * ct(5, 1) &
439 + work2(6) * ct(6, 1) &
440 + work2(7) * ct(7, 1) &
441 + work2(8) * ct(8, 1) &
442 + work2(9) * ct(9, 1) &
443 + work2(10) * ct(10, 1) &
444 + work2(11) * ct(11, 1) &
445 + work2(12) * ct(12, 1)
446
447 end subroutine tnsr3d_el_1_12_cpu
448
459 !OCL SERIAL
460 subroutine tnsr3d_el_n_cpu(v, u, A, Bt, Ct, n)
461 integer, intent(in) :: n
462 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
463 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
464 real(kind=rp) :: work(n**3), work2(n**3), tmp
465 integer :: i, j, l, k
466 integer :: ii, jj, nn
467
468 nn = n**2
469
470 do j = 1, nn
471 do i = 1, n
472 ii = i + n * (j - 1)
473 tmp = 0.0_rp
474 do k = 1, n
475 tmp = tmp + a(i,k) * u(k + n * (j - 1))
476 end do
477 work(ii) = tmp
478 end do
479 end do
480
481 do i = 1, n
482 do j = 1, n
483 do l = 1, n
484 ii = l + n * (j - 1) + nn * (i - 1)
485 tmp = 0.0_rp
486 do k = 1, n
487 tmp = tmp + work(l + n * (k - 1) + nn * (i - 1)) * bt(k,j)
488 end do
489 work2(ii) = tmp
490 end do
491 end do
492 end do
493
494 do j = 1, n
495 do i = 1, nn
496 jj = i + nn * (j - 1)
497 tmp = 0.0_rp
498 do k = 1, n
499 tmp = tmp + work2(i + nn * (k - 1)) * ct(k, j)
500 end do
501 v(jj) = tmp
502 end do
503 end do
504
505 end subroutine tnsr3d_el_n_cpu
506
510 !OCL SERIAL
511 subroutine tnsr3d_el_n14_cpu(v, u, A, Bt, Ct)
512 integer, parameter :: n = 14
513 integer, parameter :: nn = n**2
514 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
515 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
516 real(kind=rp) :: work(n**3), work2(n**3)
517 integer :: i, j, l
518 integer :: ii, jj
519
520 do j = 1, nn
521 do i = 1, n
522 ii = i + n * (j - 1)
523 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
524 + a(i,2) * u(2 + n * (j - 1)) &
525 + a(i,3) * u(3 + n * (j - 1)) &
526 + a(i,4) * u(4 + n * (j - 1)) &
527 + a(i,5) * u(5 + n * (j - 1)) &
528 + a(i,6) * u(6 + n * (j - 1)) &
529 + a(i,7) * u(7 + n * (j - 1)) &
530 + a(i,8) * u(8 + n * (j - 1)) &
531 + a(i,9) * u(9 + n * (j - 1)) &
532 + a(i,10) * u(10 + n * (j - 1)) &
533 + a(i,11) * u(11 + n * (j - 1)) &
534 + a(i,12) * u(12 + n * (j - 1)) &
535 + a(i,13) * u(13 + n * (j - 1)) &
536 + a(i,14) * u(14 + n * (j - 1))
537 end do
538 end do
539
540 do i = 1, n
541 do j = 1, n
542 do l = 1, n
543 ii = l + n * (j - 1) + nn * (i - 1)
544 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
545 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
546 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
547 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
548 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
549 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
550 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
551 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
552 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j) &
553 + work(l + n * (10 - 1) + nn * (i - 1)) * bt(10,j) &
554 + work(l + n * (11 - 1) + nn * (i - 1)) * bt(11,j) &
555 + work(l + n * (12 - 1) + nn * (i - 1)) * bt(12,j) &
556 + work(l + n * (13 - 1) + nn * (i - 1)) * bt(13,j) &
557 + work(l + n * (14 - 1) + nn * (i - 1)) * bt(14,j)
558 end do
559 end do
560 end do
561
562 do j = 1, n
563 do i = 1, nn
564 jj = i + nn * (j - 1)
565 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
566 + work2(i + nn * (2 - 1)) * ct(2, j) &
567 + work2(i + nn * (3 - 1)) * ct(3, j) &
568 + work2(i + nn * (4 - 1)) * ct(4, j) &
569 + work2(i + nn * (5 - 1)) * ct(5, j) &
570 + work2(i + nn * (6 - 1)) * ct(6, j) &
571 + work2(i + nn * (7 - 1)) * ct(7, j) &
572 + work2(i + nn * (8 - 1)) * ct(8, j) &
573 + work2(i + nn * (9 - 1)) * ct(9, j) &
574 + work2(i + nn * (10 - 1)) * ct(10, j) &
575 + work2(i + nn * (11 - 1)) * ct(11, j) &
576 + work2(i + nn * (12 - 1)) * ct(12, j) &
577 + work2(i + nn * (13 - 1)) * ct(13, j) &
578 + work2(i + nn * (14 - 1)) * ct(14, j)
579 end do
580 end do
581
582 end subroutine tnsr3d_el_n14_cpu
583
587 !OCL SERIAL
588 subroutine tnsr3d_el_n13_cpu(v, u, A, Bt, Ct)
589 integer, parameter :: n = 13
590 integer, parameter :: nn = n**2
591 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
592 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
593 real(kind=rp) :: work(n**3), work2(n**3)
594 integer :: i, j, l
595 integer :: ii, jj
596
597 do j = 1, nn
598 do i = 1, n
599 ii = i + n * (j - 1)
600 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
601 + a(i,2) * u(2 + n * (j - 1)) &
602 + a(i,3) * u(3 + n * (j - 1)) &
603 + a(i,4) * u(4 + n * (j - 1)) &
604 + a(i,5) * u(5 + n * (j - 1)) &
605 + a(i,6) * u(6 + n * (j - 1)) &
606 + a(i,7) * u(7 + n * (j - 1)) &
607 + a(i,8) * u(8 + n * (j - 1)) &
608 + a(i,9) * u(9 + n * (j - 1)) &
609 + a(i,10) * u(10 + n * (j - 1)) &
610 + a(i,11) * u(11 + n * (j - 1)) &
611 + a(i,12) * u(12 + n * (j - 1)) &
612 + a(i,13) * u(13 + n * (j - 1))
613 end do
614 end do
615
616 do i = 1, n
617 do j = 1, n
618 do l = 1, n
619 ii = l + n * (j - 1) + nn * (i - 1)
620 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
621 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
622 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
623 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
624 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
625 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
626 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
627 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
628 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j) &
629 + work(l + n * (10 - 1) + nn * (i - 1)) * bt(10,j) &
630 + work(l + n * (11 - 1) + nn * (i - 1)) * bt(11,j) &
631 + work(l + n * (12 - 1) + nn * (i - 1)) * bt(12,j) &
632 + work(l + n * (13 - 1) + nn * (i - 1)) * bt(13,j)
633 end do
634 end do
635 end do
636
637 do j = 1, n
638 do i = 1, nn
639 jj = i + nn * (j - 1)
640 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
641 + work2(i + nn * (2 - 1)) * ct(2, j) &
642 + work2(i + nn * (3 - 1)) * ct(3, j) &
643 + work2(i + nn * (4 - 1)) * ct(4, j) &
644 + work2(i + nn * (5 - 1)) * ct(5, j) &
645 + work2(i + nn * (6 - 1)) * ct(6, j) &
646 + work2(i + nn * (7 - 1)) * ct(7, j) &
647 + work2(i + nn * (8 - 1)) * ct(8, j) &
648 + work2(i + nn * (9 - 1)) * ct(9, j) &
649 + work2(i + nn * (10 - 1)) * ct(10, j) &
650 + work2(i + nn * (11 - 1)) * ct(11, j) &
651 + work2(i + nn * (12 - 1)) * ct(12, j) &
652 + work2(i + nn * (13 - 1)) * ct(13, j)
653 end do
654 end do
655
656 end subroutine tnsr3d_el_n13_cpu
657
661 !OCL SERIAL
662 subroutine tnsr3d_el_n12_cpu(v, u, A, Bt, Ct)
663 integer, parameter :: n = 12
664 integer, parameter :: nn = n**2
665 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
666 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
667 real(kind=rp) :: work(n**3), work2(n**3)
668 integer :: i, j, l
669 integer :: ii, jj
670
671 do j = 1, nn
672 do i = 1, n
673 ii = i + n * (j - 1)
674 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
675 + a(i,2) * u(2 + n * (j - 1)) &
676 + a(i,3) * u(3 + n * (j - 1)) &
677 + a(i,4) * u(4 + n * (j - 1)) &
678 + a(i,5) * u(5 + n * (j - 1)) &
679 + a(i,6) * u(6 + n * (j - 1)) &
680 + a(i,7) * u(7 + n * (j - 1)) &
681 + a(i,8) * u(8 + n * (j - 1)) &
682 + a(i,9) * u(9 + n * (j - 1)) &
683 + a(i,10) * u(10 + n * (j - 1)) &
684 + a(i,11) * u(11 + n * (j - 1)) &
685 + a(i,12) * u(12 + n * (j - 1))
686 end do
687 end do
688
689 do i = 1, n
690 do j = 1, n
691 do l = 1, n
692 ii = l + n * (j - 1) + nn * (i - 1)
693 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
694 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
695 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
696 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
697 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
698 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
699 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
700 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
701 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j) &
702 + work(l + n * (10 - 1) + nn * (i - 1)) * bt(10,j) &
703 + work(l + n * (11 - 1) + nn * (i - 1)) * bt(11,j) &
704 + work(l + n * (12 - 1) + nn * (i - 1)) * bt(12,j)
705 end do
706 end do
707 end do
708
709 do j = 1, n
710 do i = 1, nn
711 jj = i + nn * (j - 1)
712 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
713 + work2(i + nn * (2 - 1)) * ct(2, j) &
714 + work2(i + nn * (3 - 1)) * ct(3, j) &
715 + work2(i + nn * (4 - 1)) * ct(4, j) &
716 + work2(i + nn * (5 - 1)) * ct(5, j) &
717 + work2(i + nn * (6 - 1)) * ct(6, j) &
718 + work2(i + nn * (7 - 1)) * ct(7, j) &
719 + work2(i + nn * (8 - 1)) * ct(8, j) &
720 + work2(i + nn * (9 - 1)) * ct(9, j) &
721 + work2(i + nn * (10 - 1)) * ct(10, j) &
722 + work2(i + nn * (11 - 1)) * ct(11, j) &
723 + work2(i + nn * (12 - 1)) * ct(12, j)
724 end do
725 end do
726
727 end subroutine tnsr3d_el_n12_cpu
728
732 !OCL SERIAL
733 subroutine tnsr3d_el_n11_cpu(v, u, A, Bt, Ct)
734 integer, parameter :: n = 11
735 integer, parameter :: nn = n**2
736 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
737 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
738 real(kind=rp) :: work(n**3), work2(n**3)
739 integer :: i, j, l
740 integer :: ii, jj
741
742 do j = 1, nn
743 do i = 1, n
744 ii = i + n * (j - 1)
745 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
746 + a(i,2) * u(2 + n * (j - 1)) &
747 + a(i,3) * u(3 + n * (j - 1)) &
748 + a(i,4) * u(4 + n * (j - 1)) &
749 + a(i,5) * u(5 + n * (j - 1)) &
750 + a(i,6) * u(6 + n * (j - 1)) &
751 + a(i,7) * u(7 + n * (j - 1)) &
752 + a(i,8) * u(8 + n * (j - 1)) &
753 + a(i,9) * u(9 + n * (j - 1)) &
754 + a(i,10) * u(10 + n * (j - 1)) &
755 + a(i,11) * u(11 + n * (j - 1))
756 end do
757 end do
758
759 do i = 1, n
760 do j = 1, n
761 do l = 1, n
762 ii = l + n * (j - 1) + nn * (i - 1)
763 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
764 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
765 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
766 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
767 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
768 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
769 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
770 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
771 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j) &
772 + work(l + n * (10 - 1) + nn * (i - 1)) * bt(10,j) &
773 + work(l + n * (11 - 1) + nn * (i - 1)) * bt(11,j)
774 end do
775 end do
776 end do
777
778 do j = 1, n
779 do i = 1, nn
780 jj = i + nn * (j - 1)
781 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
782 + work2(i + nn * (2 - 1)) * ct(2, j) &
783 + work2(i + nn * (3 - 1)) * ct(3, j) &
784 + work2(i + nn * (4 - 1)) * ct(4, j) &
785 + work2(i + nn * (5 - 1)) * ct(5, j) &
786 + work2(i + nn * (6 - 1)) * ct(6, j) &
787 + work2(i + nn * (7 - 1)) * ct(7, j) &
788 + work2(i + nn * (8 - 1)) * ct(8, j) &
789 + work2(i + nn * (9 - 1)) * ct(9, j) &
790 + work2(i + nn * (10 - 1)) * ct(10, j) &
791 + work2(i + nn * (11 - 1)) * ct(11, j)
792 end do
793 end do
794
795 end subroutine tnsr3d_el_n11_cpu
796
800 !OCL SERIAL
801 subroutine tnsr3d_el_n10_cpu(v, u, A, Bt, Ct)
802 integer, parameter :: n = 10
803 integer, parameter :: nn = n**2
804 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
805 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
806 real(kind=rp) :: work(n**3), work2(n**3)
807 integer :: i, j, l
808 integer :: ii, jj
809
810 do j = 1, nn
811 do i = 1, n
812 ii = i + n * (j - 1)
813 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
814 + a(i,2) * u(2 + n * (j - 1)) &
815 + a(i,3) * u(3 + n * (j - 1)) &
816 + a(i,4) * u(4 + n * (j - 1)) &
817 + a(i,5) * u(5 + n * (j - 1)) &
818 + a(i,6) * u(6 + n * (j - 1)) &
819 + a(i,7) * u(7 + n * (j - 1)) &
820 + a(i,8) * u(8 + n * (j - 1)) &
821 + a(i,9) * u(9 + n * (j - 1)) &
822 + a(i,10) * u(10 + n * (j - 1))
823 end do
824 end do
825
826 do i = 1, n
827 do j = 1, n
828 do l = 1, n
829 ii = l + n * (j - 1) + nn * (i - 1)
830 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
831 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
832 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
833 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
834 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
835 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
836 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
837 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
838 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j) &
839 + work(l + n * (10 - 1) + nn * (i - 1)) * bt(10,j)
840 end do
841 end do
842 end do
843
844 do j = 1, n
845 do i = 1, nn
846 jj = i + nn * (j - 1)
847 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
848 + work2(i + nn * (2 - 1)) * ct(2, j) &
849 + work2(i + nn * (3 - 1)) * ct(3, j) &
850 + work2(i + nn * (4 - 1)) * ct(4, j) &
851 + work2(i + nn * (5 - 1)) * ct(5, j) &
852 + work2(i + nn * (6 - 1)) * ct(6, j) &
853 + work2(i + nn * (7 - 1)) * ct(7, j) &
854 + work2(i + nn * (8 - 1)) * ct(8, j) &
855 + work2(i + nn * (9 - 1)) * ct(9, j) &
856 + work2(i + nn * (10 - 1)) * ct(10, j)
857 end do
858 end do
859
860 end subroutine tnsr3d_el_n10_cpu
861
865 !OCL SERIAL
866 subroutine tnsr3d_el_n9_cpu(v, u, A, Bt, Ct)
867 integer, parameter :: n = 9
868 integer, parameter :: nn = n**2
869 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
870 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
871 real(kind=rp) :: work(n**3), work2(n**3)
872 integer :: i, j, l
873 integer :: ii, jj
874
875 do j = 1, nn
876 do i = 1, n
877 ii = i + n * (j - 1)
878 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
879 + a(i,2) * u(2 + n * (j - 1)) &
880 + a(i,3) * u(3 + n * (j - 1)) &
881 + a(i,4) * u(4 + n * (j - 1)) &
882 + a(i,5) * u(5 + n * (j - 1)) &
883 + a(i,6) * u(6 + n * (j - 1)) &
884 + a(i,7) * u(7 + n * (j - 1)) &
885 + a(i,8) * u(8 + n * (j - 1)) &
886 + a(i,9) * u(9 + n * (j - 1))
887 end do
888 end do
889
890 do i = 1, n
891 do j = 1, n
892 do l = 1, n
893 ii = l + n * (j - 1) + nn * (i - 1)
894 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
895 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
896 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
897 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
898 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
899 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
900 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
901 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j) &
902 + work(l + n * (9 - 1) + nn * (i - 1)) * bt(9,j)
903 end do
904 end do
905 end do
906
907 do j = 1, n
908 do i = 1, nn
909 jj = i + nn * (j - 1)
910 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
911 + work2(i + nn * (2 - 1)) * ct(2, j) &
912 + work2(i + nn * (3 - 1)) * ct(3, j) &
913 + work2(i + nn * (4 - 1)) * ct(4, j) &
914 + work2(i + nn * (5 - 1)) * ct(5, j) &
915 + work2(i + nn * (6 - 1)) * ct(6, j) &
916 + work2(i + nn * (7 - 1)) * ct(7, j) &
917 + work2(i + nn * (8 - 1)) * ct(8, j) &
918 + work2(i + nn * (9 - 1)) * ct(9, j)
919 end do
920 end do
921
922 end subroutine tnsr3d_el_n9_cpu
923
927 !OCL SERIAL
928 subroutine tnsr3d_el_n8_cpu(v, u, A, Bt, Ct)
929 integer, parameter :: n = 8
930 integer, parameter :: nn = n**2
931 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
932 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
933 real(kind=rp) :: work(n**3), work2(n**3)
934 integer :: i, j, l
935 integer :: ii, jj
936
937 do j = 1, nn
938 do i = 1, n
939 ii = i + n * (j - 1)
940 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
941 + a(i,2) * u(2 + n * (j - 1)) &
942 + a(i,3) * u(3 + n * (j - 1)) &
943 + a(i,4) * u(4 + n * (j - 1)) &
944 + a(i,5) * u(5 + n * (j - 1)) &
945 + a(i,6) * u(6 + n * (j - 1)) &
946 + a(i,7) * u(7 + n * (j - 1)) &
947 + a(i,8) * u(8 + n * (j - 1))
948 end do
949 end do
950
951 do i = 1, n
952 do j = 1, n
953 do l = 1, n
954 ii = l + n * (j - 1) + nn * (i - 1)
955 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
956 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
957 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
958 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
959 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
960 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
961 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j) &
962 + work(l + n * (8 - 1) + nn * (i - 1)) * bt(8,j)
963 end do
964 end do
965 end do
966
967 do j = 1, n
968 do i = 1, nn
969 jj = i + nn * (j - 1)
970 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
971 + work2(i + nn * (2 - 1)) * ct(2, j) &
972 + work2(i + nn * (3 - 1)) * ct(3, j) &
973 + work2(i + nn * (4 - 1)) * ct(4, j) &
974 + work2(i + nn * (5 - 1)) * ct(5, j) &
975 + work2(i + nn * (6 - 1)) * ct(6, j) &
976 + work2(i + nn * (7 - 1)) * ct(7, j) &
977 + work2(i + nn * (8 - 1)) * ct(8, j)
978 end do
979 end do
980
981 end subroutine tnsr3d_el_n8_cpu
982
986 !OCL SERIAL
987 subroutine tnsr3d_el_n7_cpu(v, u, A, Bt, Ct)
988 integer, parameter :: n = 7
989 integer, parameter :: nn = n**2
990 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
991 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
992 real(kind=rp) :: work(n**3), work2(n**3)
993 integer :: i, j, l
994 integer :: ii, jj
995
996 do j = 1, nn
997 do i = 1, n
998 ii = i + n * (j - 1)
999 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1000 + a(i,2) * u(2 + n * (j - 1)) &
1001 + a(i,3) * u(3 + n * (j - 1)) &
1002 + a(i,4) * u(4 + n * (j - 1)) &
1003 + a(i,5) * u(5 + n * (j - 1)) &
1004 + a(i,6) * u(6 + n * (j - 1)) &
1005 + a(i,7) * u(7 + n * (j - 1))
1006 end do
1007 end do
1008
1009 do i = 1, n
1010 do j = 1, n
1011 do l = 1, n
1012 ii = l + n * (j - 1) + nn * (i - 1)
1013 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1014 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
1015 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
1016 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
1017 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
1018 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j) &
1019 + work(l + n * (7 - 1) + nn * (i - 1)) * bt(7,j)
1020 end do
1021 end do
1022 end do
1023
1024 do j = 1, n
1025 do i = 1, nn
1026 jj = i + nn * (j - 1)
1027 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1028 + work2(i + nn * (2 - 1)) * ct(2, j) &
1029 + work2(i + nn * (3 - 1)) * ct(3, j) &
1030 + work2(i + nn * (4 - 1)) * ct(4, j) &
1031 + work2(i + nn * (5 - 1)) * ct(5, j) &
1032 + work2(i + nn * (6 - 1)) * ct(6, j) &
1033 + work2(i + nn * (7 - 1)) * ct(7, j)
1034 end do
1035 end do
1036
1037 end subroutine tnsr3d_el_n7_cpu
1038
1042 !OCL SERIAL
1043 subroutine tnsr3d_el_n6_cpu(v, u, A, Bt, Ct)
1044 integer, parameter :: n = 6
1045 integer, parameter :: nn = n**2
1046 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
1047 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
1048 real(kind=rp) :: work(n**3), work2(n**3)
1049 integer :: i, j, l
1050 integer :: ii, jj
1051
1052 do j = 1, nn
1053 do i = 1, n
1054 ii = i + n * (j - 1)
1055 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1056 + a(i,2) * u(2 + n * (j - 1)) &
1057 + a(i,3) * u(3 + n * (j - 1)) &
1058 + a(i,4) * u(4 + n * (j - 1)) &
1059 + a(i,5) * u(5 + n * (j - 1)) &
1060 + a(i,6) * u(6 + n * (j - 1))
1061 end do
1062 end do
1063
1064 do i = 1, n
1065 do j = 1, n
1066 do l = 1, n
1067 ii = l + n * (j - 1) + nn * (i - 1)
1068 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1069 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
1070 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
1071 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
1072 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j) &
1073 + work(l + n * (6 - 1) + nn * (i - 1)) * bt(6,j)
1074 end do
1075 end do
1076 end do
1077
1078 do j = 1, n
1079 do i = 1, nn
1080 jj = i + nn * (j - 1)
1081 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1082 + work2(i + nn * (2 - 1)) * ct(2, j) &
1083 + work2(i + nn * (3 - 1)) * ct(3, j) &
1084 + work2(i + nn * (4 - 1)) * ct(4, j) &
1085 + work2(i + nn * (5 - 1)) * ct(5, j) &
1086 + work2(i + nn * (6 - 1)) * ct(6, j)
1087 end do
1088 end do
1089
1090 end subroutine tnsr3d_el_n6_cpu
1091
1095 !OCL SERIAL
1096 subroutine tnsr3d_el_n5_cpu(v, u, A, Bt, Ct)
1097 integer, parameter :: n = 5
1098 integer, parameter :: nn = n**2
1099 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
1100 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
1101 real(kind=rp) :: work(n**3), work2(n**3)
1102 integer :: i, j, l
1103 integer :: ii, jj
1104
1105 do j = 1, nn
1106 do i = 1, n
1107 ii = i + n * (j - 1)
1108 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1109 + a(i,2) * u(2 + n * (j - 1)) &
1110 + a(i,3) * u(3 + n * (j - 1)) &
1111 + a(i,4) * u(4 + n * (j - 1)) &
1112 + a(i,5) * u(5 + n * (j - 1))
1113 end do
1114 end do
1115
1116 do i = 1, n
1117 do j = 1, n
1118 do l = 1, n
1119 ii = l + n * (j - 1) + nn * (i - 1)
1120 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1121 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
1122 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
1123 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j) &
1124 + work(l + n * (5 - 1) + nn * (i - 1)) * bt(5,j)
1125 end do
1126 end do
1127 end do
1128
1129 do j = 1, n
1130 do i = 1, nn
1131 jj = i + nn * (j - 1)
1132 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1133 + work2(i + nn * (2 - 1)) * ct(2, j) &
1134 + work2(i + nn * (3 - 1)) * ct(3, j) &
1135 + work2(i + nn * (4 - 1)) * ct(4, j) &
1136 + work2(i + nn * (5 - 1)) * ct(5, j)
1137 end do
1138 end do
1139
1140 end subroutine tnsr3d_el_n5_cpu
1141
1145 !OCL SERIAL
1146 subroutine tnsr3d_el_n4_cpu(v, u, A, Bt, Ct)
1147 integer, parameter :: n = 4
1148 integer, parameter :: nn = n**2
1149 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
1150 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
1151 real(kind=rp) :: work(n**3), work2(n**3)
1152 integer :: i, j, l
1153 integer :: ii, jj
1154
1155 do j = 1, nn
1156 do i = 1, n
1157 ii = i + n * (j - 1)
1158 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1159 + a(i,2) * u(2 + n * (j - 1)) &
1160 + a(i,3) * u(3 + n * (j - 1)) &
1161 + a(i,4) * u(4 + n * (j - 1))
1162 end do
1163 end do
1164
1165 do i = 1, n
1166 do j = 1, n
1167 do l = 1, n
1168 ii = l + n * (j - 1) + nn * (i - 1)
1169 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1170 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
1171 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j) &
1172 + work(l + n * (4 - 1) + nn * (i - 1)) * bt(4,j)
1173 end do
1174 end do
1175 end do
1176
1177 do j = 1, n
1178 do i = 1, nn
1179 jj = i + nn * (j - 1)
1180 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1181 + work2(i + nn * (2 - 1)) * ct(2, j) &
1182 + work2(i + nn * (3 - 1)) * ct(3, j) &
1183 + work2(i + nn * (4 - 1)) * ct(4, j)
1184 end do
1185 end do
1186
1187 end subroutine tnsr3d_el_n4_cpu
1188
1192 !OCL SERIAL
1193 subroutine tnsr3d_el_n3_cpu(v, u, A, Bt, Ct)
1194 integer, parameter :: n = 3
1195 integer, parameter :: nn = n**2
1196 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
1197 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
1198 real(kind=rp) :: work(n**3), work2(n**3)
1199 integer :: i, j, l
1200 integer :: ii, jj
1201
1202 do j = 1, nn
1203 do i = 1, n
1204 ii = i + n * (j - 1)
1205 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1206 + a(i,2) * u(2 + n * (j - 1)) &
1207 + a(i,3) * u(3 + n * (j - 1))
1208 end do
1209 end do
1210
1211 do i = 1, n
1212 do j = 1, n
1213 do l = 1, n
1214 ii = l + n * (j - 1) + nn * (i - 1)
1215 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1216 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j) &
1217 + work(l + n * (3 - 1) + nn * (i - 1)) * bt(3,j)
1218 end do
1219 end do
1220 end do
1221
1222 do j = 1, n
1223 do i = 1, nn
1224 jj = i + nn * (j - 1)
1225 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1226 + work2(i + nn * (2 - 1)) * ct(2, j) &
1227 + work2(i + nn * (3 - 1)) * ct(3, j)
1228 end do
1229 end do
1230
1231 end subroutine tnsr3d_el_n3_cpu
1232
1236 !OCL SERIAL
1237 subroutine tnsr3d_el_n2_cpu(v, u, A, Bt, Ct)
1238 integer, parameter :: n = 2
1239 integer, parameter :: nn = n**2
1240 real(kind=rp), intent(inout) :: v(n*n*n), u(n*n*n)
1241 real(kind=rp), intent(inout) :: a(n,n), bt(n,n), ct(n,n)
1242 real(kind=rp) :: work(n**3), work2(n**3)
1243 integer :: i, j, l
1244 integer :: ii, jj
1245
1246 do j = 1, nn
1247 do i = 1, n
1248 ii = i + n * (j - 1)
1249 work(ii) = a(i,1) * u(1 + n * (j - 1)) &
1250 + a(i,2) * u(2 + n * (j - 1))
1251 end do
1252 end do
1253
1254 do i = 1, n
1255 do j = 1, n
1256 do l = 1, n
1257 ii = l + n * (j - 1) + nn * (i - 1)
1258 work2(ii) = work(l + n * (1 - 1) + nn * (i - 1)) * bt(1,j) &
1259 + work(l + n * (2 - 1) + nn * (i - 1)) * bt(2,j)
1260 end do
1261 end do
1262 end do
1263
1264 do j = 1, n
1265 do i = 1, nn
1266 jj = i + nn * (j - 1)
1267 v(jj) = work2(i + nn * (1 - 1)) * ct(1, j) &
1268 + work2(i + nn * (2 - 1)) * ct(2, j)
1269 end do
1270 end do
1271
1272 end subroutine tnsr3d_el_n2_cpu
1273
1291 subroutine tnsr3d_cpu(v, nv, u, nu, A, Bt, Ct, nelv)
1292 integer, intent(in) :: nv, nu, nelv
1293 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1294 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1295 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1296
1297 if (nu .eq. 2 .and. nv .eq. 4) then
1298 call tnsr3d_nu2nv4_cpu(v, u, a, bt, ct, nelv)
1299 else if (nu .eq. 4) then
1300 call tnsr3d_nu4_cpu(v, nv, u, a, bt, ct, nelv)
1301 else if (nu .eq. 8) then
1302 call tnsr3d_nu8_cpu(v, nv, u, a, bt, ct, nelv)
1303 else if (nu .eq. 12) then
1304 call tnsr3d_nu12_cpu(v, nv, u, a, bt, ct, nelv)
1305 else
1306 call tnsr3d_nvnu_cpu(v, nv, u, nu, a, bt, ct, nelv)
1307 end if
1308
1309 end subroutine tnsr3d_cpu
1310
1323 subroutine tnsr3d_nvnu_cpu(v, nv, u, nu, A, Bt, Ct, nelv)
1324 integer, intent(in) :: nv, nu, nelv
1325 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1326 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1327 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1328 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2), c
1329 integer :: ie, i, j, k, l, ii, jj, kk
1330 integer :: nunu, nvnu, nvnv
1331
1332 nvnu = nv * nu
1333 nunu = nu * nu
1334 nvnv = nv * nv
1335
1336 !$omp parallel do private(ie, i, j, k, l, ii, jj, kk, c, work, work2)
1337 do ie = 1, nelv
1338 do j = 1, nunu
1339 ii = nv * (j - 1)
1340 kk = nu * (j - 1)
1341 c = u(1 + kk, ie)
1342 do i = 1, nv
1343 work(i + ii) = a(i,1) * c
1344 end do
1345 do k = 2, nu
1346 c = u(k + kk, ie)
1347 !OCL NORECURRENCE, NOVREC, NOALIAS
1348 !DIR$ IVDEP
1349 !GCC$ ivdep
1350 do i = 1, nv
1351 work(i + ii) = work(i + ii) + a(i,k) * c
1352 end do
1353 end do
1354 end do
1355
1356 do i = 1, nu
1357 do j = 1, nv
1358 ii = nv * (j - 1) + nvnv * (i - 1)
1359 jj = nvnu * (i - 1)
1360 c = bt(1,j)
1361 do l = 1, nv
1362 work2(l + ii) = work(l + jj) * c
1363 end do
1364 do k = 2, nu
1365 c = bt(k,j)
1366 kk = nv * (k - 1) + nvnu * (i - 1)
1367 !OCL NORECURRENCE, NOVREC, NOALIAS
1368 !DIR$ IVDEP
1369 !GCC$ ivdep
1370 do l = 1, nv
1371 work2(l + ii) = work2(l + ii) + work(l + kk) * c
1372 end do
1373 end do
1374 end do
1375 end do
1376
1377 do j = 1, nv
1378 jj = nvnv * (j - 1)
1379 c = ct(1,j)
1380 do i = 1, nvnv
1381 v(i + jj, ie) = work2(i) * c
1382 end do
1383 do k = 2, nu
1384 c = ct(k,j)
1385 ii = nvnv * (k - 1)
1386 !OCL NORECURRENCE, NOVREC, NOALIAS
1387 !DIR$ IVDEP
1388 !GCC$ ivdep
1389 do i = 1, nvnv
1390 v(i + jj, ie) = v(i + jj, ie) + work2(i + ii) * c
1391 end do
1392 end do
1393 end do
1394 end do
1395 !$omp end parallel do
1396
1397 end subroutine tnsr3d_nvnu_cpu
1398
1402 subroutine tnsr3d_nu2nv4_cpu(v, u, A, Bt, Ct, nelv)
1403 integer, parameter :: nu = 2
1404 integer, parameter :: nv = 4
1405 integer, parameter :: nunu = 4
1406 integer, parameter :: nvnu = 8
1407 integer, parameter :: nvnv = 16
1408 integer, intent(in) :: nelv
1409 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1410 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1411 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1412 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2), tmp
1413 integer :: ie, i, j, k, l, ii, jj
1414
1415 !$omp parallel do private(ie, i, j, k, l, ii, jj, tmp, work, work2)
1416 do ie = 1, nelv
1417 do j = 1, nunu
1418 do i = 1, nv
1419 ii = i + nv * (j - 1)
1420 work(ii) = a(i,1) * u(1 + nu * (j - 1), ie) &
1421 + a(i,2) * u(2 + nu * (j - 1), ie)
1422 end do
1423 end do
1424
1425 do i = 1, nu
1426 do j = 1, nv
1427 do l = 1, nv
1428 ii = l + nv * (j - 1) + nvnv * (i - 1)
1429 tmp = 0.0_rp
1430 do k = 1, nu
1431 jj = l + nv * (k - 1) + nvnu * (i - 1)
1432 tmp = tmp + work(jj) * bt(k,j)
1433 end do
1434 work2(ii) = tmp
1435 end do
1436 end do
1437 end do
1438
1439 do j = 1, nv
1440 do i = 1, nvnv
1441 jj = i + nvnv * (j - 1)
1442 v(jj, ie) = work2(i + nvnv * (1 - 1)) * ct(1, j) &
1443 + work2(i + nvnv * (2 - 1)) * ct(2, j)
1444 end do
1445 end do
1446 end do
1447 !$omp end parallel do
1448
1449 end subroutine tnsr3d_nu2nv4_cpu
1450
1454 subroutine tnsr3d_nu4_cpu(v, nv, u, A, Bt, Ct, nelv)
1455 integer, parameter :: nu = 4
1456 integer, parameter :: nunu = 16
1457 integer, intent(in) :: nv, nelv
1458 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1459 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1460 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1461 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2), tmp
1462 integer :: ie, i, j, k, l, ii, jj
1463 integer :: nvnu, nvnv
1464
1465 nvnu = nv * nu
1466 nvnv = nv * nv
1467
1468 !$omp parallel do private(ie, i, j, k, l, ii, jj, tmp, work, work2)
1469 do ie = 1, nelv
1470 do j = 1, nunu
1471 do i = 1, nv
1472 ii = i + nv * (j - 1)
1473 work(ii) = a(i,1) * u(1 + nu * (j - 1), ie) &
1474 + a(i,2) * u(2 + nu * (j - 1), ie) &
1475 + a(i,3) * u(3 + nu * (j - 1), ie) &
1476 + a(i,4) * u(4 + nu * (j - 1), ie)
1477 end do
1478 end do
1479
1480 do i = 1, nu
1481 do j = 1, nv
1482 do l = 1, nv
1483 ii = l + nv * (j - 1) + nvnv * (i - 1)
1484 tmp = 0.0_rp
1485 do k = 1, nu
1486 jj = l + nv * (k - 1) + nvnu * (i - 1)
1487 tmp = tmp + work(jj) * bt(k,j)
1488 end do
1489 work2(ii) = tmp
1490 end do
1491 end do
1492 end do
1493
1494 do j = 1, nv
1495 do i = 1, nvnv
1496 jj = i + nvnv * (j - 1)
1497 v(jj, ie) = work2(i + nvnv * (1 - 1)) * ct(1, j) &
1498 + work2(i + nvnv * (2 - 1)) * ct(2, j) &
1499 + work2(i + nvnv * (3 - 1)) * ct(3, j) &
1500 + work2(i + nvnv * (4 - 1)) * ct(4, j)
1501 end do
1502 end do
1503 end do
1504 !$omp end parallel do
1505
1506 end subroutine tnsr3d_nu4_cpu
1507
1515 subroutine tnsr3d_nu8_cpu(v, nv, u, A, Bt, Ct, nelv)
1516 integer, parameter :: nu = 8
1517 integer, parameter :: nunu = 64
1518 integer, intent(in) :: nv, nelv
1519 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1520 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1521 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1522 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2)
1523 integer :: ie, i, j, l, ii, jj
1524 integer :: nvnu, nvnv
1525
1526 nvnu = nv * nu
1527 nvnv = nv * nv
1528
1529 !$omp parallel do private(ie, i, j, l, ii, jj, work, work2)
1530 do ie = 1, nelv
1531 do j = 1, nunu
1532 do i = 1, nv
1533 ii = i + nv * (j - 1)
1534 work(ii) = a(i,1) * u(1 + nu * (j - 1), ie) &
1535 + a(i,2) * u(2 + nu * (j - 1), ie) &
1536 + a(i,3) * u(3 + nu * (j - 1), ie) &
1537 + a(i,4) * u(4 + nu * (j - 1), ie) &
1538 + a(i,5) * u(5 + nu * (j - 1), ie) &
1539 + a(i,6) * u(6 + nu * (j - 1), ie) &
1540 + a(i,7) * u(7 + nu * (j - 1), ie) &
1541 + a(i,8) * u(8 + nu * (j - 1), ie)
1542 end do
1543 end do
1544
1545 do i = 1, nu
1546 do j = 1, nv
1547 ii = nv * (j - 1) + nvnv * (i - 1)
1548 jj = nvnu * (i - 1)
1549 do l = 1, nv
1550 work2(l + ii) = work(l + jj) * bt(1,j) &
1551 + work(l + nv + jj) * bt(2,j) &
1552 + work(l + 2 * nv + jj) * bt(3,j) &
1553 + work(l + 3 * nv + jj) * bt(4,j) &
1554 + work(l + 4 * nv + jj) * bt(5,j) &
1555 + work(l + 5 * nv + jj) * bt(6,j) &
1556 + work(l + 6 * nv + jj) * bt(7,j) &
1557 + work(l + 7 * nv + jj) * bt(8,j)
1558 end do
1559 end do
1560 end do
1561
1562 do j = 1, nv
1563 do i = 1, nvnv
1564 jj = i + nvnv * (j - 1)
1565 v(jj, ie) = work2(i + nvnv * (1 - 1)) * ct(1, j) &
1566 + work2(i + nvnv * (2 - 1)) * ct(2, j) &
1567 + work2(i + nvnv * (3 - 1)) * ct(3, j) &
1568 + work2(i + nvnv * (4 - 1)) * ct(4, j) &
1569 + work2(i + nvnv * (5 - 1)) * ct(5, j) &
1570 + work2(i + nvnv * (6 - 1)) * ct(6, j) &
1571 + work2(i + nvnv * (7 - 1)) * ct(7, j) &
1572 + work2(i + nvnv * (8 - 1)) * ct(8, j)
1573 end do
1574 end do
1575 end do
1576 !$omp end parallel do
1577
1578 end subroutine tnsr3d_nu8_cpu
1579
1586 subroutine tnsr3d_nu12_cpu(v, nv, u, A, Bt, Ct, nelv)
1587 integer, parameter :: nu = 12
1588 integer, parameter :: nunu = 144
1589 integer, intent(in) :: nv, nelv
1590 real(kind=rp), intent(inout) :: v(nv*nv*nv, nelv)
1591 real(kind=rp), intent(in) :: u(nu*nu*nu, nelv)
1592 real(kind=rp), intent(in) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1593 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2)
1594 integer :: ie, i, j, l, ii, jj
1595 integer :: nvnu, nvnv
1596
1597 nvnu = nv * nu
1598 nvnv = nv * nv
1599
1600 !$omp parallel do private(ie, i, j, l, ii, jj, work, work2)
1601 do ie = 1, nelv
1602 do j = 1, nunu
1603 do i = 1, nv
1604 ii = i + nv * (j - 1)
1605 work(ii) = a(i,1) * u(1 + nu * (j - 1), ie) &
1606 + a(i,2) * u(2 + nu * (j - 1), ie) &
1607 + a(i,3) * u(3 + nu * (j - 1), ie) &
1608 + a(i,4) * u(4 + nu * (j - 1), ie) &
1609 + a(i,5) * u(5 + nu * (j - 1), ie) &
1610 + a(i,6) * u(6 + nu * (j - 1), ie) &
1611 + a(i,7) * u(7 + nu * (j - 1), ie) &
1612 + a(i,8) * u(8 + nu * (j - 1), ie) &
1613 + a(i,9) * u(9 + nu * (j - 1), ie) &
1614 + a(i,10) * u(10 + nu * (j - 1), ie) &
1615 + a(i,11) * u(11 + nu * (j - 1), ie) &
1616 + a(i,12) * u(12 + nu * (j - 1), ie)
1617 end do
1618 end do
1619
1620 do i = 1, nu
1621 do j = 1, nv
1622 ii = nv * (j - 1) + nvnv * (i - 1)
1623 jj = nvnu * (i - 1)
1624 do l = 1, nv
1625 work2(l + ii) = work(l + jj) * bt(1,j) &
1626 + work(l + nv + jj) * bt(2,j) &
1627 + work(l + 2 * nv + jj) * bt(3,j) &
1628 + work(l + 3 * nv + jj) * bt(4,j) &
1629 + work(l + 4 * nv + jj) * bt(5,j) &
1630 + work(l + 5 * nv + jj) * bt(6,j) &
1631 + work(l + 6 * nv + jj) * bt(7,j) &
1632 + work(l + 7 * nv + jj) * bt(8,j) &
1633 + work(l + 8 * nv + jj) * bt(9,j) &
1634 + work(l + 9 * nv + jj) * bt(10,j) &
1635 + work(l + 10 * nv + jj) * bt(11,j) &
1636 + work(l + 11 * nv + jj) * bt(12,j)
1637 end do
1638 end do
1639 end do
1640
1641 do j = 1, nv
1642 do i = 1, nvnv
1643 jj = i + nvnv * (j - 1)
1644 v(jj, ie) = work2(i + nvnv * (1 - 1)) * ct(1, j) &
1645 + work2(i + nvnv * (2 - 1)) * ct(2, j) &
1646 + work2(i + nvnv * (3 - 1)) * ct(3, j) &
1647 + work2(i + nvnv * (4 - 1)) * ct(4, j) &
1648 + work2(i + nvnv * (5 - 1)) * ct(5, j) &
1649 + work2(i + nvnv * (6 - 1)) * ct(6, j) &
1650 + work2(i + nvnv * (7 - 1)) * ct(7, j) &
1651 + work2(i + nvnv * (8 - 1)) * ct(8, j) &
1652 + work2(i + nvnv * (9 - 1)) * ct(9, j) &
1653 + work2(i + nvnv * (10 - 1)) * ct(10, j) &
1654 + work2(i + nvnv * (11 - 1)) * ct(11, j) &
1655 + work2(i + nvnv * (12 - 1)) * ct(12, j)
1656 end do
1657 end do
1658 end do
1659 !$omp end parallel do
1660
1661 end subroutine tnsr3d_nu12_cpu
1662
1675 subroutine tnsr1_3d_cpu(v, nv, nu, A, Bt, Ct, nelv)
1676 integer, intent(in) :: nv, nu, nelv
1677 real(kind=rp), intent(inout) :: v(nv*nv*nv*nelv)
1678 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1679
1680 if (nu .eq. 4 .and. nv .eq. 2) then
1681 call tnsr1_3d_nu4nv2_cpu(v, a, bt, ct, nelv)
1682 else
1683 call tnsr1_3d_nvnu_cpu(v, nv, nu, a, bt, ct, nelv)
1684 end if
1685
1686 end subroutine tnsr1_3d_cpu
1687
1692 subroutine tnsr1_3d_nvnu_cpu(v, nv, nu, A, Bt, Ct, nelv)
1693 integer, intent(in) :: nv, nu, nelv
1694 real(kind=rp), intent(inout) :: v(nv*nv*nv*nelv)
1695 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1696 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2)
1697 integer :: e, e0, ee, es, iu, iv, nu3, nv3
1698 integer :: i, j, k, l, ii, jj, kk
1699 integer :: nunu, nvnu, nvnv
1700 real(kind=rp) :: tmp
1701
1702 nvnu = nv * nu
1703 nunu = nu * nu
1704 nvnv = nv * nv
1705
1706 e0 = 1
1707 es = 1
1708 ee = nelv
1709
1710 if (nv .gt. nu) then
1711 e0 = nelv
1712 es = -1
1713 ee = 1
1714 end if
1715
1716 nu3 = nu**3
1717 nv3 = nv**3
1718
1719 !$omp parallel do private(e,iu,iv,i,j,k,l,ii,jj,kk,work,work2,tmp)
1720 do e = e0, ee, es
1721 iu = (e-1)*nu3
1722 iv = (e-1)*nv3
1723
1724 do j = 1, nunu
1725 do i = 1, nv
1726 ii = i + nv * (j - 1)
1727 tmp = 0.0_rp
1728 do k = 1, nu
1729 kk = k + nu * (j - 1) + iu
1730 tmp = tmp + a(i,k) * v(kk)
1731 end do
1732 work(ii) = tmp
1733 end do
1734 end do
1735
1736 do i = 1, nu
1737 do j = 1, nv
1738 do l = 1, nv
1739 ii = l + nv * (j - 1) + nvnv * (i - 1)
1740 tmp = 0.0_rp
1741 do k = 1, nu
1742 jj = l + nv * (k - 1) + nvnu * (i - 1)
1743 tmp = tmp + work(jj) * bt(k,j)
1744 end do
1745 work2(ii) = tmp
1746 end do
1747 end do
1748 end do
1749
1750 do j = 1, nv
1751 do i = 1, nvnv
1752 jj = i + nvnv * (j - 1) + iv
1753 tmp = 0.0_rp
1754 do k = 1, nu
1755 ii = i + nvnv * (k - 1)
1756 tmp = tmp + work2(ii) * ct(k, j)
1757 end do
1758 v(jj) = tmp
1759 end do
1760 end do
1761 end do
1762 !$omp end parallel do
1763 end subroutine tnsr1_3d_nvnu_cpu
1764
1767 subroutine tnsr1_3d_nu4nv2_cpu(v, A, Bt, Ct, nelv)
1768 integer, parameter :: nu = 4
1769 integer, parameter :: nv = 2
1770 integer, parameter :: nunu = 16
1771 integer, parameter :: nvnu = 8
1772 integer, parameter :: nvnv = 4
1773 integer, parameter :: nununu = 64
1774 integer, parameter :: nvnvnv = 8
1775 integer, intent(in) :: nelv
1776 real(kind=rp), intent(inout) :: v(nv*nv*nv*nelv)
1777 real(kind=rp), intent(inout) :: a(nv, nu), bt(nu, nv), ct(nu, nv)
1778 real(kind=rp) :: work(nu**2*nv), work2(nu*nv**2)
1779 integer :: e, iu, iv
1780 integer :: i, j, k, l, ii, jj
1781 real(kind=rp) :: tmp
1782
1783 !$omp parallel do private(e,iu,iv,i,j,k,l,ii,jj, work, work2, tmp)
1784 do e = 1, nelv
1785 iu = (e-1)*nununu
1786 iv = (e-1)*nvnvnv
1787
1788 do j = 1, nunu
1789 do i = 1, nv
1790 ii = i + nv * (j - 1)
1791 work(ii) = a(i,1) * v(1 + nu * (j - 1) + iu) &
1792 + a(i,2) * v(2 + nu * (j - 1) + iu) &
1793 + a(i,3) * v(3 + nu * (j - 1) + iu) &
1794 + a(i,4) * v(4 + nu * (j - 1) + iu)
1795 end do
1796 end do
1797
1798 do i = 1, nu
1799 do j = 1, nv
1800 do l = 1, nv
1801 ii = l + nv * (j - 1) + nvnv * (i - 1)
1802 tmp = 0.0_rp
1803 do k = 1, nu
1804 jj = l + nv * (k - 1) + nvnu * (i - 1)
1805 tmp = tmp + work(jj) * bt(k,j)
1806 end do
1807 work2(ii) = tmp
1808 end do
1809 end do
1810 end do
1811
1812 do j = 1, nv
1813 do i = 1, nvnv
1814 jj = i + nvnv * (j - 1) + iv
1815 v(jj) = work2(i + nvnv * (1 - 1)) * ct(1, j) &
1816 + work2(i + nvnv * (2 - 1)) * ct(2, j) &
1817 + work2(i + nvnv * (3 - 1)) * ct(3, j) &
1818 + work2(i + nvnv * (4 - 1)) * ct(4, j)
1819
1820 end do
1821 end do
1822 end do
1823 !$omp end parallel do
1824 end subroutine tnsr1_3d_nu4nv2_cpu
1825
1826end module tensor_cpu
mxm_wrapper
Wrapper for all matrix-matrix product implementations.
Definition mxm_wrapper.F90:2
mxm_wrapper::mxm
subroutine, public mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product for contiguously packed matrices A,B, and C.
Definition mxm_wrapper.F90:31
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
tensor_cpu
Definition tensor_cpu.f90:1
tensor_cpu::tnsr3d_el_1_12_cpu
subroutine tnsr3d_el_1_12_cpu(v, u, a, bt, ct)
Single-element evaluation at one point, specialised for nu = 12.
Definition tensor_cpu.f90:396
tensor_cpu::tnsr3d_el_n5_cpu
subroutine tnsr3d_el_n5_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 5.
Definition tensor_cpu.f90:1097
tensor_cpu::tnsr3d_el_1_8_cpu
subroutine tnsr3d_el_1_8_cpu(v, u, a, bt, ct)
Single-element evaluation at one point, specialised for nu = 8.
Definition tensor_cpu.f90:297
tensor_cpu::tnsr3d_el_1_6_cpu
subroutine tnsr3d_el_1_6_cpu(v, u, a, bt, ct)
Single-element evaluation at one point, specialised for nu = 6.
Definition tensor_cpu.f90:257
tensor_cpu::tnsr1_3d_nvnu_cpu
subroutine tnsr1_3d_nvnu_cpu(v, nv, nu, a, bt, ct, nelv)
In-place tensor product for arbitrary nu and nv.
Definition tensor_cpu.f90:1693
tensor_cpu::tnsr3d_nu12_cpu
subroutine tnsr3d_nu12_cpu(v, nv, u, a, bt, ct, nelv)
Tensor-product evaluation specialised for nu = 12, generic in nv.
Definition tensor_cpu.f90:1587
tensor_cpu::tnsr3d_el_n13_cpu
subroutine tnsr3d_el_n13_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 13.
Definition tensor_cpu.f90:589
tensor_cpu::tnsr3d_cpu
subroutine, public tnsr3d_cpu(v, nv, u, nu, a, bt, ct, nelv)
Three-dimensional tensor product over a list of elements.
Definition tensor_cpu.f90:1292
tensor_cpu::tnsr3d_el_n3_cpu
subroutine tnsr3d_el_n3_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 3.
Definition tensor_cpu.f90:1194
tensor_cpu::tnsr3d_nu2nv4_cpu
subroutine tnsr3d_nu2nv4_cpu(v, u, a, bt, ct, nelv)
Batched tensor product specialised for nu = 2, nv = 4.
Definition tensor_cpu.f90:1403
tensor_cpu::tnsr3d_nvnu_cpu
subroutine tnsr3d_nvnu_cpu(v, nv, u, nu, a, bt, ct, nelv)
Generic tensor-product evaluation for arbitrary nu, nv.
Definition tensor_cpu.f90:1324
tensor_cpu::tnsr3d_el_n7_cpu
subroutine tnsr3d_el_n7_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 7.
Definition tensor_cpu.f90:988
tensor_cpu::tnsr3d_el_n11_cpu
subroutine tnsr3d_el_n11_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 11.
Definition tensor_cpu.f90:734
tensor_cpu::tnsr3d_el_1_nu_cpu
subroutine tnsr3d_el_1_nu_cpu(v, u, nu, a, bt, ct)
Single-element evaluation at one point, for arbitrary nu.
Definition tensor_cpu.f90:183
tensor_cpu::tnsr3d_nu8_cpu
subroutine tnsr3d_nu8_cpu(v, nv, u, a, bt, ct, nelv)
Tensor-product evaluation specialised for nu = 8, generic in nv.
Definition tensor_cpu.f90:1516
tensor_cpu::tnsr3d_el_n9_cpu
subroutine tnsr3d_el_n9_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 9.
Definition tensor_cpu.f90:867
tensor_cpu::tnsr3d_el_n6_cpu
subroutine tnsr3d_el_n6_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 6.
Definition tensor_cpu.f90:1044
tensor_cpu::tnsr3d_el_n2_cpu
subroutine tnsr3d_el_n2_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 2.
Definition tensor_cpu.f90:1238
tensor_cpu::tnsr3d_el_n12_cpu
subroutine tnsr3d_el_n12_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 12.
Definition tensor_cpu.f90:663
tensor_cpu::tnsr3d_el_1_4_cpu
subroutine tnsr3d_el_1_4_cpu(v, u, a, bt, ct)
Single-element evaluation at one point, specialised for nu = 4.
Definition tensor_cpu.f90:223
tensor_cpu::tnsr3d_el_n4_cpu
subroutine tnsr3d_el_n4_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 4.
Definition tensor_cpu.f90:1147
tensor_cpu::tnsr3d_el_1_10_cpu
subroutine tnsr3d_el_1_10_cpu(v, u, a, bt, ct)
Single-element evaluation at one point, specialised for nu = 10.
Definition tensor_cpu.f90:344
tensor_cpu::tnsr1_3d_nu4nv2_cpu
subroutine tnsr1_3d_nu4nv2_cpu(v, a, bt, ct, nelv)
In-place tensor product specialised for nu = 4, nv = 2.
Definition tensor_cpu.f90:1768
tensor_cpu::tnsr2d_el_cpu
subroutine, public tnsr2d_el_cpu(v, nv, u, nu, a, bt)
Two-dimensional tensor product on a single element.
Definition tensor_cpu.f90:23
tensor_cpu::tnsr1_3d_cpu
subroutine, public tnsr1_3d_cpu(v, nv, nu, a, bt, ct, nelv)
Three-dimensional tensor product applied in place.
Definition tensor_cpu.f90:1676
tensor_cpu::tnsr3d_nu4_cpu
subroutine tnsr3d_nu4_cpu(v, nv, u, a, bt, ct, nelv)
Batched tensor product specialised for nu = 4, generic in nv.
Definition tensor_cpu.f90:1455
tensor_cpu::tnsr3d_el_cpu
subroutine, public tnsr3d_el_cpu(v, nv, u, nu, a, bt, ct)
Three-dimensional tensor product on a single element.
Definition tensor_cpu.f90:50
tensor_cpu::tnsr3d_el_n_cpu
subroutine tnsr3d_el_n_cpu(v, u, a, bt, ct, n)
Single-element tensor product with nv = nu, for arbitrary order.
Definition tensor_cpu.f90:461
tensor_cpu::tnsr3d_el_nvnu_cpu
subroutine tnsr3d_el_nvnu_cpu(v, nv, u, nu, a, bt, ct)
Single-element tensor product for arbitrary nu and nv.
Definition tensor_cpu.f90:121
tensor_cpu::tnsr3d_el_n14_cpu
subroutine tnsr3d_el_n14_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 14.
Definition tensor_cpu.f90:512
tensor_cpu::tnsr3d_el_n10_cpu
subroutine tnsr3d_el_n10_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 10.
Definition tensor_cpu.f90:802
tensor_cpu::tnsr3d_el_n8_cpu
subroutine tnsr3d_el_n8_cpu(v, u, a, bt, ct)
Single-element tensor product specialised for nv = nu = 8.
Definition tensor_cpu.f90:929