FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_hevi_common_linalgebra.F90
Go to the documentation of this file.
1
2
3!-------------------------------------------------------------------------------
4! Warning: This file was generated from fluid_dyn_solver/scale_atm_dyn_dgm_hevi_common_linalgebra.F90.erb.
5! Do not edit this file.
6!-------------------------------------------------------------------------------
7!-------------------------------------------------------------------------------
8!> module FElib / Fluid dyn solver / Atmosphere / HEVI / Common
9!!
10!! @par Description
11!! Computational kernels associated with linear algebra for vertical implicit solvers
12!!
13!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
14!<
15!-------------------------------------------------------------------------------
16#include "scaleFElib.h"
17module scale_atm_dyn_dgm_hevi_common_linalgebra
18 !-----------------------------------------------------------------------------
19 !
20 !++ Used modules
21 !
22 use scale_precision
23 use scale_io
24 use scale_prc
25 use scale_prof
26 !-----------------------------------------------------------------------------
27 implicit none
28 private
29 !-----------------------------------------------------------------------------
30 !
31 !++ Public procedures
32 !
33 public :: atm_dyn_dgm_hevi_common_linalgebra_get_param
34 public :: atm_dyn_dgm_hevi_common_linalgebra_solve_uv
35 public :: atm_dyn_dgm_hevi_common_linalgebra_solve_var3
36
37 !-----------------------------------------------------------------------------
38 !
39 !++ Public parameters & variables
40 !
41
42 !-----------------------------------------------------------------------------
43 !
44 !++ Private procedures & variables
45 !
46 !-------------------
47contains
48!OCL SERIAL
49 subroutine atm_dyn_dgm_hevi_common_linalgebra_get_param( im, jm, Nnode_h1D )
50 implicit none
51 integer, intent(out) :: im
52 integer, intent(out) :: jm
53 integer, intent(in) :: nnode_h1d
54 !-----------------------------------------------------------
55 select case(nnode_h1d)
56 case( 2 )
57 im = 4
58 case( 3 )
59 im = 9
60 case( 4 )
61 im = 16
62 case( 5 )
63 im = 5
64 case( 6 )
65 im = 6
66 case( 7 )
67 im = 7
68 case( 8 )
69 im = 8
70 case( 9 )
71 im = 9
72 case( 10 )
73 im = 10
74 case( 11 )
75 im = 22
76 case( 12 )
77 im = 24
78 case( 13 )
79 im = 26
80 case( 14 )
81 im = 28
82 case( 15 )
83 im = 15
84 case( 16 )
85 im = 16
86 end select
87 jm = nnode_h1d**2/im
88
89 return
90 end subroutine atm_dyn_dgm_hevi_common_linalgebra_get_param
91
92!OCL SERIAL
93 subroutine atm_dyn_dgm_hevi_common_linalgebra_solve_uv( D, b, G, Nnode_v, im, jm, Ne2D, is_top )
94 implicit none
95 integer, intent(in) :: nnode_v
96 integer, intent(in) :: im, jm
97 integer, intent(in) :: ne2d
98 real(rp), intent(in) :: d(im,nnode_v,nnode_v,jm,ne2d)
99 real(rp), intent(inout) :: b(im,nnode_v,2,jm,ne2d)
100 real(rp), intent(inout) :: g(im,nnode_v,jm,ne2d)
101 logical, intent(in) :: is_top
102 !------------------------------------------
103 select case(nnode_v)
104 case( 2 )
105 if ( im == 4 ) then
106 call solve_nnode2_uv( d, b, g, nnode_v, jm, ne2d, is_top )
107 else
108 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
109 end if
110 case( 3 )
111 if ( im == 9 ) then
112 call solve_nnode3_uv( d, b, g, nnode_v, jm, ne2d, is_top )
113 else
114 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
115 end if
116 case( 4 )
117 if ( im == 16 ) then
118 call solve_nnode4_uv( d, b, g, nnode_v, jm, ne2d, is_top )
119 else
120 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
121 end if
122 case( 5 )
123 if ( im == 5 ) then
124 call solve_nnode5_uv( d, b, g, nnode_v, jm, ne2d, is_top )
125 else
126 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
127 end if
128 case( 6 )
129 if ( im == 6 ) then
130 call solve_nnode6_uv( d, b, g, nnode_v, jm, ne2d, is_top )
131 else
132 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
133 end if
134 case( 7 )
135 if ( im == 7 ) then
136 call solve_nnode7_uv( d, b, g, nnode_v, jm, ne2d, is_top )
137 else
138 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
139 end if
140 case( 8 )
141 if ( im == 8 ) then
142 call solve_nnode8_uv( d, b, g, nnode_v, jm, ne2d, is_top )
143 else
144 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
145 end if
146 case( 9 )
147 if ( im == 9 ) then
148 call solve_nnode9_uv( d, b, g, nnode_v, jm, ne2d, is_top )
149 else
150 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
151 end if
152 case( 10 )
153 if ( im == 10 ) then
154 call solve_nnode10_uv( d, b, g, nnode_v, jm, ne2d, is_top )
155 else
156 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
157 end if
158 case( 11 )
159 if ( im == 22 ) then
160 call solve_nnode11_uv( d, b, g, nnode_v, jm, ne2d, is_top )
161 else
162 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
163 end if
164 case( 12 )
165 if ( im == 24 ) then
166 call solve_nnode12_uv( d, b, g, nnode_v, jm, ne2d, is_top )
167 else
168 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
169 end if
170 case( 13 )
171 if ( im == 26 ) then
172 call solve_nnode13_uv( d, b, g, nnode_v, jm, ne2d, is_top )
173 else
174 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
175 end if
176 case( 14 )
177 if ( im == 28 ) then
178 call solve_nnode14_uv( d, b, g, nnode_v, jm, ne2d, is_top )
179 else
180 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
181 end if
182 case( 15 )
183 if ( im == 15 ) then
184 call solve_nnode15_uv( d, b, g, nnode_v, jm, ne2d, is_top )
185 else
186 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
187 end if
188 case( 16 )
189 if ( im == 16 ) then
190 call solve_nnode16_uv( d, b, g, nnode_v, jm, ne2d, is_top )
191 else
192 call solve_uv_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
193 end if
194 case default
195 log_info('atm_dyn_dgm_hevi_common_linalgebra_solve_uv',*) "im > 16 is not supported. Check!"
196 call prc_abort
197 end select
198 return
199 end subroutine atm_dyn_dgm_hevi_common_linalgebra_solve_uv
200
201!OCL SERIAL
202 subroutine atm_dyn_dgm_hevi_common_linalgebra_solve_var3( D, b, G, Nnode_v, im, jm, Ne2D, is_top)
203 implicit none
204 integer, intent(in) :: nnode_v
205 integer, intent(in) :: im, jm
206 integer, intent(in) :: ne2d
207 real(rp), intent(in) :: d(im,3*nnode_v,3*nnode_v,jm,ne2d)
208 real(rp), intent(inout) :: b(im,3*nnode_v,jm,ne2d)
209 real(rp), intent(inout) :: g(im,3*nnode_v,3,jm,ne2d)
210 logical, intent(in) :: is_top
211 !------------------------------------------
212 select case(nnode_v)
213 case( 2 )
214 if ( im == 4 ) then
215 call solve_nnode2_var3( d, b, g, nnode_v, jm, ne2d, is_top )
216 else
217 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
218 end if
219 case( 3 )
220 if ( im == 9 ) then
221 call solve_nnode3_var3( d, b, g, nnode_v, jm, ne2d, is_top )
222 else
223 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
224 end if
225 case( 4 )
226 if ( im == 16 ) then
227 call solve_nnode4_var3( d, b, g, nnode_v, jm, ne2d, is_top )
228 else
229 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
230 end if
231 case( 5 )
232 if ( im == 5 ) then
233 call solve_nnode5_var3( d, b, g, nnode_v, jm, ne2d, is_top )
234 else
235 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
236 end if
237 case( 6 )
238 if ( im == 6 ) then
239 call solve_nnode6_var3( d, b, g, nnode_v, jm, ne2d, is_top )
240 else
241 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
242 end if
243 case( 7 )
244 if ( im == 7 ) then
245 call solve_nnode7_var3( d, b, g, nnode_v, jm, ne2d, is_top )
246 else
247 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
248 end if
249 case( 8 )
250 if ( im == 8 ) then
251 call solve_nnode8_var3( d, b, g, nnode_v, jm, ne2d, is_top )
252 else
253 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
254 end if
255 case( 9 )
256 if ( im == 9 ) then
257 call solve_nnode9_var3( d, b, g, nnode_v, jm, ne2d, is_top )
258 else
259 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
260 end if
261 case( 10 )
262 if ( im == 10 ) then
263 call solve_nnode10_var3( d, b, g, nnode_v, jm, ne2d, is_top )
264 else
265 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
266 end if
267 case( 11 )
268 if ( im == 22 ) then
269 call solve_nnode11_var3( d, b, g, nnode_v, jm, ne2d, is_top )
270 else
271 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
272 end if
273 case( 12 )
274 if ( im == 24 ) then
275 call solve_nnode12_var3( d, b, g, nnode_v, jm, ne2d, is_top )
276 else
277 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
278 end if
279 case( 13 )
280 if ( im == 26 ) then
281 call solve_nnode13_var3( d, b, g, nnode_v, jm, ne2d, is_top )
282 else
283 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
284 end if
285 case( 14 )
286 if ( im == 28 ) then
287 call solve_nnode14_var3( d, b, g, nnode_v, jm, ne2d, is_top )
288 else
289 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
290 end if
291 case( 15 )
292 if ( im == 15 ) then
293 call solve_nnode15_var3( d, b, g, nnode_v, jm, ne2d, is_top )
294 else
295 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
296 end if
297 case( 16 )
298 if ( im == 16 ) then
299 call solve_nnode16_var3( d, b, g, nnode_v, jm, ne2d, is_top )
300 else
301 call solve_var3_im( d, b, g, nnode_v, im, jm, ne2d, is_top )
302 end if
303 case default
304 log_info('atm_dyn_dgm_hevi_common_linalgebra_solve',*) "im > 16 is not supported. Check!"
305 call prc_abort
306 end select
307 return
308 end subroutine atm_dyn_dgm_hevi_common_linalgebra_solve_var3
309
310!------
311!OCL SERIAL
312 subroutine solve_nnode2_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
313 implicit none
314 integer, intent(in) :: nnode_v
315 integer, intent(in) :: jm
316 integer, intent(in) :: ne2d
317 real(rp), intent(in) :: d(4,2,2,jm,ne2d)
318 real(rp), intent(inout) :: b(4,2,2,jm,ne2d)
319 real(rp), intent(inout) :: g(4,2,jm,ne2d)
320 logical, intent(in) :: is_top
321
322 integer :: ke_xy
323 integer :: v2
324 integer :: k, i, j, v, r
325 integer :: ip
326 integer :: ipiv(4,2)
327 real(rp) :: tmpv(4)
328 real(rp) :: tmpv2(4)
329 real(rp) :: a(4,2,2)
330 real(rp) :: bb(4,3,2)
331 real(rp) :: tmpb(4,3)
332 real(rp) :: tmp
333 !------------------------------
334
335 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
336 do ke_xy=1, ne2d
337 do v2=1, jm
338 a(:,:,:) = d(:,:,:,v2,ke_xy)
339 do k=1, 2
340 tmpv(:) = abs(a(:,k,k))
341 ipiv(:,k) = k
342 do i=k+1, 2
343 do v=1, 4
344 tmp = abs(a(v,i,k))
345 if ( tmp > tmpv(v) ) then
346 ipiv(v,k) = i
347 tmpv(v) = tmp
348 end if
349 end do
350 end do
351 do v=1, 4
352 ip = ipiv(v,k)
353 if ( ip /= k ) then
354 do j=1, 2
355 tmp = a(v,k,j)
356 a(v,k,j) = a(v,ip,j)
357 a(v,ip,j) = tmp
358 end do
359 end if
360 end do
361 !--
362 tmpv(:) = 1.0_rp / a(:,k,k)
363 a(:,k,k) = tmpv(:)
364 do i=k+1, 2
365 a(:,i,k) = a(:,i,k) * tmpv(:)
366 end do
367 do j=k+1, 2
368 do i=k+1, 2
369 do v=1, 4
370 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
371 end do
372 end do
373 end do
374 end do
375 !---------------------------
376 if ( is_top ) then
377 do i=1, 2-1
378 do v=1, 4
379 ip = ipiv(v,i)
380 if (ip /= i) then
381 tmp = b(v,i,1,v2,ke_xy)
382 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
383 b(v,ip,1,v2,ke_xy) = tmp
384 tmp = b(v,i,2,v2,ke_xy)
385 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
386 b(v,ip,2,v2,ke_xy) = tmp
387 end if
388 end do
389 end do
390 do i=1, 2
391 tmpv(:) = b(:,i,1,v2,ke_xy)
392 tmpv2(:) = b(:,i,2,v2,ke_xy)
393 do j=1, i-1
394 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
395 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
396 end do
397 b(:,i,1,v2,ke_xy) = tmpv(:)
398 b(:,i,2,v2,ke_xy) = tmpv2(:)
399 end do
400 do i=2, 1, -1
401 tmpv(:) = b(:,i,1,v2,ke_xy)
402 tmpv2(:) = b(:,i,2,v2,ke_xy)
403 do j=i+1, 2
404 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
405 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
406 end do
407 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
408 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
409 end do
410 else
411 do i=1, 2-1
412 do v=1, 4
413 ip = ipiv(v,i)
414 if (ip /= i) then
415 tmp = b(v,i,1,v2,ke_xy)
416 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
417 b(v,ip,1,v2,ke_xy) = tmp
418 tmp = b(v,i,2,v2,ke_xy)
419 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
420 b(v,ip,2,v2,ke_xy) = tmp
421
422 tmp = g(v,i,v2,ke_xy)
423 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
424 g(v,ip,v2,ke_xy) = tmp
425 end if
426 end do
427 end do
428 do i=1, 2
429 do v=1,4
430 bb(v,1,i) = b(v,i,1,v2,ke_xy)
431 bb(v,2,i) = b(v,i,2,v2,ke_xy)
432 bb(v,3,i) = g(v,i,v2,ke_xy)
433 end do
434 end do
435 do i=1, 2
436 tmpb(:,:) = bb(:,:,i)
437 do j=1, i-1
438 do r=1, 3
439 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
440 end do
441 end do
442 bb(:,:,i) = tmpb(:,:)
443 end do
444 do i=2, 1, -1
445 tmpb(:,:) = bb(:,:,i)
446 do j=i+1, 2
447 do r=1, 3
448 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
449 end do
450 end do
451 do r=1, 3
452 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
453 end do
454 end do
455 do i=1, 2
456 b(:,i,1,v2,ke_xy) = bb(:,1,i)
457 b(:,i,2,v2,ke_xy) = bb(:,2,i)
458 g(:,i,v2,ke_xy) = bb(:,3,i)
459 end do
460 end if
461 end do
462 end do
463 return
464 end subroutine solve_nnode2_uv
465!OCL SERIAL
466 subroutine solve_nnode2_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
467 implicit none
468 integer, intent(in) :: nnode_v
469 integer, intent(in) :: jm
470 integer, intent(in) :: ne2d
471 real(rp), intent(in) :: d(4,6,6,jm,ne2d)
472 real(rp), intent(inout) :: b(4,6,jm,ne2d)
473 real(rp), intent(inout) :: g(4,6,3,jm,ne2d)
474 logical, intent(in) :: is_top
475
476 integer :: ke_xy
477 integer :: v2
478 integer :: k, i, j, v, r
479 integer :: ip
480 integer :: ipiv(4,6)
481 real(rp) :: tmpv(4)
482 real(rp) :: a(4,6,6)
483 real(rp) :: bb(4,4,6)
484 real(rp) :: tmpb(4,4)
485 real(rp) :: tmp
486 !------------------------------
487
488 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
489 do ke_xy=1, ne2d
490 do v2=1, jm
491 a(:,:,:) = d(:,:,:,v2,ke_xy)
492 do k=1, 6
493 tmpv(:) = abs(a(:,k,k))
494 ipiv(:,k) = k
495 do i=k+1, 6
496 do v=1, 4
497 tmp = abs(a(v,i,k))
498 if ( tmp > tmpv(v) ) then
499 ipiv(v,k) = i
500 tmpv(v) = tmp
501 end if
502 end do
503 end do
504 do v=1, 4
505 ip = ipiv(v,k)
506 if ( ip /= k ) then
507 do j=1, 6
508 tmp = a(v,k,j)
509 a(v,k,j) = a(v,ip,j)
510 a(v,ip,j) = tmp
511 end do
512 end if
513 end do
514 !--
515 tmpv(:) = 1.0_rp / a(:,k,k)
516 a(:,k,k) = tmpv(:)
517 do i=k+1, 6
518 do v=1, 4
519 a(v,i,k) = a(v,i,k) * tmpv(v)
520 end do
521 end do
522 do j=k+1, 6
523 do i=k+1, 6
524 do v=1, 4
525 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
526 end do
527 end do
528 end do
529 end do
530 !---------------------------
531 if ( is_top ) then
532 do i=1, 6-1
533 do v=1, 4
534 ip = ipiv(v,i)
535 if (ip /= i) then
536 tmp = b(v,i,v2,ke_xy)
537 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
538 b(v,ip,v2,ke_xy) = tmp
539 end if
540 end do
541 end do
542 do i=1, 6
543 tmpv(:) = b(:,i,v2,ke_xy)
544 do j=1, i-1
545 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
546 end do
547 b(:,i,v2,ke_xy) = tmpv(:)
548 end do
549 do i=6, 1, -1
550 tmpv(:) = b(:,i,v2,ke_xy)
551 do j=i+1, 6
552 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
553 end do
554 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
555 end do
556 else
557 do i=1, 6-1
558 do v=1, 4
559 ip = ipiv(v,i)
560 if (ip /= i) then
561 tmp = b(v,i,v2,ke_xy)
562 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
563 b(v,ip,v2,ke_xy) = tmp
564
565 tmp = g(v,i,1,v2,ke_xy)
566 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
567 g(v,ip,1,v2,ke_xy) = tmp
568 tmp = g(v,i,2,v2,ke_xy)
569 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
570 g(v,ip,2,v2,ke_xy) = tmp
571 tmp = g(v,i,3,v2,ke_xy)
572 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
573 g(v,ip,3,v2,ke_xy) = tmp
574 end if
575 end do
576 end do
577 do i=1, 6
578 do v=1,4
579 bb(v,1,i) = b(v,i,v2,ke_xy)
580 bb(v,2,i) = g(v,i,1,v2,ke_xy)
581 bb(v,3,i) = g(v,i,2,v2,ke_xy)
582 bb(v,4,i) = g(v,i,3,v2,ke_xy)
583 end do
584 end do
585 do i=1, 6
586 tmpb(:,:) = bb(:,:,i)
587 do j=1, i-1
588 do r=1, 4
589 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
590 end do
591 end do
592 bb(:,:,i) = tmpb(:,:)
593 end do
594 do i=6, 1, -1
595 tmpb(:,:) = bb(:,:,i)
596 do j=i+1, 6
597 do r=1, 4
598 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
599 end do
600 end do
601 do r=1, 4
602 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
603 end do
604 end do
605 do i=1, 6
606 b(:,i,v2,ke_xy) = bb(:,1,i)
607 g(:,i,1,v2,ke_xy) = bb(:,2,i)
608 g(:,i,2,v2,ke_xy) = bb(:,3,i)
609 g(:,i,3,v2,ke_xy) = bb(:,4,i)
610 end do
611 end if
612 end do
613 end do
614 return
615 end subroutine solve_nnode2_var3
616!OCL SERIAL
617 subroutine solve_nnode3_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
618 implicit none
619 integer, intent(in) :: nnode_v
620 integer, intent(in) :: jm
621 integer, intent(in) :: ne2d
622 real(rp), intent(in) :: d(9,3,3,jm,ne2d)
623 real(rp), intent(inout) :: b(9,3,2,jm,ne2d)
624 real(rp), intent(inout) :: g(9,3,jm,ne2d)
625 logical, intent(in) :: is_top
626
627 integer :: ke_xy
628 integer :: v2
629 integer :: k, i, j, v, r
630 integer :: ip
631 integer :: ipiv(9,3)
632 real(rp) :: tmpv(9)
633 real(rp) :: tmpv2(9)
634 real(rp) :: a(9,3,3)
635 real(rp) :: bb(9,3,3)
636 real(rp) :: tmpb(9,3)
637 real(rp) :: tmp
638 !------------------------------
639
640 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
641 do ke_xy=1, ne2d
642 do v2=1, jm
643 a(:,:,:) = d(:,:,:,v2,ke_xy)
644 do k=1, 3
645 tmpv(:) = abs(a(:,k,k))
646 ipiv(:,k) = k
647 do i=k+1, 3
648 do v=1, 9
649 tmp = abs(a(v,i,k))
650 if ( tmp > tmpv(v) ) then
651 ipiv(v,k) = i
652 tmpv(v) = tmp
653 end if
654 end do
655 end do
656 do v=1, 9
657 ip = ipiv(v,k)
658 if ( ip /= k ) then
659 do j=1, 3
660 tmp = a(v,k,j)
661 a(v,k,j) = a(v,ip,j)
662 a(v,ip,j) = tmp
663 end do
664 end if
665 end do
666 !--
667 tmpv(:) = 1.0_rp / a(:,k,k)
668 a(:,k,k) = tmpv(:)
669 do i=k+1, 3
670 a(:,i,k) = a(:,i,k) * tmpv(:)
671 end do
672 do j=k+1, 3
673 do i=k+1, 3
674 do v=1, 9
675 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
676 end do
677 end do
678 end do
679 end do
680 !---------------------------
681 if ( is_top ) then
682 do i=1, 3-1
683 do v=1, 9
684 ip = ipiv(v,i)
685 if (ip /= i) then
686 tmp = b(v,i,1,v2,ke_xy)
687 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
688 b(v,ip,1,v2,ke_xy) = tmp
689 tmp = b(v,i,2,v2,ke_xy)
690 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
691 b(v,ip,2,v2,ke_xy) = tmp
692 end if
693 end do
694 end do
695 do i=1, 3
696 tmpv(:) = b(:,i,1,v2,ke_xy)
697 tmpv2(:) = b(:,i,2,v2,ke_xy)
698 do j=1, i-1
699 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
700 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
701 end do
702 b(:,i,1,v2,ke_xy) = tmpv(:)
703 b(:,i,2,v2,ke_xy) = tmpv2(:)
704 end do
705 do i=3, 1, -1
706 tmpv(:) = b(:,i,1,v2,ke_xy)
707 tmpv2(:) = b(:,i,2,v2,ke_xy)
708 do j=i+1, 3
709 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
710 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
711 end do
712 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
713 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
714 end do
715 else
716 do i=1, 3-1
717 do v=1, 9
718 ip = ipiv(v,i)
719 if (ip /= i) then
720 tmp = b(v,i,1,v2,ke_xy)
721 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
722 b(v,ip,1,v2,ke_xy) = tmp
723 tmp = b(v,i,2,v2,ke_xy)
724 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
725 b(v,ip,2,v2,ke_xy) = tmp
726
727 tmp = g(v,i,v2,ke_xy)
728 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
729 g(v,ip,v2,ke_xy) = tmp
730 end if
731 end do
732 end do
733 do i=1, 3
734 do v=1,9
735 bb(v,1,i) = b(v,i,1,v2,ke_xy)
736 bb(v,2,i) = b(v,i,2,v2,ke_xy)
737 bb(v,3,i) = g(v,i,v2,ke_xy)
738 end do
739 end do
740 do i=1, 3
741 tmpb(:,:) = bb(:,:,i)
742 do j=1, i-1
743 do r=1, 3
744 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
745 end do
746 end do
747 bb(:,:,i) = tmpb(:,:)
748 end do
749 do i=3, 1, -1
750 tmpb(:,:) = bb(:,:,i)
751 do j=i+1, 3
752 do r=1, 3
753 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
754 end do
755 end do
756 do r=1, 3
757 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
758 end do
759 end do
760 do i=1, 3
761 b(:,i,1,v2,ke_xy) = bb(:,1,i)
762 b(:,i,2,v2,ke_xy) = bb(:,2,i)
763 g(:,i,v2,ke_xy) = bb(:,3,i)
764 end do
765 end if
766 end do
767 end do
768 return
769 end subroutine solve_nnode3_uv
770!OCL SERIAL
771 subroutine solve_nnode3_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
772 implicit none
773 integer, intent(in) :: nnode_v
774 integer, intent(in) :: jm
775 integer, intent(in) :: ne2d
776 real(rp), intent(in) :: d(9,9,9,jm,ne2d)
777 real(rp), intent(inout) :: b(9,9,jm,ne2d)
778 real(rp), intent(inout) :: g(9,9,3,jm,ne2d)
779 logical, intent(in) :: is_top
780
781 integer :: ke_xy
782 integer :: v2
783 integer :: k, i, j, v, r
784 integer :: ip
785 integer :: ipiv(9,9)
786 real(rp) :: tmpv(9)
787 real(rp) :: a(9,9,9)
788 real(rp) :: bb(9,4,9)
789 real(rp) :: tmpb(9,4)
790 real(rp) :: tmp
791 !------------------------------
792
793 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
794 do ke_xy=1, ne2d
795 do v2=1, jm
796 a(:,:,:) = d(:,:,:,v2,ke_xy)
797 do k=1, 9
798 tmpv(:) = abs(a(:,k,k))
799 ipiv(:,k) = k
800 do i=k+1, 9
801 do v=1, 9
802 tmp = abs(a(v,i,k))
803 if ( tmp > tmpv(v) ) then
804 ipiv(v,k) = i
805 tmpv(v) = tmp
806 end if
807 end do
808 end do
809 do v=1, 9
810 ip = ipiv(v,k)
811 if ( ip /= k ) then
812 do j=1, 9
813 tmp = a(v,k,j)
814 a(v,k,j) = a(v,ip,j)
815 a(v,ip,j) = tmp
816 end do
817 end if
818 end do
819 !--
820 tmpv(:) = 1.0_rp / a(:,k,k)
821 a(:,k,k) = tmpv(:)
822 do i=k+1, 9
823 do v=1, 9
824 a(v,i,k) = a(v,i,k) * tmpv(v)
825 end do
826 end do
827 do j=k+1, 9
828 do i=k+1, 9
829 do v=1, 9
830 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
831 end do
832 end do
833 end do
834 end do
835 !---------------------------
836 if ( is_top ) then
837 do i=1, 9-1
838 do v=1, 9
839 ip = ipiv(v,i)
840 if (ip /= i) then
841 tmp = b(v,i,v2,ke_xy)
842 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
843 b(v,ip,v2,ke_xy) = tmp
844 end if
845 end do
846 end do
847 do i=1, 9
848 tmpv(:) = b(:,i,v2,ke_xy)
849 do j=1, i-1
850 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
851 end do
852 b(:,i,v2,ke_xy) = tmpv(:)
853 end do
854 do i=9, 1, -1
855 tmpv(:) = b(:,i,v2,ke_xy)
856 do j=i+1, 9
857 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
858 end do
859 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
860 end do
861 else
862 do i=1, 9-1
863 do v=1, 9
864 ip = ipiv(v,i)
865 if (ip /= i) then
866 tmp = b(v,i,v2,ke_xy)
867 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
868 b(v,ip,v2,ke_xy) = tmp
869
870 tmp = g(v,i,1,v2,ke_xy)
871 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
872 g(v,ip,1,v2,ke_xy) = tmp
873 tmp = g(v,i,2,v2,ke_xy)
874 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
875 g(v,ip,2,v2,ke_xy) = tmp
876 tmp = g(v,i,3,v2,ke_xy)
877 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
878 g(v,ip,3,v2,ke_xy) = tmp
879 end if
880 end do
881 end do
882 do i=1, 9
883 do v=1,9
884 bb(v,1,i) = b(v,i,v2,ke_xy)
885 bb(v,2,i) = g(v,i,1,v2,ke_xy)
886 bb(v,3,i) = g(v,i,2,v2,ke_xy)
887 bb(v,4,i) = g(v,i,3,v2,ke_xy)
888 end do
889 end do
890 do i=1, 9
891 tmpb(:,:) = bb(:,:,i)
892 do j=1, i-1
893 do r=1, 4
894 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
895 end do
896 end do
897 bb(:,:,i) = tmpb(:,:)
898 end do
899 do i=9, 1, -1
900 tmpb(:,:) = bb(:,:,i)
901 do j=i+1, 9
902 do r=1, 4
903 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
904 end do
905 end do
906 do r=1, 4
907 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
908 end do
909 end do
910 do i=1, 9
911 b(:,i,v2,ke_xy) = bb(:,1,i)
912 g(:,i,1,v2,ke_xy) = bb(:,2,i)
913 g(:,i,2,v2,ke_xy) = bb(:,3,i)
914 g(:,i,3,v2,ke_xy) = bb(:,4,i)
915 end do
916 end if
917 end do
918 end do
919 return
920 end subroutine solve_nnode3_var3
921!OCL SERIAL
922 subroutine solve_nnode4_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
923 implicit none
924 integer, intent(in) :: nnode_v
925 integer, intent(in) :: jm
926 integer, intent(in) :: ne2d
927 real(rp), intent(in) :: d(16,4,4,jm,ne2d)
928 real(rp), intent(inout) :: b(16,4,2,jm,ne2d)
929 real(rp), intent(inout) :: g(16,4,jm,ne2d)
930 logical, intent(in) :: is_top
931
932 integer :: ke_xy
933 integer :: v2
934 integer :: k, i, j, v, r
935 integer :: ip
936 integer :: ipiv(16,4)
937 real(rp) :: tmpv(16)
938 real(rp) :: tmpv2(16)
939 real(rp) :: a(16,4,4)
940 real(rp) :: bb(16,3,4)
941 real(rp) :: tmpb(16,3)
942 real(rp) :: tmp
943 !------------------------------
944
945 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
946 do ke_xy=1, ne2d
947 do v2=1, jm
948 a(:,:,:) = d(:,:,:,v2,ke_xy)
949 do k=1, 4
950 tmpv(:) = abs(a(:,k,k))
951 ipiv(:,k) = k
952 do i=k+1, 4
953 do v=1, 16
954 tmp = abs(a(v,i,k))
955 if ( tmp > tmpv(v) ) then
956 ipiv(v,k) = i
957 tmpv(v) = tmp
958 end if
959 end do
960 end do
961 do v=1, 16
962 ip = ipiv(v,k)
963 if ( ip /= k ) then
964 do j=1, 4
965 tmp = a(v,k,j)
966 a(v,k,j) = a(v,ip,j)
967 a(v,ip,j) = tmp
968 end do
969 end if
970 end do
971 !--
972 tmpv(:) = 1.0_rp / a(:,k,k)
973 a(:,k,k) = tmpv(:)
974 do i=k+1, 4
975 a(:,i,k) = a(:,i,k) * tmpv(:)
976 end do
977 do j=k+1, 4
978 do i=k+1, 4
979 do v=1, 16
980 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
981 end do
982 end do
983 end do
984 end do
985 !---------------------------
986 if ( is_top ) then
987 do i=1, 4-1
988 do v=1, 16
989 ip = ipiv(v,i)
990 if (ip /= i) then
991 tmp = b(v,i,1,v2,ke_xy)
992 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
993 b(v,ip,1,v2,ke_xy) = tmp
994 tmp = b(v,i,2,v2,ke_xy)
995 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
996 b(v,ip,2,v2,ke_xy) = tmp
997 end if
998 end do
999 end do
1000 do i=1, 4
1001 tmpv(:) = b(:,i,1,v2,ke_xy)
1002 tmpv2(:) = b(:,i,2,v2,ke_xy)
1003 do j=1, i-1
1004 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
1005 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
1006 end do
1007 b(:,i,1,v2,ke_xy) = tmpv(:)
1008 b(:,i,2,v2,ke_xy) = tmpv2(:)
1009 end do
1010 do i=4, 1, -1
1011 tmpv(:) = b(:,i,1,v2,ke_xy)
1012 tmpv2(:) = b(:,i,2,v2,ke_xy)
1013 do j=i+1, 4
1014 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
1015 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
1016 end do
1017 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
1018 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
1019 end do
1020 else
1021 do i=1, 4-1
1022 do v=1, 16
1023 ip = ipiv(v,i)
1024 if (ip /= i) then
1025 tmp = b(v,i,1,v2,ke_xy)
1026 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1027 b(v,ip,1,v2,ke_xy) = tmp
1028 tmp = b(v,i,2,v2,ke_xy)
1029 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1030 b(v,ip,2,v2,ke_xy) = tmp
1031
1032 tmp = g(v,i,v2,ke_xy)
1033 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
1034 g(v,ip,v2,ke_xy) = tmp
1035 end if
1036 end do
1037 end do
1038 do i=1, 4
1039 do v=1,16
1040 bb(v,1,i) = b(v,i,1,v2,ke_xy)
1041 bb(v,2,i) = b(v,i,2,v2,ke_xy)
1042 bb(v,3,i) = g(v,i,v2,ke_xy)
1043 end do
1044 end do
1045 do i=1, 4
1046 tmpb(:,:) = bb(:,:,i)
1047 do j=1, i-1
1048 do r=1, 3
1049 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1050 end do
1051 end do
1052 bb(:,:,i) = tmpb(:,:)
1053 end do
1054 do i=4, 1, -1
1055 tmpb(:,:) = bb(:,:,i)
1056 do j=i+1, 4
1057 do r=1, 3
1058 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1059 end do
1060 end do
1061 do r=1, 3
1062 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1063 end do
1064 end do
1065 do i=1, 4
1066 b(:,i,1,v2,ke_xy) = bb(:,1,i)
1067 b(:,i,2,v2,ke_xy) = bb(:,2,i)
1068 g(:,i,v2,ke_xy) = bb(:,3,i)
1069 end do
1070 end if
1071 end do
1072 end do
1073 return
1074 end subroutine solve_nnode4_uv
1075!OCL SERIAL
1076 subroutine solve_nnode4_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
1077 implicit none
1078 integer, intent(in) :: nnode_v
1079 integer, intent(in) :: jm
1080 integer, intent(in) :: ne2d
1081 real(rp), intent(in) :: d(16,12,12,jm,ne2d)
1082 real(rp), intent(inout) :: b(16,12,jm,ne2d)
1083 real(rp), intent(inout) :: g(16,12,3,jm,ne2d)
1084 logical, intent(in) :: is_top
1085
1086 integer :: ke_xy
1087 integer :: v2
1088 integer :: k, i, j, v, r
1089 integer :: ip
1090 integer :: ipiv(16,12)
1091 real(rp) :: tmpv(16)
1092 real(rp) :: a(16,12,12)
1093 real(rp) :: bb(16,4,12)
1094 real(rp) :: tmpb(16,4)
1095 real(rp) :: tmp
1096 !------------------------------
1097
1098 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
1099 do ke_xy=1, ne2d
1100 do v2=1, jm
1101 a(:,:,:) = d(:,:,:,v2,ke_xy)
1102 do k=1, 12
1103 tmpv(:) = abs(a(:,k,k))
1104 ipiv(:,k) = k
1105 do i=k+1, 12
1106 do v=1, 16
1107 tmp = abs(a(v,i,k))
1108 if ( tmp > tmpv(v) ) then
1109 ipiv(v,k) = i
1110 tmpv(v) = tmp
1111 end if
1112 end do
1113 end do
1114 do v=1, 16
1115 ip = ipiv(v,k)
1116 if ( ip /= k ) then
1117 do j=1, 12
1118 tmp = a(v,k,j)
1119 a(v,k,j) = a(v,ip,j)
1120 a(v,ip,j) = tmp
1121 end do
1122 end if
1123 end do
1124 !--
1125 tmpv(:) = 1.0_rp / a(:,k,k)
1126 a(:,k,k) = tmpv(:)
1127 do i=k+1, 12
1128 do v=1, 16
1129 a(v,i,k) = a(v,i,k) * tmpv(v)
1130 end do
1131 end do
1132 do j=k+1, 12
1133 do i=k+1, 12
1134 do v=1, 16
1135 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1136 end do
1137 end do
1138 end do
1139 end do
1140 !---------------------------
1141 if ( is_top ) then
1142 do i=1, 12-1
1143 do v=1, 16
1144 ip = ipiv(v,i)
1145 if (ip /= i) then
1146 tmp = b(v,i,v2,ke_xy)
1147 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1148 b(v,ip,v2,ke_xy) = tmp
1149 end if
1150 end do
1151 end do
1152 do i=1, 12
1153 tmpv(:) = b(:,i,v2,ke_xy)
1154 do j=1, i-1
1155 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
1156 end do
1157 b(:,i,v2,ke_xy) = tmpv(:)
1158 end do
1159 do i=12, 1, -1
1160 tmpv(:) = b(:,i,v2,ke_xy)
1161 do j=i+1, 12
1162 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
1163 end do
1164 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
1165 end do
1166 else
1167 do i=1, 12-1
1168 do v=1, 16
1169 ip = ipiv(v,i)
1170 if (ip /= i) then
1171 tmp = b(v,i,v2,ke_xy)
1172 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1173 b(v,ip,v2,ke_xy) = tmp
1174
1175 tmp = g(v,i,1,v2,ke_xy)
1176 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
1177 g(v,ip,1,v2,ke_xy) = tmp
1178 tmp = g(v,i,2,v2,ke_xy)
1179 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
1180 g(v,ip,2,v2,ke_xy) = tmp
1181 tmp = g(v,i,3,v2,ke_xy)
1182 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
1183 g(v,ip,3,v2,ke_xy) = tmp
1184 end if
1185 end do
1186 end do
1187 do i=1, 12
1188 do v=1,16
1189 bb(v,1,i) = b(v,i,v2,ke_xy)
1190 bb(v,2,i) = g(v,i,1,v2,ke_xy)
1191 bb(v,3,i) = g(v,i,2,v2,ke_xy)
1192 bb(v,4,i) = g(v,i,3,v2,ke_xy)
1193 end do
1194 end do
1195 do i=1, 12
1196 tmpb(:,:) = bb(:,:,i)
1197 do j=1, i-1
1198 do r=1, 4
1199 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1200 end do
1201 end do
1202 bb(:,:,i) = tmpb(:,:)
1203 end do
1204 do i=12, 1, -1
1205 tmpb(:,:) = bb(:,:,i)
1206 do j=i+1, 12
1207 do r=1, 4
1208 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1209 end do
1210 end do
1211 do r=1, 4
1212 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1213 end do
1214 end do
1215 do i=1, 12
1216 b(:,i,v2,ke_xy) = bb(:,1,i)
1217 g(:,i,1,v2,ke_xy) = bb(:,2,i)
1218 g(:,i,2,v2,ke_xy) = bb(:,3,i)
1219 g(:,i,3,v2,ke_xy) = bb(:,4,i)
1220 end do
1221 end if
1222 end do
1223 end do
1224 return
1225 end subroutine solve_nnode4_var3
1226!OCL SERIAL
1227 subroutine solve_nnode5_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
1228 implicit none
1229 integer, intent(in) :: nnode_v
1230 integer, intent(in) :: jm
1231 integer, intent(in) :: ne2d
1232 real(rp), intent(in) :: d(5,5,5,jm,ne2d)
1233 real(rp), intent(inout) :: b(5,5,2,jm,ne2d)
1234 real(rp), intent(inout) :: g(5,5,jm,ne2d)
1235 logical, intent(in) :: is_top
1236
1237 integer :: ke_xy
1238 integer :: v2
1239 integer :: k, i, j, v, r
1240 integer :: ip
1241 integer :: ipiv(5,5)
1242 real(rp) :: tmpv(5)
1243 real(rp) :: tmpv2(5)
1244 real(rp) :: a(5,5,5)
1245 real(rp) :: bb(5,3,5)
1246 real(rp) :: tmpb(5,3)
1247 real(rp) :: tmp
1248 !------------------------------
1249
1250 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
1251 do ke_xy=1, ne2d
1252 do v2=1, jm
1253 a(:,:,:) = d(:,:,:,v2,ke_xy)
1254 do k=1, 5
1255 tmpv(:) = abs(a(:,k,k))
1256 ipiv(:,k) = k
1257 do i=k+1, 5
1258 do v=1, 5
1259 tmp = abs(a(v,i,k))
1260 if ( tmp > tmpv(v) ) then
1261 ipiv(v,k) = i
1262 tmpv(v) = tmp
1263 end if
1264 end do
1265 end do
1266 do v=1, 5
1267 ip = ipiv(v,k)
1268 if ( ip /= k ) then
1269 do j=1, 5
1270 tmp = a(v,k,j)
1271 a(v,k,j) = a(v,ip,j)
1272 a(v,ip,j) = tmp
1273 end do
1274 end if
1275 end do
1276 !--
1277 tmpv(:) = 1.0_rp / a(:,k,k)
1278 a(:,k,k) = tmpv(:)
1279 do i=k+1, 5
1280 a(:,i,k) = a(:,i,k) * tmpv(:)
1281 end do
1282 do j=k+1, 5
1283 do i=k+1, 5
1284 do v=1, 5
1285 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1286 end do
1287 end do
1288 end do
1289 end do
1290 !---------------------------
1291 if ( is_top ) then
1292 do i=1, 5-1
1293 do v=1, 5
1294 ip = ipiv(v,i)
1295 if (ip /= i) then
1296 tmp = b(v,i,1,v2,ke_xy)
1297 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1298 b(v,ip,1,v2,ke_xy) = tmp
1299 tmp = b(v,i,2,v2,ke_xy)
1300 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1301 b(v,ip,2,v2,ke_xy) = tmp
1302 end if
1303 end do
1304 end do
1305 do i=1, 5
1306 tmpv(:) = b(:,i,1,v2,ke_xy)
1307 tmpv2(:) = b(:,i,2,v2,ke_xy)
1308 do j=1, i-1
1309 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
1310 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
1311 end do
1312 b(:,i,1,v2,ke_xy) = tmpv(:)
1313 b(:,i,2,v2,ke_xy) = tmpv2(:)
1314 end do
1315 do i=5, 1, -1
1316 tmpv(:) = b(:,i,1,v2,ke_xy)
1317 tmpv2(:) = b(:,i,2,v2,ke_xy)
1318 do j=i+1, 5
1319 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
1320 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
1321 end do
1322 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
1323 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
1324 end do
1325 else
1326 do i=1, 5-1
1327 do v=1, 5
1328 ip = ipiv(v,i)
1329 if (ip /= i) then
1330 tmp = b(v,i,1,v2,ke_xy)
1331 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1332 b(v,ip,1,v2,ke_xy) = tmp
1333 tmp = b(v,i,2,v2,ke_xy)
1334 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1335 b(v,ip,2,v2,ke_xy) = tmp
1336
1337 tmp = g(v,i,v2,ke_xy)
1338 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
1339 g(v,ip,v2,ke_xy) = tmp
1340 end if
1341 end do
1342 end do
1343 do i=1, 5
1344 do v=1,5
1345 bb(v,1,i) = b(v,i,1,v2,ke_xy)
1346 bb(v,2,i) = b(v,i,2,v2,ke_xy)
1347 bb(v,3,i) = g(v,i,v2,ke_xy)
1348 end do
1349 end do
1350 do i=1, 5
1351 tmpb(:,:) = bb(:,:,i)
1352 do j=1, i-1
1353 do r=1, 3
1354 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1355 end do
1356 end do
1357 bb(:,:,i) = tmpb(:,:)
1358 end do
1359 do i=5, 1, -1
1360 tmpb(:,:) = bb(:,:,i)
1361 do j=i+1, 5
1362 do r=1, 3
1363 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1364 end do
1365 end do
1366 do r=1, 3
1367 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1368 end do
1369 end do
1370 do i=1, 5
1371 b(:,i,1,v2,ke_xy) = bb(:,1,i)
1372 b(:,i,2,v2,ke_xy) = bb(:,2,i)
1373 g(:,i,v2,ke_xy) = bb(:,3,i)
1374 end do
1375 end if
1376 end do
1377 end do
1378 return
1379 end subroutine solve_nnode5_uv
1380!OCL SERIAL
1381 subroutine solve_nnode5_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
1382 implicit none
1383 integer, intent(in) :: nnode_v
1384 integer, intent(in) :: jm
1385 integer, intent(in) :: ne2d
1386 real(rp), intent(in) :: d(5,15,15,jm,ne2d)
1387 real(rp), intent(inout) :: b(5,15,jm,ne2d)
1388 real(rp), intent(inout) :: g(5,15,3,jm,ne2d)
1389 logical, intent(in) :: is_top
1390
1391 integer :: ke_xy
1392 integer :: v2
1393 integer :: k, i, j, v, r
1394 integer :: ip
1395 integer :: ipiv(5,15)
1396 real(rp) :: tmpv(5)
1397 real(rp) :: a(5,15,15)
1398 real(rp) :: bb(5,4,15)
1399 real(rp) :: tmpb(5,4)
1400 real(rp) :: tmp
1401 !------------------------------
1402
1403 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
1404 do ke_xy=1, ne2d
1405 do v2=1, jm
1406 a(:,:,:) = d(:,:,:,v2,ke_xy)
1407 do k=1, 15
1408 tmpv(:) = abs(a(:,k,k))
1409 ipiv(:,k) = k
1410 do i=k+1, 15
1411 do v=1, 5
1412 tmp = abs(a(v,i,k))
1413 if ( tmp > tmpv(v) ) then
1414 ipiv(v,k) = i
1415 tmpv(v) = tmp
1416 end if
1417 end do
1418 end do
1419 do v=1, 5
1420 ip = ipiv(v,k)
1421 if ( ip /= k ) then
1422 do j=1, 15
1423 tmp = a(v,k,j)
1424 a(v,k,j) = a(v,ip,j)
1425 a(v,ip,j) = tmp
1426 end do
1427 end if
1428 end do
1429 !--
1430 tmpv(:) = 1.0_rp / a(:,k,k)
1431 a(:,k,k) = tmpv(:)
1432 do i=k+1, 15
1433 do v=1, 5
1434 a(v,i,k) = a(v,i,k) * tmpv(v)
1435 end do
1436 end do
1437 do j=k+1, 15
1438 do i=k+1, 15
1439 do v=1, 5
1440 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1441 end do
1442 end do
1443 end do
1444 end do
1445 !---------------------------
1446 if ( is_top ) then
1447 do i=1, 15-1
1448 do v=1, 5
1449 ip = ipiv(v,i)
1450 if (ip /= i) then
1451 tmp = b(v,i,v2,ke_xy)
1452 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1453 b(v,ip,v2,ke_xy) = tmp
1454 end if
1455 end do
1456 end do
1457 do i=1, 15
1458 tmpv(:) = b(:,i,v2,ke_xy)
1459 do j=1, i-1
1460 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
1461 end do
1462 b(:,i,v2,ke_xy) = tmpv(:)
1463 end do
1464 do i=15, 1, -1
1465 tmpv(:) = b(:,i,v2,ke_xy)
1466 do j=i+1, 15
1467 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
1468 end do
1469 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
1470 end do
1471 else
1472 do i=1, 15-1
1473 do v=1, 5
1474 ip = ipiv(v,i)
1475 if (ip /= i) then
1476 tmp = b(v,i,v2,ke_xy)
1477 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1478 b(v,ip,v2,ke_xy) = tmp
1479
1480 tmp = g(v,i,1,v2,ke_xy)
1481 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
1482 g(v,ip,1,v2,ke_xy) = tmp
1483 tmp = g(v,i,2,v2,ke_xy)
1484 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
1485 g(v,ip,2,v2,ke_xy) = tmp
1486 tmp = g(v,i,3,v2,ke_xy)
1487 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
1488 g(v,ip,3,v2,ke_xy) = tmp
1489 end if
1490 end do
1491 end do
1492 do i=1, 15
1493 do v=1,5
1494 bb(v,1,i) = b(v,i,v2,ke_xy)
1495 bb(v,2,i) = g(v,i,1,v2,ke_xy)
1496 bb(v,3,i) = g(v,i,2,v2,ke_xy)
1497 bb(v,4,i) = g(v,i,3,v2,ke_xy)
1498 end do
1499 end do
1500 do i=1, 15
1501 tmpb(:,:) = bb(:,:,i)
1502 do j=1, i-1
1503 do r=1, 4
1504 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1505 end do
1506 end do
1507 bb(:,:,i) = tmpb(:,:)
1508 end do
1509 do i=15, 1, -1
1510 tmpb(:,:) = bb(:,:,i)
1511 do j=i+1, 15
1512 do r=1, 4
1513 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1514 end do
1515 end do
1516 do r=1, 4
1517 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1518 end do
1519 end do
1520 do i=1, 15
1521 b(:,i,v2,ke_xy) = bb(:,1,i)
1522 g(:,i,1,v2,ke_xy) = bb(:,2,i)
1523 g(:,i,2,v2,ke_xy) = bb(:,3,i)
1524 g(:,i,3,v2,ke_xy) = bb(:,4,i)
1525 end do
1526 end if
1527 end do
1528 end do
1529 return
1530 end subroutine solve_nnode5_var3
1531!OCL SERIAL
1532 subroutine solve_nnode6_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
1533 implicit none
1534 integer, intent(in) :: nnode_v
1535 integer, intent(in) :: jm
1536 integer, intent(in) :: ne2d
1537 real(rp), intent(in) :: d(6,6,6,jm,ne2d)
1538 real(rp), intent(inout) :: b(6,6,2,jm,ne2d)
1539 real(rp), intent(inout) :: g(6,6,jm,ne2d)
1540 logical, intent(in) :: is_top
1541
1542 integer :: ke_xy
1543 integer :: v2
1544 integer :: k, i, j, v, r
1545 integer :: ip
1546 integer :: ipiv(6,6)
1547 real(rp) :: tmpv(6)
1548 real(rp) :: tmpv2(6)
1549 real(rp) :: a(6,6,6)
1550 real(rp) :: bb(6,3,6)
1551 real(rp) :: tmpb(6,3)
1552 real(rp) :: tmp
1553 !------------------------------
1554
1555 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
1556 do ke_xy=1, ne2d
1557 do v2=1, jm
1558 a(:,:,:) = d(:,:,:,v2,ke_xy)
1559 do k=1, 6
1560 tmpv(:) = abs(a(:,k,k))
1561 ipiv(:,k) = k
1562 do i=k+1, 6
1563 do v=1, 6
1564 tmp = abs(a(v,i,k))
1565 if ( tmp > tmpv(v) ) then
1566 ipiv(v,k) = i
1567 tmpv(v) = tmp
1568 end if
1569 end do
1570 end do
1571 do v=1, 6
1572 ip = ipiv(v,k)
1573 if ( ip /= k ) then
1574 do j=1, 6
1575 tmp = a(v,k,j)
1576 a(v,k,j) = a(v,ip,j)
1577 a(v,ip,j) = tmp
1578 end do
1579 end if
1580 end do
1581 !--
1582 tmpv(:) = 1.0_rp / a(:,k,k)
1583 a(:,k,k) = tmpv(:)
1584 do i=k+1, 6
1585 a(:,i,k) = a(:,i,k) * tmpv(:)
1586 end do
1587 do j=k+1, 6
1588 do i=k+1, 6
1589 do v=1, 6
1590 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1591 end do
1592 end do
1593 end do
1594 end do
1595 !---------------------------
1596 if ( is_top ) then
1597 do i=1, 6-1
1598 do v=1, 6
1599 ip = ipiv(v,i)
1600 if (ip /= i) then
1601 tmp = b(v,i,1,v2,ke_xy)
1602 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1603 b(v,ip,1,v2,ke_xy) = tmp
1604 tmp = b(v,i,2,v2,ke_xy)
1605 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1606 b(v,ip,2,v2,ke_xy) = tmp
1607 end if
1608 end do
1609 end do
1610 do i=1, 6
1611 tmpv(:) = b(:,i,1,v2,ke_xy)
1612 tmpv2(:) = b(:,i,2,v2,ke_xy)
1613 do j=1, i-1
1614 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
1615 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
1616 end do
1617 b(:,i,1,v2,ke_xy) = tmpv(:)
1618 b(:,i,2,v2,ke_xy) = tmpv2(:)
1619 end do
1620 do i=6, 1, -1
1621 tmpv(:) = b(:,i,1,v2,ke_xy)
1622 tmpv2(:) = b(:,i,2,v2,ke_xy)
1623 do j=i+1, 6
1624 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
1625 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
1626 end do
1627 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
1628 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
1629 end do
1630 else
1631 do i=1, 6-1
1632 do v=1, 6
1633 ip = ipiv(v,i)
1634 if (ip /= i) then
1635 tmp = b(v,i,1,v2,ke_xy)
1636 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1637 b(v,ip,1,v2,ke_xy) = tmp
1638 tmp = b(v,i,2,v2,ke_xy)
1639 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1640 b(v,ip,2,v2,ke_xy) = tmp
1641
1642 tmp = g(v,i,v2,ke_xy)
1643 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
1644 g(v,ip,v2,ke_xy) = tmp
1645 end if
1646 end do
1647 end do
1648 do i=1, 6
1649 do v=1,6
1650 bb(v,1,i) = b(v,i,1,v2,ke_xy)
1651 bb(v,2,i) = b(v,i,2,v2,ke_xy)
1652 bb(v,3,i) = g(v,i,v2,ke_xy)
1653 end do
1654 end do
1655 do i=1, 6
1656 tmpb(:,:) = bb(:,:,i)
1657 do j=1, i-1
1658 do r=1, 3
1659 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1660 end do
1661 end do
1662 bb(:,:,i) = tmpb(:,:)
1663 end do
1664 do i=6, 1, -1
1665 tmpb(:,:) = bb(:,:,i)
1666 do j=i+1, 6
1667 do r=1, 3
1668 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1669 end do
1670 end do
1671 do r=1, 3
1672 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1673 end do
1674 end do
1675 do i=1, 6
1676 b(:,i,1,v2,ke_xy) = bb(:,1,i)
1677 b(:,i,2,v2,ke_xy) = bb(:,2,i)
1678 g(:,i,v2,ke_xy) = bb(:,3,i)
1679 end do
1680 end if
1681 end do
1682 end do
1683 return
1684 end subroutine solve_nnode6_uv
1685!OCL SERIAL
1686 subroutine solve_nnode6_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
1687 implicit none
1688 integer, intent(in) :: nnode_v
1689 integer, intent(in) :: jm
1690 integer, intent(in) :: ne2d
1691 real(rp), intent(in) :: d(6,18,18,jm,ne2d)
1692 real(rp), intent(inout) :: b(6,18,jm,ne2d)
1693 real(rp), intent(inout) :: g(6,18,3,jm,ne2d)
1694 logical, intent(in) :: is_top
1695
1696 integer :: ke_xy
1697 integer :: v2
1698 integer :: k, i, j, v, r
1699 integer :: ip
1700 integer :: ipiv(6,18)
1701 real(rp) :: tmpv(6)
1702 real(rp) :: a(6,18,18)
1703 real(rp) :: bb(6,4,18)
1704 real(rp) :: tmpb(6,4)
1705 real(rp) :: tmp
1706 !------------------------------
1707
1708 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
1709 do ke_xy=1, ne2d
1710 do v2=1, jm
1711 a(:,:,:) = d(:,:,:,v2,ke_xy)
1712 do k=1, 18
1713 tmpv(:) = abs(a(:,k,k))
1714 ipiv(:,k) = k
1715 do i=k+1, 18
1716 do v=1, 6
1717 tmp = abs(a(v,i,k))
1718 if ( tmp > tmpv(v) ) then
1719 ipiv(v,k) = i
1720 tmpv(v) = tmp
1721 end if
1722 end do
1723 end do
1724 do v=1, 6
1725 ip = ipiv(v,k)
1726 if ( ip /= k ) then
1727 do j=1, 18
1728 tmp = a(v,k,j)
1729 a(v,k,j) = a(v,ip,j)
1730 a(v,ip,j) = tmp
1731 end do
1732 end if
1733 end do
1734 !--
1735 tmpv(:) = 1.0_rp / a(:,k,k)
1736 a(:,k,k) = tmpv(:)
1737 do i=k+1, 18
1738 do v=1, 6
1739 a(v,i,k) = a(v,i,k) * tmpv(v)
1740 end do
1741 end do
1742 do j=k+1, 18
1743 do i=k+1, 18
1744 do v=1, 6
1745 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1746 end do
1747 end do
1748 end do
1749 end do
1750 !---------------------------
1751 if ( is_top ) then
1752 do i=1, 18-1
1753 do v=1, 6
1754 ip = ipiv(v,i)
1755 if (ip /= i) then
1756 tmp = b(v,i,v2,ke_xy)
1757 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1758 b(v,ip,v2,ke_xy) = tmp
1759 end if
1760 end do
1761 end do
1762 do i=1, 18
1763 tmpv(:) = b(:,i,v2,ke_xy)
1764 do j=1, i-1
1765 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
1766 end do
1767 b(:,i,v2,ke_xy) = tmpv(:)
1768 end do
1769 do i=18, 1, -1
1770 tmpv(:) = b(:,i,v2,ke_xy)
1771 do j=i+1, 18
1772 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
1773 end do
1774 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
1775 end do
1776 else
1777 do i=1, 18-1
1778 do v=1, 6
1779 ip = ipiv(v,i)
1780 if (ip /= i) then
1781 tmp = b(v,i,v2,ke_xy)
1782 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
1783 b(v,ip,v2,ke_xy) = tmp
1784
1785 tmp = g(v,i,1,v2,ke_xy)
1786 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
1787 g(v,ip,1,v2,ke_xy) = tmp
1788 tmp = g(v,i,2,v2,ke_xy)
1789 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
1790 g(v,ip,2,v2,ke_xy) = tmp
1791 tmp = g(v,i,3,v2,ke_xy)
1792 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
1793 g(v,ip,3,v2,ke_xy) = tmp
1794 end if
1795 end do
1796 end do
1797 do i=1, 18
1798 do v=1,6
1799 bb(v,1,i) = b(v,i,v2,ke_xy)
1800 bb(v,2,i) = g(v,i,1,v2,ke_xy)
1801 bb(v,3,i) = g(v,i,2,v2,ke_xy)
1802 bb(v,4,i) = g(v,i,3,v2,ke_xy)
1803 end do
1804 end do
1805 do i=1, 18
1806 tmpb(:,:) = bb(:,:,i)
1807 do j=1, i-1
1808 do r=1, 4
1809 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1810 end do
1811 end do
1812 bb(:,:,i) = tmpb(:,:)
1813 end do
1814 do i=18, 1, -1
1815 tmpb(:,:) = bb(:,:,i)
1816 do j=i+1, 18
1817 do r=1, 4
1818 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1819 end do
1820 end do
1821 do r=1, 4
1822 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1823 end do
1824 end do
1825 do i=1, 18
1826 b(:,i,v2,ke_xy) = bb(:,1,i)
1827 g(:,i,1,v2,ke_xy) = bb(:,2,i)
1828 g(:,i,2,v2,ke_xy) = bb(:,3,i)
1829 g(:,i,3,v2,ke_xy) = bb(:,4,i)
1830 end do
1831 end if
1832 end do
1833 end do
1834 return
1835 end subroutine solve_nnode6_var3
1836!OCL SERIAL
1837 subroutine solve_nnode7_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
1838 implicit none
1839 integer, intent(in) :: nnode_v
1840 integer, intent(in) :: jm
1841 integer, intent(in) :: ne2d
1842 real(rp), intent(in) :: d(7,7,7,jm,ne2d)
1843 real(rp), intent(inout) :: b(7,7,2,jm,ne2d)
1844 real(rp), intent(inout) :: g(7,7,jm,ne2d)
1845 logical, intent(in) :: is_top
1846
1847 integer :: ke_xy
1848 integer :: v2
1849 integer :: k, i, j, v, r
1850 integer :: ip
1851 integer :: ipiv(7,7)
1852 real(rp) :: tmpv(7)
1853 real(rp) :: tmpv2(7)
1854 real(rp) :: a(7,7,7)
1855 real(rp) :: bb(7,3,7)
1856 real(rp) :: tmpb(7,3)
1857 real(rp) :: tmp
1858 !------------------------------
1859
1860 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
1861 do ke_xy=1, ne2d
1862 do v2=1, jm
1863 a(:,:,:) = d(:,:,:,v2,ke_xy)
1864 do k=1, 7
1865 tmpv(:) = abs(a(:,k,k))
1866 ipiv(:,k) = k
1867 do i=k+1, 7
1868 do v=1, 7
1869 tmp = abs(a(v,i,k))
1870 if ( tmp > tmpv(v) ) then
1871 ipiv(v,k) = i
1872 tmpv(v) = tmp
1873 end if
1874 end do
1875 end do
1876 do v=1, 7
1877 ip = ipiv(v,k)
1878 if ( ip /= k ) then
1879 do j=1, 7
1880 tmp = a(v,k,j)
1881 a(v,k,j) = a(v,ip,j)
1882 a(v,ip,j) = tmp
1883 end do
1884 end if
1885 end do
1886 !--
1887 tmpv(:) = 1.0_rp / a(:,k,k)
1888 a(:,k,k) = tmpv(:)
1889 do i=k+1, 7
1890 a(:,i,k) = a(:,i,k) * tmpv(:)
1891 end do
1892 do j=k+1, 7
1893 do i=k+1, 7
1894 do v=1, 7
1895 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
1896 end do
1897 end do
1898 end do
1899 end do
1900 !---------------------------
1901 if ( is_top ) then
1902 do i=1, 7-1
1903 do v=1, 7
1904 ip = ipiv(v,i)
1905 if (ip /= i) then
1906 tmp = b(v,i,1,v2,ke_xy)
1907 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1908 b(v,ip,1,v2,ke_xy) = tmp
1909 tmp = b(v,i,2,v2,ke_xy)
1910 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1911 b(v,ip,2,v2,ke_xy) = tmp
1912 end if
1913 end do
1914 end do
1915 do i=1, 7
1916 tmpv(:) = b(:,i,1,v2,ke_xy)
1917 tmpv2(:) = b(:,i,2,v2,ke_xy)
1918 do j=1, i-1
1919 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
1920 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
1921 end do
1922 b(:,i,1,v2,ke_xy) = tmpv(:)
1923 b(:,i,2,v2,ke_xy) = tmpv2(:)
1924 end do
1925 do i=7, 1, -1
1926 tmpv(:) = b(:,i,1,v2,ke_xy)
1927 tmpv2(:) = b(:,i,2,v2,ke_xy)
1928 do j=i+1, 7
1929 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
1930 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
1931 end do
1932 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
1933 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
1934 end do
1935 else
1936 do i=1, 7-1
1937 do v=1, 7
1938 ip = ipiv(v,i)
1939 if (ip /= i) then
1940 tmp = b(v,i,1,v2,ke_xy)
1941 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
1942 b(v,ip,1,v2,ke_xy) = tmp
1943 tmp = b(v,i,2,v2,ke_xy)
1944 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
1945 b(v,ip,2,v2,ke_xy) = tmp
1946
1947 tmp = g(v,i,v2,ke_xy)
1948 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
1949 g(v,ip,v2,ke_xy) = tmp
1950 end if
1951 end do
1952 end do
1953 do i=1, 7
1954 do v=1,7
1955 bb(v,1,i) = b(v,i,1,v2,ke_xy)
1956 bb(v,2,i) = b(v,i,2,v2,ke_xy)
1957 bb(v,3,i) = g(v,i,v2,ke_xy)
1958 end do
1959 end do
1960 do i=1, 7
1961 tmpb(:,:) = bb(:,:,i)
1962 do j=1, i-1
1963 do r=1, 3
1964 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
1965 end do
1966 end do
1967 bb(:,:,i) = tmpb(:,:)
1968 end do
1969 do i=7, 1, -1
1970 tmpb(:,:) = bb(:,:,i)
1971 do j=i+1, 7
1972 do r=1, 3
1973 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
1974 end do
1975 end do
1976 do r=1, 3
1977 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
1978 end do
1979 end do
1980 do i=1, 7
1981 b(:,i,1,v2,ke_xy) = bb(:,1,i)
1982 b(:,i,2,v2,ke_xy) = bb(:,2,i)
1983 g(:,i,v2,ke_xy) = bb(:,3,i)
1984 end do
1985 end if
1986 end do
1987 end do
1988 return
1989 end subroutine solve_nnode7_uv
1990!OCL SERIAL
1991 subroutine solve_nnode7_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
1992 implicit none
1993 integer, intent(in) :: nnode_v
1994 integer, intent(in) :: jm
1995 integer, intent(in) :: ne2d
1996 real(rp), intent(in) :: d(7,21,21,jm,ne2d)
1997 real(rp), intent(inout) :: b(7,21,jm,ne2d)
1998 real(rp), intent(inout) :: g(7,21,3,jm,ne2d)
1999 logical, intent(in) :: is_top
2000
2001 integer :: ke_xy
2002 integer :: v2
2003 integer :: k, i, j, v, r
2004 integer :: ip
2005 integer :: ipiv(7,21)
2006 real(rp) :: tmpv(7)
2007 real(rp) :: a(7,21,21)
2008 real(rp) :: bb(7,4,21)
2009 real(rp) :: tmpb(7,4)
2010 real(rp) :: tmp
2011 !------------------------------
2012
2013 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
2014 do ke_xy=1, ne2d
2015 do v2=1, jm
2016 a(:,:,:) = d(:,:,:,v2,ke_xy)
2017 do k=1, 21
2018 tmpv(:) = abs(a(:,k,k))
2019 ipiv(:,k) = k
2020 do i=k+1, 21
2021 do v=1, 7
2022 tmp = abs(a(v,i,k))
2023 if ( tmp > tmpv(v) ) then
2024 ipiv(v,k) = i
2025 tmpv(v) = tmp
2026 end if
2027 end do
2028 end do
2029 do v=1, 7
2030 ip = ipiv(v,k)
2031 if ( ip /= k ) then
2032 do j=1, 21
2033 tmp = a(v,k,j)
2034 a(v,k,j) = a(v,ip,j)
2035 a(v,ip,j) = tmp
2036 end do
2037 end if
2038 end do
2039 !--
2040 tmpv(:) = 1.0_rp / a(:,k,k)
2041 a(:,k,k) = tmpv(:)
2042 do i=k+1, 21
2043 do v=1, 7
2044 a(v,i,k) = a(v,i,k) * tmpv(v)
2045 end do
2046 end do
2047 do j=k+1, 21
2048 do i=k+1, 21
2049 do v=1, 7
2050 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2051 end do
2052 end do
2053 end do
2054 end do
2055 !---------------------------
2056 if ( is_top ) then
2057 do i=1, 21-1
2058 do v=1, 7
2059 ip = ipiv(v,i)
2060 if (ip /= i) then
2061 tmp = b(v,i,v2,ke_xy)
2062 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2063 b(v,ip,v2,ke_xy) = tmp
2064 end if
2065 end do
2066 end do
2067 do i=1, 21
2068 tmpv(:) = b(:,i,v2,ke_xy)
2069 do j=1, i-1
2070 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
2071 end do
2072 b(:,i,v2,ke_xy) = tmpv(:)
2073 end do
2074 do i=21, 1, -1
2075 tmpv(:) = b(:,i,v2,ke_xy)
2076 do j=i+1, 21
2077 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
2078 end do
2079 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
2080 end do
2081 else
2082 do i=1, 21-1
2083 do v=1, 7
2084 ip = ipiv(v,i)
2085 if (ip /= i) then
2086 tmp = b(v,i,v2,ke_xy)
2087 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2088 b(v,ip,v2,ke_xy) = tmp
2089
2090 tmp = g(v,i,1,v2,ke_xy)
2091 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
2092 g(v,ip,1,v2,ke_xy) = tmp
2093 tmp = g(v,i,2,v2,ke_xy)
2094 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
2095 g(v,ip,2,v2,ke_xy) = tmp
2096 tmp = g(v,i,3,v2,ke_xy)
2097 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
2098 g(v,ip,3,v2,ke_xy) = tmp
2099 end if
2100 end do
2101 end do
2102 do i=1, 21
2103 do v=1,7
2104 bb(v,1,i) = b(v,i,v2,ke_xy)
2105 bb(v,2,i) = g(v,i,1,v2,ke_xy)
2106 bb(v,3,i) = g(v,i,2,v2,ke_xy)
2107 bb(v,4,i) = g(v,i,3,v2,ke_xy)
2108 end do
2109 end do
2110 do i=1, 21
2111 tmpb(:,:) = bb(:,:,i)
2112 do j=1, i-1
2113 do r=1, 4
2114 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2115 end do
2116 end do
2117 bb(:,:,i) = tmpb(:,:)
2118 end do
2119 do i=21, 1, -1
2120 tmpb(:,:) = bb(:,:,i)
2121 do j=i+1, 21
2122 do r=1, 4
2123 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2124 end do
2125 end do
2126 do r=1, 4
2127 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2128 end do
2129 end do
2130 do i=1, 21
2131 b(:,i,v2,ke_xy) = bb(:,1,i)
2132 g(:,i,1,v2,ke_xy) = bb(:,2,i)
2133 g(:,i,2,v2,ke_xy) = bb(:,3,i)
2134 g(:,i,3,v2,ke_xy) = bb(:,4,i)
2135 end do
2136 end if
2137 end do
2138 end do
2139 return
2140 end subroutine solve_nnode7_var3
2141!OCL SERIAL
2142 subroutine solve_nnode8_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
2143 implicit none
2144 integer, intent(in) :: nnode_v
2145 integer, intent(in) :: jm
2146 integer, intent(in) :: ne2d
2147 real(rp), intent(in) :: d(8,8,8,jm,ne2d)
2148 real(rp), intent(inout) :: b(8,8,2,jm,ne2d)
2149 real(rp), intent(inout) :: g(8,8,jm,ne2d)
2150 logical, intent(in) :: is_top
2151
2152 integer :: ke_xy
2153 integer :: v2
2154 integer :: k, i, j, v, r
2155 integer :: ip
2156 integer :: ipiv(8,8)
2157 real(rp) :: tmpv(8)
2158 real(rp) :: tmpv2(8)
2159 real(rp) :: a(8,8,8)
2160 real(rp) :: bb(8,3,8)
2161 real(rp) :: tmpb(8,3)
2162 real(rp) :: tmp
2163 !------------------------------
2164
2165 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
2166 do ke_xy=1, ne2d
2167 do v2=1, jm
2168 a(:,:,:) = d(:,:,:,v2,ke_xy)
2169 do k=1, 8
2170 tmpv(:) = abs(a(:,k,k))
2171 ipiv(:,k) = k
2172 do i=k+1, 8
2173 do v=1, 8
2174 tmp = abs(a(v,i,k))
2175 if ( tmp > tmpv(v) ) then
2176 ipiv(v,k) = i
2177 tmpv(v) = tmp
2178 end if
2179 end do
2180 end do
2181 do v=1, 8
2182 ip = ipiv(v,k)
2183 if ( ip /= k ) then
2184 do j=1, 8
2185 tmp = a(v,k,j)
2186 a(v,k,j) = a(v,ip,j)
2187 a(v,ip,j) = tmp
2188 end do
2189 end if
2190 end do
2191 !--
2192 tmpv(:) = 1.0_rp / a(:,k,k)
2193 a(:,k,k) = tmpv(:)
2194 do i=k+1, 8
2195 a(:,i,k) = a(:,i,k) * tmpv(:)
2196 end do
2197 do j=k+1, 8
2198 do i=k+1, 8
2199 do v=1, 8
2200 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2201 end do
2202 end do
2203 end do
2204 end do
2205 !---------------------------
2206 if ( is_top ) then
2207 do i=1, 8-1
2208 do v=1, 8
2209 ip = ipiv(v,i)
2210 if (ip /= i) then
2211 tmp = b(v,i,1,v2,ke_xy)
2212 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2213 b(v,ip,1,v2,ke_xy) = tmp
2214 tmp = b(v,i,2,v2,ke_xy)
2215 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2216 b(v,ip,2,v2,ke_xy) = tmp
2217 end if
2218 end do
2219 end do
2220 do i=1, 8
2221 tmpv(:) = b(:,i,1,v2,ke_xy)
2222 tmpv2(:) = b(:,i,2,v2,ke_xy)
2223 do j=1, i-1
2224 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
2225 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
2226 end do
2227 b(:,i,1,v2,ke_xy) = tmpv(:)
2228 b(:,i,2,v2,ke_xy) = tmpv2(:)
2229 end do
2230 do i=8, 1, -1
2231 tmpv(:) = b(:,i,1,v2,ke_xy)
2232 tmpv2(:) = b(:,i,2,v2,ke_xy)
2233 do j=i+1, 8
2234 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
2235 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
2236 end do
2237 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
2238 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
2239 end do
2240 else
2241 do i=1, 8-1
2242 do v=1, 8
2243 ip = ipiv(v,i)
2244 if (ip /= i) then
2245 tmp = b(v,i,1,v2,ke_xy)
2246 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2247 b(v,ip,1,v2,ke_xy) = tmp
2248 tmp = b(v,i,2,v2,ke_xy)
2249 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2250 b(v,ip,2,v2,ke_xy) = tmp
2251
2252 tmp = g(v,i,v2,ke_xy)
2253 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
2254 g(v,ip,v2,ke_xy) = tmp
2255 end if
2256 end do
2257 end do
2258 do i=1, 8
2259 do v=1,8
2260 bb(v,1,i) = b(v,i,1,v2,ke_xy)
2261 bb(v,2,i) = b(v,i,2,v2,ke_xy)
2262 bb(v,3,i) = g(v,i,v2,ke_xy)
2263 end do
2264 end do
2265 do i=1, 8
2266 tmpb(:,:) = bb(:,:,i)
2267 do j=1, i-1
2268 do r=1, 3
2269 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2270 end do
2271 end do
2272 bb(:,:,i) = tmpb(:,:)
2273 end do
2274 do i=8, 1, -1
2275 tmpb(:,:) = bb(:,:,i)
2276 do j=i+1, 8
2277 do r=1, 3
2278 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2279 end do
2280 end do
2281 do r=1, 3
2282 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2283 end do
2284 end do
2285 do i=1, 8
2286 b(:,i,1,v2,ke_xy) = bb(:,1,i)
2287 b(:,i,2,v2,ke_xy) = bb(:,2,i)
2288 g(:,i,v2,ke_xy) = bb(:,3,i)
2289 end do
2290 end if
2291 end do
2292 end do
2293 return
2294 end subroutine solve_nnode8_uv
2295!OCL SERIAL
2296 subroutine solve_nnode8_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
2297 implicit none
2298 integer, intent(in) :: nnode_v
2299 integer, intent(in) :: jm
2300 integer, intent(in) :: ne2d
2301 real(rp), intent(in) :: d(8,24,24,jm,ne2d)
2302 real(rp), intent(inout) :: b(8,24,jm,ne2d)
2303 real(rp), intent(inout) :: g(8,24,3,jm,ne2d)
2304 logical, intent(in) :: is_top
2305
2306 integer :: ke_xy
2307 integer :: v2
2308 integer :: k, i, j, v, r
2309 integer :: ip
2310 integer :: ipiv(8,24)
2311 real(rp) :: tmpv(8)
2312 real(rp) :: a(8,24,24)
2313 real(rp) :: bb(8,4,24)
2314 real(rp) :: tmpb(8,4)
2315 real(rp) :: tmp
2316 !------------------------------
2317
2318 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
2319 do ke_xy=1, ne2d
2320 do v2=1, jm
2321 a(:,:,:) = d(:,:,:,v2,ke_xy)
2322 do k=1, 24
2323 tmpv(:) = abs(a(:,k,k))
2324 ipiv(:,k) = k
2325 do i=k+1, 24
2326 do v=1, 8
2327 tmp = abs(a(v,i,k))
2328 if ( tmp > tmpv(v) ) then
2329 ipiv(v,k) = i
2330 tmpv(v) = tmp
2331 end if
2332 end do
2333 end do
2334 do v=1, 8
2335 ip = ipiv(v,k)
2336 if ( ip /= k ) then
2337 do j=1, 24
2338 tmp = a(v,k,j)
2339 a(v,k,j) = a(v,ip,j)
2340 a(v,ip,j) = tmp
2341 end do
2342 end if
2343 end do
2344 !--
2345 tmpv(:) = 1.0_rp / a(:,k,k)
2346 a(:,k,k) = tmpv(:)
2347 do i=k+1, 24
2348 do v=1, 8
2349 a(v,i,k) = a(v,i,k) * tmpv(v)
2350 end do
2351 end do
2352 do j=k+1, 24
2353 do i=k+1, 24
2354 do v=1, 8
2355 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2356 end do
2357 end do
2358 end do
2359 end do
2360 !---------------------------
2361 if ( is_top ) then
2362 do i=1, 24-1
2363 do v=1, 8
2364 ip = ipiv(v,i)
2365 if (ip /= i) then
2366 tmp = b(v,i,v2,ke_xy)
2367 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2368 b(v,ip,v2,ke_xy) = tmp
2369 end if
2370 end do
2371 end do
2372 do i=1, 24
2373 tmpv(:) = b(:,i,v2,ke_xy)
2374 do j=1, i-1
2375 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
2376 end do
2377 b(:,i,v2,ke_xy) = tmpv(:)
2378 end do
2379 do i=24, 1, -1
2380 tmpv(:) = b(:,i,v2,ke_xy)
2381 do j=i+1, 24
2382 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
2383 end do
2384 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
2385 end do
2386 else
2387 do i=1, 24-1
2388 do v=1, 8
2389 ip = ipiv(v,i)
2390 if (ip /= i) then
2391 tmp = b(v,i,v2,ke_xy)
2392 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2393 b(v,ip,v2,ke_xy) = tmp
2394
2395 tmp = g(v,i,1,v2,ke_xy)
2396 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
2397 g(v,ip,1,v2,ke_xy) = tmp
2398 tmp = g(v,i,2,v2,ke_xy)
2399 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
2400 g(v,ip,2,v2,ke_xy) = tmp
2401 tmp = g(v,i,3,v2,ke_xy)
2402 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
2403 g(v,ip,3,v2,ke_xy) = tmp
2404 end if
2405 end do
2406 end do
2407 do i=1, 24
2408 do v=1,8
2409 bb(v,1,i) = b(v,i,v2,ke_xy)
2410 bb(v,2,i) = g(v,i,1,v2,ke_xy)
2411 bb(v,3,i) = g(v,i,2,v2,ke_xy)
2412 bb(v,4,i) = g(v,i,3,v2,ke_xy)
2413 end do
2414 end do
2415 do i=1, 24
2416 tmpb(:,:) = bb(:,:,i)
2417 do j=1, i-1
2418 do r=1, 4
2419 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2420 end do
2421 end do
2422 bb(:,:,i) = tmpb(:,:)
2423 end do
2424 do i=24, 1, -1
2425 tmpb(:,:) = bb(:,:,i)
2426 do j=i+1, 24
2427 do r=1, 4
2428 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2429 end do
2430 end do
2431 do r=1, 4
2432 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2433 end do
2434 end do
2435 do i=1, 24
2436 b(:,i,v2,ke_xy) = bb(:,1,i)
2437 g(:,i,1,v2,ke_xy) = bb(:,2,i)
2438 g(:,i,2,v2,ke_xy) = bb(:,3,i)
2439 g(:,i,3,v2,ke_xy) = bb(:,4,i)
2440 end do
2441 end if
2442 end do
2443 end do
2444 return
2445 end subroutine solve_nnode8_var3
2446!OCL SERIAL
2447 subroutine solve_nnode9_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
2448 implicit none
2449 integer, intent(in) :: nnode_v
2450 integer, intent(in) :: jm
2451 integer, intent(in) :: ne2d
2452 real(rp), intent(in) :: d(9,9,9,jm,ne2d)
2453 real(rp), intent(inout) :: b(9,9,2,jm,ne2d)
2454 real(rp), intent(inout) :: g(9,9,jm,ne2d)
2455 logical, intent(in) :: is_top
2456
2457 integer :: ke_xy
2458 integer :: v2
2459 integer :: k, i, j, v, r
2460 integer :: ip
2461 integer :: ipiv(9,9)
2462 real(rp) :: tmpv(9)
2463 real(rp) :: tmpv2(9)
2464 real(rp) :: a(9,9,9)
2465 real(rp) :: bb(9,3,9)
2466 real(rp) :: tmpb(9,3)
2467 real(rp) :: tmp
2468 !------------------------------
2469
2470 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
2471 do ke_xy=1, ne2d
2472 do v2=1, jm
2473 a(:,:,:) = d(:,:,:,v2,ke_xy)
2474 do k=1, 9
2475 tmpv(:) = abs(a(:,k,k))
2476 ipiv(:,k) = k
2477 do i=k+1, 9
2478 do v=1, 9
2479 tmp = abs(a(v,i,k))
2480 if ( tmp > tmpv(v) ) then
2481 ipiv(v,k) = i
2482 tmpv(v) = tmp
2483 end if
2484 end do
2485 end do
2486 do v=1, 9
2487 ip = ipiv(v,k)
2488 if ( ip /= k ) then
2489 do j=1, 9
2490 tmp = a(v,k,j)
2491 a(v,k,j) = a(v,ip,j)
2492 a(v,ip,j) = tmp
2493 end do
2494 end if
2495 end do
2496 !--
2497 tmpv(:) = 1.0_rp / a(:,k,k)
2498 a(:,k,k) = tmpv(:)
2499 do i=k+1, 9
2500 a(:,i,k) = a(:,i,k) * tmpv(:)
2501 end do
2502 do j=k+1, 9
2503 do i=k+1, 9
2504 do v=1, 9
2505 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2506 end do
2507 end do
2508 end do
2509 end do
2510 !---------------------------
2511 if ( is_top ) then
2512 do i=1, 9-1
2513 do v=1, 9
2514 ip = ipiv(v,i)
2515 if (ip /= i) then
2516 tmp = b(v,i,1,v2,ke_xy)
2517 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2518 b(v,ip,1,v2,ke_xy) = tmp
2519 tmp = b(v,i,2,v2,ke_xy)
2520 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2521 b(v,ip,2,v2,ke_xy) = tmp
2522 end if
2523 end do
2524 end do
2525 do i=1, 9
2526 tmpv(:) = b(:,i,1,v2,ke_xy)
2527 tmpv2(:) = b(:,i,2,v2,ke_xy)
2528 do j=1, i-1
2529 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
2530 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
2531 end do
2532 b(:,i,1,v2,ke_xy) = tmpv(:)
2533 b(:,i,2,v2,ke_xy) = tmpv2(:)
2534 end do
2535 do i=9, 1, -1
2536 tmpv(:) = b(:,i,1,v2,ke_xy)
2537 tmpv2(:) = b(:,i,2,v2,ke_xy)
2538 do j=i+1, 9
2539 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
2540 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
2541 end do
2542 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
2543 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
2544 end do
2545 else
2546 do i=1, 9-1
2547 do v=1, 9
2548 ip = ipiv(v,i)
2549 if (ip /= i) then
2550 tmp = b(v,i,1,v2,ke_xy)
2551 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2552 b(v,ip,1,v2,ke_xy) = tmp
2553 tmp = b(v,i,2,v2,ke_xy)
2554 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2555 b(v,ip,2,v2,ke_xy) = tmp
2556
2557 tmp = g(v,i,v2,ke_xy)
2558 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
2559 g(v,ip,v2,ke_xy) = tmp
2560 end if
2561 end do
2562 end do
2563 do i=1, 9
2564 do v=1,9
2565 bb(v,1,i) = b(v,i,1,v2,ke_xy)
2566 bb(v,2,i) = b(v,i,2,v2,ke_xy)
2567 bb(v,3,i) = g(v,i,v2,ke_xy)
2568 end do
2569 end do
2570 do i=1, 9
2571 tmpb(:,:) = bb(:,:,i)
2572 do j=1, i-1
2573 do r=1, 3
2574 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2575 end do
2576 end do
2577 bb(:,:,i) = tmpb(:,:)
2578 end do
2579 do i=9, 1, -1
2580 tmpb(:,:) = bb(:,:,i)
2581 do j=i+1, 9
2582 do r=1, 3
2583 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2584 end do
2585 end do
2586 do r=1, 3
2587 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2588 end do
2589 end do
2590 do i=1, 9
2591 b(:,i,1,v2,ke_xy) = bb(:,1,i)
2592 b(:,i,2,v2,ke_xy) = bb(:,2,i)
2593 g(:,i,v2,ke_xy) = bb(:,3,i)
2594 end do
2595 end if
2596 end do
2597 end do
2598 return
2599 end subroutine solve_nnode9_uv
2600!OCL SERIAL
2601 subroutine solve_nnode9_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
2602 implicit none
2603 integer, intent(in) :: nnode_v
2604 integer, intent(in) :: jm
2605 integer, intent(in) :: ne2d
2606 real(rp), intent(in) :: d(9,27,27,jm,ne2d)
2607 real(rp), intent(inout) :: b(9,27,jm,ne2d)
2608 real(rp), intent(inout) :: g(9,27,3,jm,ne2d)
2609 logical, intent(in) :: is_top
2610
2611 integer :: ke_xy
2612 integer :: v2
2613 integer :: k, i, j, v, r
2614 integer :: ip
2615 integer :: ipiv(9,27)
2616 real(rp) :: tmpv(9)
2617 real(rp) :: a(9,27,27)
2618 real(rp) :: bb(9,4,27)
2619 real(rp) :: tmpb(9,4)
2620 real(rp) :: tmp
2621 !------------------------------
2622
2623 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
2624 do ke_xy=1, ne2d
2625 do v2=1, jm
2626 a(:,:,:) = d(:,:,:,v2,ke_xy)
2627 do k=1, 27
2628 tmpv(:) = abs(a(:,k,k))
2629 ipiv(:,k) = k
2630 do i=k+1, 27
2631 do v=1, 9
2632 tmp = abs(a(v,i,k))
2633 if ( tmp > tmpv(v) ) then
2634 ipiv(v,k) = i
2635 tmpv(v) = tmp
2636 end if
2637 end do
2638 end do
2639 do v=1, 9
2640 ip = ipiv(v,k)
2641 if ( ip /= k ) then
2642 do j=1, 27
2643 tmp = a(v,k,j)
2644 a(v,k,j) = a(v,ip,j)
2645 a(v,ip,j) = tmp
2646 end do
2647 end if
2648 end do
2649 !--
2650 tmpv(:) = 1.0_rp / a(:,k,k)
2651 a(:,k,k) = tmpv(:)
2652 do i=k+1, 27
2653 do v=1, 9
2654 a(v,i,k) = a(v,i,k) * tmpv(v)
2655 end do
2656 end do
2657 do j=k+1, 27
2658 do i=k+1, 27
2659 do v=1, 9
2660 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2661 end do
2662 end do
2663 end do
2664 end do
2665 !---------------------------
2666 if ( is_top ) then
2667 do i=1, 27-1
2668 do v=1, 9
2669 ip = ipiv(v,i)
2670 if (ip /= i) then
2671 tmp = b(v,i,v2,ke_xy)
2672 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2673 b(v,ip,v2,ke_xy) = tmp
2674 end if
2675 end do
2676 end do
2677 do i=1, 27
2678 tmpv(:) = b(:,i,v2,ke_xy)
2679 do j=1, i-1
2680 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
2681 end do
2682 b(:,i,v2,ke_xy) = tmpv(:)
2683 end do
2684 do i=27, 1, -1
2685 tmpv(:) = b(:,i,v2,ke_xy)
2686 do j=i+1, 27
2687 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
2688 end do
2689 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
2690 end do
2691 else
2692 do i=1, 27-1
2693 do v=1, 9
2694 ip = ipiv(v,i)
2695 if (ip /= i) then
2696 tmp = b(v,i,v2,ke_xy)
2697 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2698 b(v,ip,v2,ke_xy) = tmp
2699
2700 tmp = g(v,i,1,v2,ke_xy)
2701 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
2702 g(v,ip,1,v2,ke_xy) = tmp
2703 tmp = g(v,i,2,v2,ke_xy)
2704 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
2705 g(v,ip,2,v2,ke_xy) = tmp
2706 tmp = g(v,i,3,v2,ke_xy)
2707 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
2708 g(v,ip,3,v2,ke_xy) = tmp
2709 end if
2710 end do
2711 end do
2712 do i=1, 27
2713 do v=1,9
2714 bb(v,1,i) = b(v,i,v2,ke_xy)
2715 bb(v,2,i) = g(v,i,1,v2,ke_xy)
2716 bb(v,3,i) = g(v,i,2,v2,ke_xy)
2717 bb(v,4,i) = g(v,i,3,v2,ke_xy)
2718 end do
2719 end do
2720 do i=1, 27
2721 tmpb(:,:) = bb(:,:,i)
2722 do j=1, i-1
2723 do r=1, 4
2724 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2725 end do
2726 end do
2727 bb(:,:,i) = tmpb(:,:)
2728 end do
2729 do i=27, 1, -1
2730 tmpb(:,:) = bb(:,:,i)
2731 do j=i+1, 27
2732 do r=1, 4
2733 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2734 end do
2735 end do
2736 do r=1, 4
2737 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2738 end do
2739 end do
2740 do i=1, 27
2741 b(:,i,v2,ke_xy) = bb(:,1,i)
2742 g(:,i,1,v2,ke_xy) = bb(:,2,i)
2743 g(:,i,2,v2,ke_xy) = bb(:,3,i)
2744 g(:,i,3,v2,ke_xy) = bb(:,4,i)
2745 end do
2746 end if
2747 end do
2748 end do
2749 return
2750 end subroutine solve_nnode9_var3
2751!OCL SERIAL
2752 subroutine solve_nnode10_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
2753 implicit none
2754 integer, intent(in) :: nnode_v
2755 integer, intent(in) :: jm
2756 integer, intent(in) :: ne2d
2757 real(rp), intent(in) :: d(10,10,10,jm,ne2d)
2758 real(rp), intent(inout) :: b(10,10,2,jm,ne2d)
2759 real(rp), intent(inout) :: g(10,10,jm,ne2d)
2760 logical, intent(in) :: is_top
2761
2762 integer :: ke_xy
2763 integer :: v2
2764 integer :: k, i, j, v, r
2765 integer :: ip
2766 integer :: ipiv(10,10)
2767 real(rp) :: tmpv(10)
2768 real(rp) :: tmpv2(10)
2769 real(rp) :: a(10,10,10)
2770 real(rp) :: bb(10,3,10)
2771 real(rp) :: tmpb(10,3)
2772 real(rp) :: tmp
2773 !------------------------------
2774
2775 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
2776 do ke_xy=1, ne2d
2777 do v2=1, jm
2778 a(:,:,:) = d(:,:,:,v2,ke_xy)
2779 do k=1, 10
2780 tmpv(:) = abs(a(:,k,k))
2781 ipiv(:,k) = k
2782 do i=k+1, 10
2783 do v=1, 10
2784 tmp = abs(a(v,i,k))
2785 if ( tmp > tmpv(v) ) then
2786 ipiv(v,k) = i
2787 tmpv(v) = tmp
2788 end if
2789 end do
2790 end do
2791 do v=1, 10
2792 ip = ipiv(v,k)
2793 if ( ip /= k ) then
2794 do j=1, 10
2795 tmp = a(v,k,j)
2796 a(v,k,j) = a(v,ip,j)
2797 a(v,ip,j) = tmp
2798 end do
2799 end if
2800 end do
2801 !--
2802 tmpv(:) = 1.0_rp / a(:,k,k)
2803 a(:,k,k) = tmpv(:)
2804 do i=k+1, 10
2805 a(:,i,k) = a(:,i,k) * tmpv(:)
2806 end do
2807 do j=k+1, 10
2808 do i=k+1, 10
2809 do v=1, 10
2810 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2811 end do
2812 end do
2813 end do
2814 end do
2815 !---------------------------
2816 if ( is_top ) then
2817 do i=1, 10-1
2818 do v=1, 10
2819 ip = ipiv(v,i)
2820 if (ip /= i) then
2821 tmp = b(v,i,1,v2,ke_xy)
2822 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2823 b(v,ip,1,v2,ke_xy) = tmp
2824 tmp = b(v,i,2,v2,ke_xy)
2825 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2826 b(v,ip,2,v2,ke_xy) = tmp
2827 end if
2828 end do
2829 end do
2830 do i=1, 10
2831 tmpv(:) = b(:,i,1,v2,ke_xy)
2832 tmpv2(:) = b(:,i,2,v2,ke_xy)
2833 do j=1, i-1
2834 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
2835 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
2836 end do
2837 b(:,i,1,v2,ke_xy) = tmpv(:)
2838 b(:,i,2,v2,ke_xy) = tmpv2(:)
2839 end do
2840 do i=10, 1, -1
2841 tmpv(:) = b(:,i,1,v2,ke_xy)
2842 tmpv2(:) = b(:,i,2,v2,ke_xy)
2843 do j=i+1, 10
2844 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
2845 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
2846 end do
2847 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
2848 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
2849 end do
2850 else
2851 do i=1, 10-1
2852 do v=1, 10
2853 ip = ipiv(v,i)
2854 if (ip /= i) then
2855 tmp = b(v,i,1,v2,ke_xy)
2856 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
2857 b(v,ip,1,v2,ke_xy) = tmp
2858 tmp = b(v,i,2,v2,ke_xy)
2859 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
2860 b(v,ip,2,v2,ke_xy) = tmp
2861
2862 tmp = g(v,i,v2,ke_xy)
2863 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
2864 g(v,ip,v2,ke_xy) = tmp
2865 end if
2866 end do
2867 end do
2868 do i=1, 10
2869 do v=1,10
2870 bb(v,1,i) = b(v,i,1,v2,ke_xy)
2871 bb(v,2,i) = b(v,i,2,v2,ke_xy)
2872 bb(v,3,i) = g(v,i,v2,ke_xy)
2873 end do
2874 end do
2875 do i=1, 10
2876 tmpb(:,:) = bb(:,:,i)
2877 do j=1, i-1
2878 do r=1, 3
2879 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
2880 end do
2881 end do
2882 bb(:,:,i) = tmpb(:,:)
2883 end do
2884 do i=10, 1, -1
2885 tmpb(:,:) = bb(:,:,i)
2886 do j=i+1, 10
2887 do r=1, 3
2888 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
2889 end do
2890 end do
2891 do r=1, 3
2892 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
2893 end do
2894 end do
2895 do i=1, 10
2896 b(:,i,1,v2,ke_xy) = bb(:,1,i)
2897 b(:,i,2,v2,ke_xy) = bb(:,2,i)
2898 g(:,i,v2,ke_xy) = bb(:,3,i)
2899 end do
2900 end if
2901 end do
2902 end do
2903 return
2904 end subroutine solve_nnode10_uv
2905!OCL SERIAL
2906 subroutine solve_nnode10_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
2907 implicit none
2908 integer, intent(in) :: nnode_v
2909 integer, intent(in) :: jm
2910 integer, intent(in) :: ne2d
2911 real(rp), intent(in) :: d(10,30,30,jm,ne2d)
2912 real(rp), intent(inout) :: b(10,30,jm,ne2d)
2913 real(rp), intent(inout) :: g(10,30,3,jm,ne2d)
2914 logical, intent(in) :: is_top
2915
2916 integer :: ke_xy
2917 integer :: v2
2918 integer :: k, i, j, v, r
2919 integer :: ip
2920 integer :: ipiv(10,30)
2921 real(rp) :: tmpv(10)
2922 real(rp) :: a(10,30,30)
2923 real(rp) :: bb(10,4,30)
2924 real(rp) :: tmpb(10,4)
2925 real(rp) :: tmp
2926 !------------------------------
2927
2928 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
2929 do ke_xy=1, ne2d
2930 do v2=1, jm
2931 a(:,:,:) = d(:,:,:,v2,ke_xy)
2932 do k=1, 30
2933 tmpv(:) = abs(a(:,k,k))
2934 ipiv(:,k) = k
2935 do i=k+1, 30
2936 do v=1, 10
2937 tmp = abs(a(v,i,k))
2938 if ( tmp > tmpv(v) ) then
2939 ipiv(v,k) = i
2940 tmpv(v) = tmp
2941 end if
2942 end do
2943 end do
2944 do v=1, 10
2945 ip = ipiv(v,k)
2946 if ( ip /= k ) then
2947 do j=1, 30
2948 tmp = a(v,k,j)
2949 a(v,k,j) = a(v,ip,j)
2950 a(v,ip,j) = tmp
2951 end do
2952 end if
2953 end do
2954 !--
2955 tmpv(:) = 1.0_rp / a(:,k,k)
2956 a(:,k,k) = tmpv(:)
2957 do i=k+1, 30
2958 do v=1, 10
2959 a(v,i,k) = a(v,i,k) * tmpv(v)
2960 end do
2961 end do
2962 do j=k+1, 30
2963 do i=k+1, 30
2964 do v=1, 10
2965 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
2966 end do
2967 end do
2968 end do
2969 end do
2970 !---------------------------
2971 if ( is_top ) then
2972 do i=1, 30-1
2973 do v=1, 10
2974 ip = ipiv(v,i)
2975 if (ip /= i) then
2976 tmp = b(v,i,v2,ke_xy)
2977 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
2978 b(v,ip,v2,ke_xy) = tmp
2979 end if
2980 end do
2981 end do
2982 do i=1, 30
2983 tmpv(:) = b(:,i,v2,ke_xy)
2984 do j=1, i-1
2985 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
2986 end do
2987 b(:,i,v2,ke_xy) = tmpv(:)
2988 end do
2989 do i=30, 1, -1
2990 tmpv(:) = b(:,i,v2,ke_xy)
2991 do j=i+1, 30
2992 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
2993 end do
2994 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
2995 end do
2996 else
2997 do i=1, 30-1
2998 do v=1, 10
2999 ip = ipiv(v,i)
3000 if (ip /= i) then
3001 tmp = b(v,i,v2,ke_xy)
3002 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3003 b(v,ip,v2,ke_xy) = tmp
3004
3005 tmp = g(v,i,1,v2,ke_xy)
3006 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
3007 g(v,ip,1,v2,ke_xy) = tmp
3008 tmp = g(v,i,2,v2,ke_xy)
3009 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
3010 g(v,ip,2,v2,ke_xy) = tmp
3011 tmp = g(v,i,3,v2,ke_xy)
3012 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
3013 g(v,ip,3,v2,ke_xy) = tmp
3014 end if
3015 end do
3016 end do
3017 do i=1, 30
3018 do v=1,10
3019 bb(v,1,i) = b(v,i,v2,ke_xy)
3020 bb(v,2,i) = g(v,i,1,v2,ke_xy)
3021 bb(v,3,i) = g(v,i,2,v2,ke_xy)
3022 bb(v,4,i) = g(v,i,3,v2,ke_xy)
3023 end do
3024 end do
3025 do i=1, 30
3026 tmpb(:,:) = bb(:,:,i)
3027 do j=1, i-1
3028 do r=1, 4
3029 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3030 end do
3031 end do
3032 bb(:,:,i) = tmpb(:,:)
3033 end do
3034 do i=30, 1, -1
3035 tmpb(:,:) = bb(:,:,i)
3036 do j=i+1, 30
3037 do r=1, 4
3038 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3039 end do
3040 end do
3041 do r=1, 4
3042 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3043 end do
3044 end do
3045 do i=1, 30
3046 b(:,i,v2,ke_xy) = bb(:,1,i)
3047 g(:,i,1,v2,ke_xy) = bb(:,2,i)
3048 g(:,i,2,v2,ke_xy) = bb(:,3,i)
3049 g(:,i,3,v2,ke_xy) = bb(:,4,i)
3050 end do
3051 end if
3052 end do
3053 end do
3054 return
3055 end subroutine solve_nnode10_var3
3056!OCL SERIAL
3057 subroutine solve_nnode11_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
3058 implicit none
3059 integer, intent(in) :: nnode_v
3060 integer, intent(in) :: jm
3061 integer, intent(in) :: ne2d
3062 real(rp), intent(in) :: d(22,11,11,jm,ne2d)
3063 real(rp), intent(inout) :: b(22,11,2,jm,ne2d)
3064 real(rp), intent(inout) :: g(22,11,jm,ne2d)
3065 logical, intent(in) :: is_top
3066
3067 integer :: ke_xy
3068 integer :: v2
3069 integer :: k, i, j, v, r
3070 integer :: ip
3071 integer :: ipiv(22,11)
3072 real(rp) :: tmpv(22)
3073 real(rp) :: tmpv2(22)
3074 real(rp) :: a(22,11,11)
3075 real(rp) :: bb(22,3,11)
3076 real(rp) :: tmpb(22,3)
3077 real(rp) :: tmp
3078 !------------------------------
3079
3080 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
3081 do ke_xy=1, ne2d
3082 do v2=1, jm
3083 a(:,:,:) = d(:,:,:,v2,ke_xy)
3084 do k=1, 11
3085 tmpv(:) = abs(a(:,k,k))
3086 ipiv(:,k) = k
3087 do i=k+1, 11
3088 do v=1, 22
3089 tmp = abs(a(v,i,k))
3090 if ( tmp > tmpv(v) ) then
3091 ipiv(v,k) = i
3092 tmpv(v) = tmp
3093 end if
3094 end do
3095 end do
3096 do v=1, 22
3097 ip = ipiv(v,k)
3098 if ( ip /= k ) then
3099 do j=1, 11
3100 tmp = a(v,k,j)
3101 a(v,k,j) = a(v,ip,j)
3102 a(v,ip,j) = tmp
3103 end do
3104 end if
3105 end do
3106 !--
3107 tmpv(:) = 1.0_rp / a(:,k,k)
3108 a(:,k,k) = tmpv(:)
3109 do i=k+1, 11
3110 a(:,i,k) = a(:,i,k) * tmpv(:)
3111 end do
3112 do j=k+1, 11
3113 do i=k+1, 11
3114 do v=1, 22
3115 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3116 end do
3117 end do
3118 end do
3119 end do
3120 !---------------------------
3121 if ( is_top ) then
3122 do i=1, 11-1
3123 do v=1, 22
3124 ip = ipiv(v,i)
3125 if (ip /= i) then
3126 tmp = b(v,i,1,v2,ke_xy)
3127 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3128 b(v,ip,1,v2,ke_xy) = tmp
3129 tmp = b(v,i,2,v2,ke_xy)
3130 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3131 b(v,ip,2,v2,ke_xy) = tmp
3132 end if
3133 end do
3134 end do
3135 do i=1, 11
3136 tmpv(:) = b(:,i,1,v2,ke_xy)
3137 tmpv2(:) = b(:,i,2,v2,ke_xy)
3138 do j=1, i-1
3139 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
3140 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
3141 end do
3142 b(:,i,1,v2,ke_xy) = tmpv(:)
3143 b(:,i,2,v2,ke_xy) = tmpv2(:)
3144 end do
3145 do i=11, 1, -1
3146 tmpv(:) = b(:,i,1,v2,ke_xy)
3147 tmpv2(:) = b(:,i,2,v2,ke_xy)
3148 do j=i+1, 11
3149 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
3150 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
3151 end do
3152 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
3153 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
3154 end do
3155 else
3156 do i=1, 11-1
3157 do v=1, 22
3158 ip = ipiv(v,i)
3159 if (ip /= i) then
3160 tmp = b(v,i,1,v2,ke_xy)
3161 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3162 b(v,ip,1,v2,ke_xy) = tmp
3163 tmp = b(v,i,2,v2,ke_xy)
3164 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3165 b(v,ip,2,v2,ke_xy) = tmp
3166
3167 tmp = g(v,i,v2,ke_xy)
3168 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
3169 g(v,ip,v2,ke_xy) = tmp
3170 end if
3171 end do
3172 end do
3173 do i=1, 11
3174 do v=1,22
3175 bb(v,1,i) = b(v,i,1,v2,ke_xy)
3176 bb(v,2,i) = b(v,i,2,v2,ke_xy)
3177 bb(v,3,i) = g(v,i,v2,ke_xy)
3178 end do
3179 end do
3180 do i=1, 11
3181 tmpb(:,:) = bb(:,:,i)
3182 do j=1, i-1
3183 do r=1, 3
3184 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3185 end do
3186 end do
3187 bb(:,:,i) = tmpb(:,:)
3188 end do
3189 do i=11, 1, -1
3190 tmpb(:,:) = bb(:,:,i)
3191 do j=i+1, 11
3192 do r=1, 3
3193 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3194 end do
3195 end do
3196 do r=1, 3
3197 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3198 end do
3199 end do
3200 do i=1, 11
3201 b(:,i,1,v2,ke_xy) = bb(:,1,i)
3202 b(:,i,2,v2,ke_xy) = bb(:,2,i)
3203 g(:,i,v2,ke_xy) = bb(:,3,i)
3204 end do
3205 end if
3206 end do
3207 end do
3208 return
3209 end subroutine solve_nnode11_uv
3210!OCL SERIAL
3211 subroutine solve_nnode11_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
3212 implicit none
3213 integer, intent(in) :: nnode_v
3214 integer, intent(in) :: jm
3215 integer, intent(in) :: ne2d
3216 real(rp), intent(in) :: d(22,33,33,jm,ne2d)
3217 real(rp), intent(inout) :: b(22,33,jm,ne2d)
3218 real(rp), intent(inout) :: g(22,33,3,jm,ne2d)
3219 logical, intent(in) :: is_top
3220
3221 integer :: ke_xy
3222 integer :: v2
3223 integer :: k, i, j, v, r
3224 integer :: ip
3225 integer :: ipiv(22,33)
3226 real(rp) :: tmpv(22)
3227 real(rp) :: a(22,33,33)
3228 real(rp) :: bb(22,4,33)
3229 real(rp) :: tmpb(22,4)
3230 real(rp) :: tmp
3231 !------------------------------
3232
3233 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
3234 do ke_xy=1, ne2d
3235 do v2=1, jm
3236 a(:,:,:) = d(:,:,:,v2,ke_xy)
3237 do k=1, 33
3238 tmpv(:) = abs(a(:,k,k))
3239 ipiv(:,k) = k
3240 do i=k+1, 33
3241 do v=1, 22
3242 tmp = abs(a(v,i,k))
3243 if ( tmp > tmpv(v) ) then
3244 ipiv(v,k) = i
3245 tmpv(v) = tmp
3246 end if
3247 end do
3248 end do
3249 do v=1, 22
3250 ip = ipiv(v,k)
3251 if ( ip /= k ) then
3252 do j=1, 33
3253 tmp = a(v,k,j)
3254 a(v,k,j) = a(v,ip,j)
3255 a(v,ip,j) = tmp
3256 end do
3257 end if
3258 end do
3259 !--
3260 tmpv(:) = 1.0_rp / a(:,k,k)
3261 a(:,k,k) = tmpv(:)
3262 do i=k+1, 33
3263 do v=1, 22
3264 a(v,i,k) = a(v,i,k) * tmpv(v)
3265 end do
3266 end do
3267 do j=k+1, 33
3268 do i=k+1, 33
3269 do v=1, 22
3270 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3271 end do
3272 end do
3273 end do
3274 end do
3275 !---------------------------
3276 if ( is_top ) then
3277 do i=1, 33-1
3278 do v=1, 22
3279 ip = ipiv(v,i)
3280 if (ip /= i) then
3281 tmp = b(v,i,v2,ke_xy)
3282 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3283 b(v,ip,v2,ke_xy) = tmp
3284 end if
3285 end do
3286 end do
3287 do i=1, 33
3288 tmpv(:) = b(:,i,v2,ke_xy)
3289 do j=1, i-1
3290 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
3291 end do
3292 b(:,i,v2,ke_xy) = tmpv(:)
3293 end do
3294 do i=33, 1, -1
3295 tmpv(:) = b(:,i,v2,ke_xy)
3296 do j=i+1, 33
3297 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
3298 end do
3299 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
3300 end do
3301 else
3302 do i=1, 33-1
3303 do v=1, 22
3304 ip = ipiv(v,i)
3305 if (ip /= i) then
3306 tmp = b(v,i,v2,ke_xy)
3307 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3308 b(v,ip,v2,ke_xy) = tmp
3309
3310 tmp = g(v,i,1,v2,ke_xy)
3311 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
3312 g(v,ip,1,v2,ke_xy) = tmp
3313 tmp = g(v,i,2,v2,ke_xy)
3314 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
3315 g(v,ip,2,v2,ke_xy) = tmp
3316 tmp = g(v,i,3,v2,ke_xy)
3317 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
3318 g(v,ip,3,v2,ke_xy) = tmp
3319 end if
3320 end do
3321 end do
3322 do i=1, 33
3323 do v=1,22
3324 bb(v,1,i) = b(v,i,v2,ke_xy)
3325 bb(v,2,i) = g(v,i,1,v2,ke_xy)
3326 bb(v,3,i) = g(v,i,2,v2,ke_xy)
3327 bb(v,4,i) = g(v,i,3,v2,ke_xy)
3328 end do
3329 end do
3330 do i=1, 33
3331 tmpb(:,:) = bb(:,:,i)
3332 do j=1, i-1
3333 do r=1, 4
3334 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3335 end do
3336 end do
3337 bb(:,:,i) = tmpb(:,:)
3338 end do
3339 do i=33, 1, -1
3340 tmpb(:,:) = bb(:,:,i)
3341 do j=i+1, 33
3342 do r=1, 4
3343 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3344 end do
3345 end do
3346 do r=1, 4
3347 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3348 end do
3349 end do
3350 do i=1, 33
3351 b(:,i,v2,ke_xy) = bb(:,1,i)
3352 g(:,i,1,v2,ke_xy) = bb(:,2,i)
3353 g(:,i,2,v2,ke_xy) = bb(:,3,i)
3354 g(:,i,3,v2,ke_xy) = bb(:,4,i)
3355 end do
3356 end if
3357 end do
3358 end do
3359 return
3360 end subroutine solve_nnode11_var3
3361!OCL SERIAL
3362 subroutine solve_nnode12_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
3363 implicit none
3364 integer, intent(in) :: nnode_v
3365 integer, intent(in) :: jm
3366 integer, intent(in) :: ne2d
3367 real(rp), intent(in) :: d(24,12,12,jm,ne2d)
3368 real(rp), intent(inout) :: b(24,12,2,jm,ne2d)
3369 real(rp), intent(inout) :: g(24,12,jm,ne2d)
3370 logical, intent(in) :: is_top
3371
3372 integer :: ke_xy
3373 integer :: v2
3374 integer :: k, i, j, v, r
3375 integer :: ip
3376 integer :: ipiv(24,12)
3377 real(rp) :: tmpv(24)
3378 real(rp) :: tmpv2(24)
3379 real(rp) :: a(24,12,12)
3380 real(rp) :: bb(24,3,12)
3381 real(rp) :: tmpb(24,3)
3382 real(rp) :: tmp
3383 !------------------------------
3384
3385 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
3386 do ke_xy=1, ne2d
3387 do v2=1, jm
3388 a(:,:,:) = d(:,:,:,v2,ke_xy)
3389 do k=1, 12
3390 tmpv(:) = abs(a(:,k,k))
3391 ipiv(:,k) = k
3392 do i=k+1, 12
3393 do v=1, 24
3394 tmp = abs(a(v,i,k))
3395 if ( tmp > tmpv(v) ) then
3396 ipiv(v,k) = i
3397 tmpv(v) = tmp
3398 end if
3399 end do
3400 end do
3401 do v=1, 24
3402 ip = ipiv(v,k)
3403 if ( ip /= k ) then
3404 do j=1, 12
3405 tmp = a(v,k,j)
3406 a(v,k,j) = a(v,ip,j)
3407 a(v,ip,j) = tmp
3408 end do
3409 end if
3410 end do
3411 !--
3412 tmpv(:) = 1.0_rp / a(:,k,k)
3413 a(:,k,k) = tmpv(:)
3414 do i=k+1, 12
3415 a(:,i,k) = a(:,i,k) * tmpv(:)
3416 end do
3417 do j=k+1, 12
3418 do i=k+1, 12
3419 do v=1, 24
3420 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3421 end do
3422 end do
3423 end do
3424 end do
3425 !---------------------------
3426 if ( is_top ) then
3427 do i=1, 12-1
3428 do v=1, 24
3429 ip = ipiv(v,i)
3430 if (ip /= i) then
3431 tmp = b(v,i,1,v2,ke_xy)
3432 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3433 b(v,ip,1,v2,ke_xy) = tmp
3434 tmp = b(v,i,2,v2,ke_xy)
3435 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3436 b(v,ip,2,v2,ke_xy) = tmp
3437 end if
3438 end do
3439 end do
3440 do i=1, 12
3441 tmpv(:) = b(:,i,1,v2,ke_xy)
3442 tmpv2(:) = b(:,i,2,v2,ke_xy)
3443 do j=1, i-1
3444 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
3445 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
3446 end do
3447 b(:,i,1,v2,ke_xy) = tmpv(:)
3448 b(:,i,2,v2,ke_xy) = tmpv2(:)
3449 end do
3450 do i=12, 1, -1
3451 tmpv(:) = b(:,i,1,v2,ke_xy)
3452 tmpv2(:) = b(:,i,2,v2,ke_xy)
3453 do j=i+1, 12
3454 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
3455 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
3456 end do
3457 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
3458 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
3459 end do
3460 else
3461 do i=1, 12-1
3462 do v=1, 24
3463 ip = ipiv(v,i)
3464 if (ip /= i) then
3465 tmp = b(v,i,1,v2,ke_xy)
3466 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3467 b(v,ip,1,v2,ke_xy) = tmp
3468 tmp = b(v,i,2,v2,ke_xy)
3469 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3470 b(v,ip,2,v2,ke_xy) = tmp
3471
3472 tmp = g(v,i,v2,ke_xy)
3473 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
3474 g(v,ip,v2,ke_xy) = tmp
3475 end if
3476 end do
3477 end do
3478 do i=1, 12
3479 do v=1,24
3480 bb(v,1,i) = b(v,i,1,v2,ke_xy)
3481 bb(v,2,i) = b(v,i,2,v2,ke_xy)
3482 bb(v,3,i) = g(v,i,v2,ke_xy)
3483 end do
3484 end do
3485 do i=1, 12
3486 tmpb(:,:) = bb(:,:,i)
3487 do j=1, i-1
3488 do r=1, 3
3489 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3490 end do
3491 end do
3492 bb(:,:,i) = tmpb(:,:)
3493 end do
3494 do i=12, 1, -1
3495 tmpb(:,:) = bb(:,:,i)
3496 do j=i+1, 12
3497 do r=1, 3
3498 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3499 end do
3500 end do
3501 do r=1, 3
3502 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3503 end do
3504 end do
3505 do i=1, 12
3506 b(:,i,1,v2,ke_xy) = bb(:,1,i)
3507 b(:,i,2,v2,ke_xy) = bb(:,2,i)
3508 g(:,i,v2,ke_xy) = bb(:,3,i)
3509 end do
3510 end if
3511 end do
3512 end do
3513 return
3514 end subroutine solve_nnode12_uv
3515!OCL SERIAL
3516 subroutine solve_nnode12_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
3517 implicit none
3518 integer, intent(in) :: nnode_v
3519 integer, intent(in) :: jm
3520 integer, intent(in) :: ne2d
3521 real(rp), intent(in) :: d(24,36,36,jm,ne2d)
3522 real(rp), intent(inout) :: b(24,36,jm,ne2d)
3523 real(rp), intent(inout) :: g(24,36,3,jm,ne2d)
3524 logical, intent(in) :: is_top
3525
3526 integer :: ke_xy
3527 integer :: v2
3528 integer :: k, i, j, v, r
3529 integer :: ip
3530 integer :: ipiv(24,36)
3531 real(rp) :: tmpv(24)
3532 real(rp) :: a(24,36,36)
3533 real(rp) :: bb(24,4,36)
3534 real(rp) :: tmpb(24,4)
3535 real(rp) :: tmp
3536 !------------------------------
3537
3538 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
3539 do ke_xy=1, ne2d
3540 do v2=1, jm
3541 a(:,:,:) = d(:,:,:,v2,ke_xy)
3542 do k=1, 36
3543 tmpv(:) = abs(a(:,k,k))
3544 ipiv(:,k) = k
3545 do i=k+1, 36
3546 do v=1, 24
3547 tmp = abs(a(v,i,k))
3548 if ( tmp > tmpv(v) ) then
3549 ipiv(v,k) = i
3550 tmpv(v) = tmp
3551 end if
3552 end do
3553 end do
3554 do v=1, 24
3555 ip = ipiv(v,k)
3556 if ( ip /= k ) then
3557 do j=1, 36
3558 tmp = a(v,k,j)
3559 a(v,k,j) = a(v,ip,j)
3560 a(v,ip,j) = tmp
3561 end do
3562 end if
3563 end do
3564 !--
3565 tmpv(:) = 1.0_rp / a(:,k,k)
3566 a(:,k,k) = tmpv(:)
3567 do i=k+1, 36
3568 do v=1, 24
3569 a(v,i,k) = a(v,i,k) * tmpv(v)
3570 end do
3571 end do
3572 do j=k+1, 36
3573 do i=k+1, 36
3574 do v=1, 24
3575 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3576 end do
3577 end do
3578 end do
3579 end do
3580 !---------------------------
3581 if ( is_top ) then
3582 do i=1, 36-1
3583 do v=1, 24
3584 ip = ipiv(v,i)
3585 if (ip /= i) then
3586 tmp = b(v,i,v2,ke_xy)
3587 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3588 b(v,ip,v2,ke_xy) = tmp
3589 end if
3590 end do
3591 end do
3592 do i=1, 36
3593 tmpv(:) = b(:,i,v2,ke_xy)
3594 do j=1, i-1
3595 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
3596 end do
3597 b(:,i,v2,ke_xy) = tmpv(:)
3598 end do
3599 do i=36, 1, -1
3600 tmpv(:) = b(:,i,v2,ke_xy)
3601 do j=i+1, 36
3602 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
3603 end do
3604 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
3605 end do
3606 else
3607 do i=1, 36-1
3608 do v=1, 24
3609 ip = ipiv(v,i)
3610 if (ip /= i) then
3611 tmp = b(v,i,v2,ke_xy)
3612 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3613 b(v,ip,v2,ke_xy) = tmp
3614
3615 tmp = g(v,i,1,v2,ke_xy)
3616 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
3617 g(v,ip,1,v2,ke_xy) = tmp
3618 tmp = g(v,i,2,v2,ke_xy)
3619 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
3620 g(v,ip,2,v2,ke_xy) = tmp
3621 tmp = g(v,i,3,v2,ke_xy)
3622 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
3623 g(v,ip,3,v2,ke_xy) = tmp
3624 end if
3625 end do
3626 end do
3627 do i=1, 36
3628 do v=1,24
3629 bb(v,1,i) = b(v,i,v2,ke_xy)
3630 bb(v,2,i) = g(v,i,1,v2,ke_xy)
3631 bb(v,3,i) = g(v,i,2,v2,ke_xy)
3632 bb(v,4,i) = g(v,i,3,v2,ke_xy)
3633 end do
3634 end do
3635 do i=1, 36
3636 tmpb(:,:) = bb(:,:,i)
3637 do j=1, i-1
3638 do r=1, 4
3639 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3640 end do
3641 end do
3642 bb(:,:,i) = tmpb(:,:)
3643 end do
3644 do i=36, 1, -1
3645 tmpb(:,:) = bb(:,:,i)
3646 do j=i+1, 36
3647 do r=1, 4
3648 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3649 end do
3650 end do
3651 do r=1, 4
3652 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3653 end do
3654 end do
3655 do i=1, 36
3656 b(:,i,v2,ke_xy) = bb(:,1,i)
3657 g(:,i,1,v2,ke_xy) = bb(:,2,i)
3658 g(:,i,2,v2,ke_xy) = bb(:,3,i)
3659 g(:,i,3,v2,ke_xy) = bb(:,4,i)
3660 end do
3661 end if
3662 end do
3663 end do
3664 return
3665 end subroutine solve_nnode12_var3
3666!OCL SERIAL
3667 subroutine solve_nnode13_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
3668 implicit none
3669 integer, intent(in) :: nnode_v
3670 integer, intent(in) :: jm
3671 integer, intent(in) :: ne2d
3672 real(rp), intent(in) :: d(26,13,13,jm,ne2d)
3673 real(rp), intent(inout) :: b(26,13,2,jm,ne2d)
3674 real(rp), intent(inout) :: g(26,13,jm,ne2d)
3675 logical, intent(in) :: is_top
3676
3677 integer :: ke_xy
3678 integer :: v2
3679 integer :: k, i, j, v, r
3680 integer :: ip
3681 integer :: ipiv(26,13)
3682 real(rp) :: tmpv(26)
3683 real(rp) :: tmpv2(26)
3684 real(rp) :: a(26,13,13)
3685 real(rp) :: bb(26,3,13)
3686 real(rp) :: tmpb(26,3)
3687 real(rp) :: tmp
3688 !------------------------------
3689
3690 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
3691 do ke_xy=1, ne2d
3692 do v2=1, jm
3693 a(:,:,:) = d(:,:,:,v2,ke_xy)
3694 do k=1, 13
3695 tmpv(:) = abs(a(:,k,k))
3696 ipiv(:,k) = k
3697 do i=k+1, 13
3698 do v=1, 26
3699 tmp = abs(a(v,i,k))
3700 if ( tmp > tmpv(v) ) then
3701 ipiv(v,k) = i
3702 tmpv(v) = tmp
3703 end if
3704 end do
3705 end do
3706 do v=1, 26
3707 ip = ipiv(v,k)
3708 if ( ip /= k ) then
3709 do j=1, 13
3710 tmp = a(v,k,j)
3711 a(v,k,j) = a(v,ip,j)
3712 a(v,ip,j) = tmp
3713 end do
3714 end if
3715 end do
3716 !--
3717 tmpv(:) = 1.0_rp / a(:,k,k)
3718 a(:,k,k) = tmpv(:)
3719 do i=k+1, 13
3720 a(:,i,k) = a(:,i,k) * tmpv(:)
3721 end do
3722 do j=k+1, 13
3723 do i=k+1, 13
3724 do v=1, 26
3725 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3726 end do
3727 end do
3728 end do
3729 end do
3730 !---------------------------
3731 if ( is_top ) then
3732 do i=1, 13-1
3733 do v=1, 26
3734 ip = ipiv(v,i)
3735 if (ip /= i) then
3736 tmp = b(v,i,1,v2,ke_xy)
3737 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3738 b(v,ip,1,v2,ke_xy) = tmp
3739 tmp = b(v,i,2,v2,ke_xy)
3740 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3741 b(v,ip,2,v2,ke_xy) = tmp
3742 end if
3743 end do
3744 end do
3745 do i=1, 13
3746 tmpv(:) = b(:,i,1,v2,ke_xy)
3747 tmpv2(:) = b(:,i,2,v2,ke_xy)
3748 do j=1, i-1
3749 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
3750 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
3751 end do
3752 b(:,i,1,v2,ke_xy) = tmpv(:)
3753 b(:,i,2,v2,ke_xy) = tmpv2(:)
3754 end do
3755 do i=13, 1, -1
3756 tmpv(:) = b(:,i,1,v2,ke_xy)
3757 tmpv2(:) = b(:,i,2,v2,ke_xy)
3758 do j=i+1, 13
3759 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
3760 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
3761 end do
3762 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
3763 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
3764 end do
3765 else
3766 do i=1, 13-1
3767 do v=1, 26
3768 ip = ipiv(v,i)
3769 if (ip /= i) then
3770 tmp = b(v,i,1,v2,ke_xy)
3771 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
3772 b(v,ip,1,v2,ke_xy) = tmp
3773 tmp = b(v,i,2,v2,ke_xy)
3774 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
3775 b(v,ip,2,v2,ke_xy) = tmp
3776
3777 tmp = g(v,i,v2,ke_xy)
3778 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
3779 g(v,ip,v2,ke_xy) = tmp
3780 end if
3781 end do
3782 end do
3783 do i=1, 13
3784 do v=1,26
3785 bb(v,1,i) = b(v,i,1,v2,ke_xy)
3786 bb(v,2,i) = b(v,i,2,v2,ke_xy)
3787 bb(v,3,i) = g(v,i,v2,ke_xy)
3788 end do
3789 end do
3790 do i=1, 13
3791 tmpb(:,:) = bb(:,:,i)
3792 do j=1, i-1
3793 do r=1, 3
3794 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3795 end do
3796 end do
3797 bb(:,:,i) = tmpb(:,:)
3798 end do
3799 do i=13, 1, -1
3800 tmpb(:,:) = bb(:,:,i)
3801 do j=i+1, 13
3802 do r=1, 3
3803 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3804 end do
3805 end do
3806 do r=1, 3
3807 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3808 end do
3809 end do
3810 do i=1, 13
3811 b(:,i,1,v2,ke_xy) = bb(:,1,i)
3812 b(:,i,2,v2,ke_xy) = bb(:,2,i)
3813 g(:,i,v2,ke_xy) = bb(:,3,i)
3814 end do
3815 end if
3816 end do
3817 end do
3818 return
3819 end subroutine solve_nnode13_uv
3820!OCL SERIAL
3821 subroutine solve_nnode13_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
3822 implicit none
3823 integer, intent(in) :: nnode_v
3824 integer, intent(in) :: jm
3825 integer, intent(in) :: ne2d
3826 real(rp), intent(in) :: d(26,39,39,jm,ne2d)
3827 real(rp), intent(inout) :: b(26,39,jm,ne2d)
3828 real(rp), intent(inout) :: g(26,39,3,jm,ne2d)
3829 logical, intent(in) :: is_top
3830
3831 integer :: ke_xy
3832 integer :: v2
3833 integer :: k, i, j, v, r
3834 integer :: ip
3835 integer :: ipiv(26,39)
3836 real(rp) :: tmpv(26)
3837 real(rp) :: a(26,39,39)
3838 real(rp) :: bb(26,4,39)
3839 real(rp) :: tmpb(26,4)
3840 real(rp) :: tmp
3841 !------------------------------
3842
3843 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
3844 do ke_xy=1, ne2d
3845 do v2=1, jm
3846 a(:,:,:) = d(:,:,:,v2,ke_xy)
3847 do k=1, 39
3848 tmpv(:) = abs(a(:,k,k))
3849 ipiv(:,k) = k
3850 do i=k+1, 39
3851 do v=1, 26
3852 tmp = abs(a(v,i,k))
3853 if ( tmp > tmpv(v) ) then
3854 ipiv(v,k) = i
3855 tmpv(v) = tmp
3856 end if
3857 end do
3858 end do
3859 do v=1, 26
3860 ip = ipiv(v,k)
3861 if ( ip /= k ) then
3862 do j=1, 39
3863 tmp = a(v,k,j)
3864 a(v,k,j) = a(v,ip,j)
3865 a(v,ip,j) = tmp
3866 end do
3867 end if
3868 end do
3869 !--
3870 tmpv(:) = 1.0_rp / a(:,k,k)
3871 a(:,k,k) = tmpv(:)
3872 do i=k+1, 39
3873 do v=1, 26
3874 a(v,i,k) = a(v,i,k) * tmpv(v)
3875 end do
3876 end do
3877 do j=k+1, 39
3878 do i=k+1, 39
3879 do v=1, 26
3880 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
3881 end do
3882 end do
3883 end do
3884 end do
3885 !---------------------------
3886 if ( is_top ) then
3887 do i=1, 39-1
3888 do v=1, 26
3889 ip = ipiv(v,i)
3890 if (ip /= i) then
3891 tmp = b(v,i,v2,ke_xy)
3892 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3893 b(v,ip,v2,ke_xy) = tmp
3894 end if
3895 end do
3896 end do
3897 do i=1, 39
3898 tmpv(:) = b(:,i,v2,ke_xy)
3899 do j=1, i-1
3900 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
3901 end do
3902 b(:,i,v2,ke_xy) = tmpv(:)
3903 end do
3904 do i=39, 1, -1
3905 tmpv(:) = b(:,i,v2,ke_xy)
3906 do j=i+1, 39
3907 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
3908 end do
3909 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
3910 end do
3911 else
3912 do i=1, 39-1
3913 do v=1, 26
3914 ip = ipiv(v,i)
3915 if (ip /= i) then
3916 tmp = b(v,i,v2,ke_xy)
3917 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
3918 b(v,ip,v2,ke_xy) = tmp
3919
3920 tmp = g(v,i,1,v2,ke_xy)
3921 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
3922 g(v,ip,1,v2,ke_xy) = tmp
3923 tmp = g(v,i,2,v2,ke_xy)
3924 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
3925 g(v,ip,2,v2,ke_xy) = tmp
3926 tmp = g(v,i,3,v2,ke_xy)
3927 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
3928 g(v,ip,3,v2,ke_xy) = tmp
3929 end if
3930 end do
3931 end do
3932 do i=1, 39
3933 do v=1,26
3934 bb(v,1,i) = b(v,i,v2,ke_xy)
3935 bb(v,2,i) = g(v,i,1,v2,ke_xy)
3936 bb(v,3,i) = g(v,i,2,v2,ke_xy)
3937 bb(v,4,i) = g(v,i,3,v2,ke_xy)
3938 end do
3939 end do
3940 do i=1, 39
3941 tmpb(:,:) = bb(:,:,i)
3942 do j=1, i-1
3943 do r=1, 4
3944 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
3945 end do
3946 end do
3947 bb(:,:,i) = tmpb(:,:)
3948 end do
3949 do i=39, 1, -1
3950 tmpb(:,:) = bb(:,:,i)
3951 do j=i+1, 39
3952 do r=1, 4
3953 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
3954 end do
3955 end do
3956 do r=1, 4
3957 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
3958 end do
3959 end do
3960 do i=1, 39
3961 b(:,i,v2,ke_xy) = bb(:,1,i)
3962 g(:,i,1,v2,ke_xy) = bb(:,2,i)
3963 g(:,i,2,v2,ke_xy) = bb(:,3,i)
3964 g(:,i,3,v2,ke_xy) = bb(:,4,i)
3965 end do
3966 end if
3967 end do
3968 end do
3969 return
3970 end subroutine solve_nnode13_var3
3971!OCL SERIAL
3972 subroutine solve_nnode14_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
3973 implicit none
3974 integer, intent(in) :: nnode_v
3975 integer, intent(in) :: jm
3976 integer, intent(in) :: ne2d
3977 real(rp), intent(in) :: d(28,14,14,jm,ne2d)
3978 real(rp), intent(inout) :: b(28,14,2,jm,ne2d)
3979 real(rp), intent(inout) :: g(28,14,jm,ne2d)
3980 logical, intent(in) :: is_top
3981
3982 integer :: ke_xy
3983 integer :: v2
3984 integer :: k, i, j, v, r
3985 integer :: ip
3986 integer :: ipiv(28,14)
3987 real(rp) :: tmpv(28)
3988 real(rp) :: tmpv2(28)
3989 real(rp) :: a(28,14,14)
3990 real(rp) :: bb(28,3,14)
3991 real(rp) :: tmpb(28,3)
3992 real(rp) :: tmp
3993 !------------------------------
3994
3995 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
3996 do ke_xy=1, ne2d
3997 do v2=1, jm
3998 a(:,:,:) = d(:,:,:,v2,ke_xy)
3999 do k=1, 14
4000 tmpv(:) = abs(a(:,k,k))
4001 ipiv(:,k) = k
4002 do i=k+1, 14
4003 do v=1, 28
4004 tmp = abs(a(v,i,k))
4005 if ( tmp > tmpv(v) ) then
4006 ipiv(v,k) = i
4007 tmpv(v) = tmp
4008 end if
4009 end do
4010 end do
4011 do v=1, 28
4012 ip = ipiv(v,k)
4013 if ( ip /= k ) then
4014 do j=1, 14
4015 tmp = a(v,k,j)
4016 a(v,k,j) = a(v,ip,j)
4017 a(v,ip,j) = tmp
4018 end do
4019 end if
4020 end do
4021 !--
4022 tmpv(:) = 1.0_rp / a(:,k,k)
4023 a(:,k,k) = tmpv(:)
4024 do i=k+1, 14
4025 a(:,i,k) = a(:,i,k) * tmpv(:)
4026 end do
4027 do j=k+1, 14
4028 do i=k+1, 14
4029 do v=1, 28
4030 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4031 end do
4032 end do
4033 end do
4034 end do
4035 !---------------------------
4036 if ( is_top ) then
4037 do i=1, 14-1
4038 do v=1, 28
4039 ip = ipiv(v,i)
4040 if (ip /= i) then
4041 tmp = b(v,i,1,v2,ke_xy)
4042 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4043 b(v,ip,1,v2,ke_xy) = tmp
4044 tmp = b(v,i,2,v2,ke_xy)
4045 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4046 b(v,ip,2,v2,ke_xy) = tmp
4047 end if
4048 end do
4049 end do
4050 do i=1, 14
4051 tmpv(:) = b(:,i,1,v2,ke_xy)
4052 tmpv2(:) = b(:,i,2,v2,ke_xy)
4053 do j=1, i-1
4054 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
4055 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
4056 end do
4057 b(:,i,1,v2,ke_xy) = tmpv(:)
4058 b(:,i,2,v2,ke_xy) = tmpv2(:)
4059 end do
4060 do i=14, 1, -1
4061 tmpv(:) = b(:,i,1,v2,ke_xy)
4062 tmpv2(:) = b(:,i,2,v2,ke_xy)
4063 do j=i+1, 14
4064 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
4065 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
4066 end do
4067 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
4068 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
4069 end do
4070 else
4071 do i=1, 14-1
4072 do v=1, 28
4073 ip = ipiv(v,i)
4074 if (ip /= i) then
4075 tmp = b(v,i,1,v2,ke_xy)
4076 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4077 b(v,ip,1,v2,ke_xy) = tmp
4078 tmp = b(v,i,2,v2,ke_xy)
4079 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4080 b(v,ip,2,v2,ke_xy) = tmp
4081
4082 tmp = g(v,i,v2,ke_xy)
4083 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
4084 g(v,ip,v2,ke_xy) = tmp
4085 end if
4086 end do
4087 end do
4088 do i=1, 14
4089 do v=1,28
4090 bb(v,1,i) = b(v,i,1,v2,ke_xy)
4091 bb(v,2,i) = b(v,i,2,v2,ke_xy)
4092 bb(v,3,i) = g(v,i,v2,ke_xy)
4093 end do
4094 end do
4095 do i=1, 14
4096 tmpb(:,:) = bb(:,:,i)
4097 do j=1, i-1
4098 do r=1, 3
4099 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4100 end do
4101 end do
4102 bb(:,:,i) = tmpb(:,:)
4103 end do
4104 do i=14, 1, -1
4105 tmpb(:,:) = bb(:,:,i)
4106 do j=i+1, 14
4107 do r=1, 3
4108 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4109 end do
4110 end do
4111 do r=1, 3
4112 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4113 end do
4114 end do
4115 do i=1, 14
4116 b(:,i,1,v2,ke_xy) = bb(:,1,i)
4117 b(:,i,2,v2,ke_xy) = bb(:,2,i)
4118 g(:,i,v2,ke_xy) = bb(:,3,i)
4119 end do
4120 end if
4121 end do
4122 end do
4123 return
4124 end subroutine solve_nnode14_uv
4125!OCL SERIAL
4126 subroutine solve_nnode14_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
4127 implicit none
4128 integer, intent(in) :: nnode_v
4129 integer, intent(in) :: jm
4130 integer, intent(in) :: ne2d
4131 real(rp), intent(in) :: d(28,42,42,jm,ne2d)
4132 real(rp), intent(inout) :: b(28,42,jm,ne2d)
4133 real(rp), intent(inout) :: g(28,42,3,jm,ne2d)
4134 logical, intent(in) :: is_top
4135
4136 integer :: ke_xy
4137 integer :: v2
4138 integer :: k, i, j, v, r
4139 integer :: ip
4140 integer :: ipiv(28,42)
4141 real(rp) :: tmpv(28)
4142 real(rp) :: a(28,42,42)
4143 real(rp) :: bb(28,4,42)
4144 real(rp) :: tmpb(28,4)
4145 real(rp) :: tmp
4146 !------------------------------
4147
4148 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
4149 do ke_xy=1, ne2d
4150 do v2=1, jm
4151 a(:,:,:) = d(:,:,:,v2,ke_xy)
4152 do k=1, 42
4153 tmpv(:) = abs(a(:,k,k))
4154 ipiv(:,k) = k
4155 do i=k+1, 42
4156 do v=1, 28
4157 tmp = abs(a(v,i,k))
4158 if ( tmp > tmpv(v) ) then
4159 ipiv(v,k) = i
4160 tmpv(v) = tmp
4161 end if
4162 end do
4163 end do
4164 do v=1, 28
4165 ip = ipiv(v,k)
4166 if ( ip /= k ) then
4167 do j=1, 42
4168 tmp = a(v,k,j)
4169 a(v,k,j) = a(v,ip,j)
4170 a(v,ip,j) = tmp
4171 end do
4172 end if
4173 end do
4174 !--
4175 tmpv(:) = 1.0_rp / a(:,k,k)
4176 a(:,k,k) = tmpv(:)
4177 do i=k+1, 42
4178 do v=1, 28
4179 a(v,i,k) = a(v,i,k) * tmpv(v)
4180 end do
4181 end do
4182 do j=k+1, 42
4183 do i=k+1, 42
4184 do v=1, 28
4185 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4186 end do
4187 end do
4188 end do
4189 end do
4190 !---------------------------
4191 if ( is_top ) then
4192 do i=1, 42-1
4193 do v=1, 28
4194 ip = ipiv(v,i)
4195 if (ip /= i) then
4196 tmp = b(v,i,v2,ke_xy)
4197 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4198 b(v,ip,v2,ke_xy) = tmp
4199 end if
4200 end do
4201 end do
4202 do i=1, 42
4203 tmpv(:) = b(:,i,v2,ke_xy)
4204 do j=1, i-1
4205 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
4206 end do
4207 b(:,i,v2,ke_xy) = tmpv(:)
4208 end do
4209 do i=42, 1, -1
4210 tmpv(:) = b(:,i,v2,ke_xy)
4211 do j=i+1, 42
4212 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
4213 end do
4214 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
4215 end do
4216 else
4217 do i=1, 42-1
4218 do v=1, 28
4219 ip = ipiv(v,i)
4220 if (ip /= i) then
4221 tmp = b(v,i,v2,ke_xy)
4222 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4223 b(v,ip,v2,ke_xy) = tmp
4224
4225 tmp = g(v,i,1,v2,ke_xy)
4226 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
4227 g(v,ip,1,v2,ke_xy) = tmp
4228 tmp = g(v,i,2,v2,ke_xy)
4229 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
4230 g(v,ip,2,v2,ke_xy) = tmp
4231 tmp = g(v,i,3,v2,ke_xy)
4232 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
4233 g(v,ip,3,v2,ke_xy) = tmp
4234 end if
4235 end do
4236 end do
4237 do i=1, 42
4238 do v=1,28
4239 bb(v,1,i) = b(v,i,v2,ke_xy)
4240 bb(v,2,i) = g(v,i,1,v2,ke_xy)
4241 bb(v,3,i) = g(v,i,2,v2,ke_xy)
4242 bb(v,4,i) = g(v,i,3,v2,ke_xy)
4243 end do
4244 end do
4245 do i=1, 42
4246 tmpb(:,:) = bb(:,:,i)
4247 do j=1, i-1
4248 do r=1, 4
4249 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4250 end do
4251 end do
4252 bb(:,:,i) = tmpb(:,:)
4253 end do
4254 do i=42, 1, -1
4255 tmpb(:,:) = bb(:,:,i)
4256 do j=i+1, 42
4257 do r=1, 4
4258 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4259 end do
4260 end do
4261 do r=1, 4
4262 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4263 end do
4264 end do
4265 do i=1, 42
4266 b(:,i,v2,ke_xy) = bb(:,1,i)
4267 g(:,i,1,v2,ke_xy) = bb(:,2,i)
4268 g(:,i,2,v2,ke_xy) = bb(:,3,i)
4269 g(:,i,3,v2,ke_xy) = bb(:,4,i)
4270 end do
4271 end if
4272 end do
4273 end do
4274 return
4275 end subroutine solve_nnode14_var3
4276!OCL SERIAL
4277 subroutine solve_nnode15_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
4278 implicit none
4279 integer, intent(in) :: nnode_v
4280 integer, intent(in) :: jm
4281 integer, intent(in) :: ne2d
4282 real(rp), intent(in) :: d(15,15,15,jm,ne2d)
4283 real(rp), intent(inout) :: b(15,15,2,jm,ne2d)
4284 real(rp), intent(inout) :: g(15,15,jm,ne2d)
4285 logical, intent(in) :: is_top
4286
4287 integer :: ke_xy
4288 integer :: v2
4289 integer :: k, i, j, v, r
4290 integer :: ip
4291 integer :: ipiv(15,15)
4292 real(rp) :: tmpv(15)
4293 real(rp) :: tmpv2(15)
4294 real(rp) :: a(15,15,15)
4295 real(rp) :: bb(15,3,15)
4296 real(rp) :: tmpb(15,3)
4297 real(rp) :: tmp
4298 !------------------------------
4299
4300 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
4301 do ke_xy=1, ne2d
4302 do v2=1, jm
4303 a(:,:,:) = d(:,:,:,v2,ke_xy)
4304 do k=1, 15
4305 tmpv(:) = abs(a(:,k,k))
4306 ipiv(:,k) = k
4307 do i=k+1, 15
4308 do v=1, 15
4309 tmp = abs(a(v,i,k))
4310 if ( tmp > tmpv(v) ) then
4311 ipiv(v,k) = i
4312 tmpv(v) = tmp
4313 end if
4314 end do
4315 end do
4316 do v=1, 15
4317 ip = ipiv(v,k)
4318 if ( ip /= k ) then
4319 do j=1, 15
4320 tmp = a(v,k,j)
4321 a(v,k,j) = a(v,ip,j)
4322 a(v,ip,j) = tmp
4323 end do
4324 end if
4325 end do
4326 !--
4327 tmpv(:) = 1.0_rp / a(:,k,k)
4328 a(:,k,k) = tmpv(:)
4329 do i=k+1, 15
4330 a(:,i,k) = a(:,i,k) * tmpv(:)
4331 end do
4332 do j=k+1, 15
4333 do i=k+1, 15
4334 do v=1, 15
4335 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4336 end do
4337 end do
4338 end do
4339 end do
4340 !---------------------------
4341 if ( is_top ) then
4342 do i=1, 15-1
4343 do v=1, 15
4344 ip = ipiv(v,i)
4345 if (ip /= i) then
4346 tmp = b(v,i,1,v2,ke_xy)
4347 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4348 b(v,ip,1,v2,ke_xy) = tmp
4349 tmp = b(v,i,2,v2,ke_xy)
4350 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4351 b(v,ip,2,v2,ke_xy) = tmp
4352 end if
4353 end do
4354 end do
4355 do i=1, 15
4356 tmpv(:) = b(:,i,1,v2,ke_xy)
4357 tmpv2(:) = b(:,i,2,v2,ke_xy)
4358 do j=1, i-1
4359 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
4360 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
4361 end do
4362 b(:,i,1,v2,ke_xy) = tmpv(:)
4363 b(:,i,2,v2,ke_xy) = tmpv2(:)
4364 end do
4365 do i=15, 1, -1
4366 tmpv(:) = b(:,i,1,v2,ke_xy)
4367 tmpv2(:) = b(:,i,2,v2,ke_xy)
4368 do j=i+1, 15
4369 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
4370 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
4371 end do
4372 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
4373 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
4374 end do
4375 else
4376 do i=1, 15-1
4377 do v=1, 15
4378 ip = ipiv(v,i)
4379 if (ip /= i) then
4380 tmp = b(v,i,1,v2,ke_xy)
4381 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4382 b(v,ip,1,v2,ke_xy) = tmp
4383 tmp = b(v,i,2,v2,ke_xy)
4384 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4385 b(v,ip,2,v2,ke_xy) = tmp
4386
4387 tmp = g(v,i,v2,ke_xy)
4388 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
4389 g(v,ip,v2,ke_xy) = tmp
4390 end if
4391 end do
4392 end do
4393 do i=1, 15
4394 do v=1,15
4395 bb(v,1,i) = b(v,i,1,v2,ke_xy)
4396 bb(v,2,i) = b(v,i,2,v2,ke_xy)
4397 bb(v,3,i) = g(v,i,v2,ke_xy)
4398 end do
4399 end do
4400 do i=1, 15
4401 tmpb(:,:) = bb(:,:,i)
4402 do j=1, i-1
4403 do r=1, 3
4404 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4405 end do
4406 end do
4407 bb(:,:,i) = tmpb(:,:)
4408 end do
4409 do i=15, 1, -1
4410 tmpb(:,:) = bb(:,:,i)
4411 do j=i+1, 15
4412 do r=1, 3
4413 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4414 end do
4415 end do
4416 do r=1, 3
4417 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4418 end do
4419 end do
4420 do i=1, 15
4421 b(:,i,1,v2,ke_xy) = bb(:,1,i)
4422 b(:,i,2,v2,ke_xy) = bb(:,2,i)
4423 g(:,i,v2,ke_xy) = bb(:,3,i)
4424 end do
4425 end if
4426 end do
4427 end do
4428 return
4429 end subroutine solve_nnode15_uv
4430!OCL SERIAL
4431 subroutine solve_nnode15_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
4432 implicit none
4433 integer, intent(in) :: nnode_v
4434 integer, intent(in) :: jm
4435 integer, intent(in) :: ne2d
4436 real(rp), intent(in) :: d(15,45,45,jm,ne2d)
4437 real(rp), intent(inout) :: b(15,45,jm,ne2d)
4438 real(rp), intent(inout) :: g(15,45,3,jm,ne2d)
4439 logical, intent(in) :: is_top
4440
4441 integer :: ke_xy
4442 integer :: v2
4443 integer :: k, i, j, v, r
4444 integer :: ip
4445 integer :: ipiv(15,45)
4446 real(rp) :: tmpv(15)
4447 real(rp) :: a(15,45,45)
4448 real(rp) :: bb(15,4,45)
4449 real(rp) :: tmpb(15,4)
4450 real(rp) :: tmp
4451 !------------------------------
4452
4453 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
4454 do ke_xy=1, ne2d
4455 do v2=1, jm
4456 a(:,:,:) = d(:,:,:,v2,ke_xy)
4457 do k=1, 45
4458 tmpv(:) = abs(a(:,k,k))
4459 ipiv(:,k) = k
4460 do i=k+1, 45
4461 do v=1, 15
4462 tmp = abs(a(v,i,k))
4463 if ( tmp > tmpv(v) ) then
4464 ipiv(v,k) = i
4465 tmpv(v) = tmp
4466 end if
4467 end do
4468 end do
4469 do v=1, 15
4470 ip = ipiv(v,k)
4471 if ( ip /= k ) then
4472 do j=1, 45
4473 tmp = a(v,k,j)
4474 a(v,k,j) = a(v,ip,j)
4475 a(v,ip,j) = tmp
4476 end do
4477 end if
4478 end do
4479 !--
4480 tmpv(:) = 1.0_rp / a(:,k,k)
4481 a(:,k,k) = tmpv(:)
4482 do i=k+1, 45
4483 do v=1, 15
4484 a(v,i,k) = a(v,i,k) * tmpv(v)
4485 end do
4486 end do
4487 do j=k+1, 45
4488 do i=k+1, 45
4489 do v=1, 15
4490 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4491 end do
4492 end do
4493 end do
4494 end do
4495 !---------------------------
4496 if ( is_top ) then
4497 do i=1, 45-1
4498 do v=1, 15
4499 ip = ipiv(v,i)
4500 if (ip /= i) then
4501 tmp = b(v,i,v2,ke_xy)
4502 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4503 b(v,ip,v2,ke_xy) = tmp
4504 end if
4505 end do
4506 end do
4507 do i=1, 45
4508 tmpv(:) = b(:,i,v2,ke_xy)
4509 do j=1, i-1
4510 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
4511 end do
4512 b(:,i,v2,ke_xy) = tmpv(:)
4513 end do
4514 do i=45, 1, -1
4515 tmpv(:) = b(:,i,v2,ke_xy)
4516 do j=i+1, 45
4517 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
4518 end do
4519 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
4520 end do
4521 else
4522 do i=1, 45-1
4523 do v=1, 15
4524 ip = ipiv(v,i)
4525 if (ip /= i) then
4526 tmp = b(v,i,v2,ke_xy)
4527 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4528 b(v,ip,v2,ke_xy) = tmp
4529
4530 tmp = g(v,i,1,v2,ke_xy)
4531 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
4532 g(v,ip,1,v2,ke_xy) = tmp
4533 tmp = g(v,i,2,v2,ke_xy)
4534 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
4535 g(v,ip,2,v2,ke_xy) = tmp
4536 tmp = g(v,i,3,v2,ke_xy)
4537 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
4538 g(v,ip,3,v2,ke_xy) = tmp
4539 end if
4540 end do
4541 end do
4542 do i=1, 45
4543 do v=1,15
4544 bb(v,1,i) = b(v,i,v2,ke_xy)
4545 bb(v,2,i) = g(v,i,1,v2,ke_xy)
4546 bb(v,3,i) = g(v,i,2,v2,ke_xy)
4547 bb(v,4,i) = g(v,i,3,v2,ke_xy)
4548 end do
4549 end do
4550 do i=1, 45
4551 tmpb(:,:) = bb(:,:,i)
4552 do j=1, i-1
4553 do r=1, 4
4554 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4555 end do
4556 end do
4557 bb(:,:,i) = tmpb(:,:)
4558 end do
4559 do i=45, 1, -1
4560 tmpb(:,:) = bb(:,:,i)
4561 do j=i+1, 45
4562 do r=1, 4
4563 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4564 end do
4565 end do
4566 do r=1, 4
4567 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4568 end do
4569 end do
4570 do i=1, 45
4571 b(:,i,v2,ke_xy) = bb(:,1,i)
4572 g(:,i,1,v2,ke_xy) = bb(:,2,i)
4573 g(:,i,2,v2,ke_xy) = bb(:,3,i)
4574 g(:,i,3,v2,ke_xy) = bb(:,4,i)
4575 end do
4576 end if
4577 end do
4578 end do
4579 return
4580 end subroutine solve_nnode15_var3
4581!OCL SERIAL
4582 subroutine solve_nnode16_uv( D, b, G, Nnode_v, jm, Ne2D, is_top )
4583 implicit none
4584 integer, intent(in) :: nnode_v
4585 integer, intent(in) :: jm
4586 integer, intent(in) :: ne2d
4587 real(rp), intent(in) :: d(16,16,16,jm,ne2d)
4588 real(rp), intent(inout) :: b(16,16,2,jm,ne2d)
4589 real(rp), intent(inout) :: g(16,16,jm,ne2d)
4590 logical, intent(in) :: is_top
4591
4592 integer :: ke_xy
4593 integer :: v2
4594 integer :: k, i, j, v, r
4595 integer :: ip
4596 integer :: ipiv(16,16)
4597 real(rp) :: tmpv(16)
4598 real(rp) :: tmpv2(16)
4599 real(rp) :: a(16,16,16)
4600 real(rp) :: bb(16,3,16)
4601 real(rp) :: tmpb(16,3)
4602 real(rp) :: tmp
4603 !------------------------------
4604
4605 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
4606 do ke_xy=1, ne2d
4607 do v2=1, jm
4608 a(:,:,:) = d(:,:,:,v2,ke_xy)
4609 do k=1, 16
4610 tmpv(:) = abs(a(:,k,k))
4611 ipiv(:,k) = k
4612 do i=k+1, 16
4613 do v=1, 16
4614 tmp = abs(a(v,i,k))
4615 if ( tmp > tmpv(v) ) then
4616 ipiv(v,k) = i
4617 tmpv(v) = tmp
4618 end if
4619 end do
4620 end do
4621 do v=1, 16
4622 ip = ipiv(v,k)
4623 if ( ip /= k ) then
4624 do j=1, 16
4625 tmp = a(v,k,j)
4626 a(v,k,j) = a(v,ip,j)
4627 a(v,ip,j) = tmp
4628 end do
4629 end if
4630 end do
4631 !--
4632 tmpv(:) = 1.0_rp / a(:,k,k)
4633 a(:,k,k) = tmpv(:)
4634 do i=k+1, 16
4635 a(:,i,k) = a(:,i,k) * tmpv(:)
4636 end do
4637 do j=k+1, 16
4638 do i=k+1, 16
4639 do v=1, 16
4640 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4641 end do
4642 end do
4643 end do
4644 end do
4645 !---------------------------
4646 if ( is_top ) then
4647 do i=1, 16-1
4648 do v=1, 16
4649 ip = ipiv(v,i)
4650 if (ip /= i) then
4651 tmp = b(v,i,1,v2,ke_xy)
4652 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4653 b(v,ip,1,v2,ke_xy) = tmp
4654 tmp = b(v,i,2,v2,ke_xy)
4655 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4656 b(v,ip,2,v2,ke_xy) = tmp
4657 end if
4658 end do
4659 end do
4660 do i=1, 16
4661 tmpv(:) = b(:,i,1,v2,ke_xy)
4662 tmpv2(:) = b(:,i,2,v2,ke_xy)
4663 do j=1, i-1
4664 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
4665 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
4666 end do
4667 b(:,i,1,v2,ke_xy) = tmpv(:)
4668 b(:,i,2,v2,ke_xy) = tmpv2(:)
4669 end do
4670 do i=16, 1, -1
4671 tmpv(:) = b(:,i,1,v2,ke_xy)
4672 tmpv2(:) = b(:,i,2,v2,ke_xy)
4673 do j=i+1, 16
4674 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
4675 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
4676 end do
4677 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
4678 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
4679 end do
4680 else
4681 do i=1, 16-1
4682 do v=1, 16
4683 ip = ipiv(v,i)
4684 if (ip /= i) then
4685 tmp = b(v,i,1,v2,ke_xy)
4686 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4687 b(v,ip,1,v2,ke_xy) = tmp
4688 tmp = b(v,i,2,v2,ke_xy)
4689 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4690 b(v,ip,2,v2,ke_xy) = tmp
4691
4692 tmp = g(v,i,v2,ke_xy)
4693 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
4694 g(v,ip,v2,ke_xy) = tmp
4695 end if
4696 end do
4697 end do
4698 do i=1, 16
4699 do v=1,16
4700 bb(v,1,i) = b(v,i,1,v2,ke_xy)
4701 bb(v,2,i) = b(v,i,2,v2,ke_xy)
4702 bb(v,3,i) = g(v,i,v2,ke_xy)
4703 end do
4704 end do
4705 do i=1, 16
4706 tmpb(:,:) = bb(:,:,i)
4707 do j=1, i-1
4708 do r=1, 3
4709 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4710 end do
4711 end do
4712 bb(:,:,i) = tmpb(:,:)
4713 end do
4714 do i=16, 1, -1
4715 tmpb(:,:) = bb(:,:,i)
4716 do j=i+1, 16
4717 do r=1, 3
4718 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4719 end do
4720 end do
4721 do r=1, 3
4722 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4723 end do
4724 end do
4725 do i=1, 16
4726 b(:,i,1,v2,ke_xy) = bb(:,1,i)
4727 b(:,i,2,v2,ke_xy) = bb(:,2,i)
4728 g(:,i,v2,ke_xy) = bb(:,3,i)
4729 end do
4730 end if
4731 end do
4732 end do
4733 return
4734 end subroutine solve_nnode16_uv
4735!OCL SERIAL
4736 subroutine solve_nnode16_var3( D, b, G, Nnode_v, jm, Ne2D, is_top )
4737 implicit none
4738 integer, intent(in) :: nnode_v
4739 integer, intent(in) :: jm
4740 integer, intent(in) :: ne2d
4741 real(rp), intent(in) :: d(16,48,48,jm,ne2d)
4742 real(rp), intent(inout) :: b(16,48,jm,ne2d)
4743 real(rp), intent(inout) :: g(16,48,3,jm,ne2d)
4744 logical, intent(in) :: is_top
4745
4746 integer :: ke_xy
4747 integer :: v2
4748 integer :: k, i, j, v, r
4749 integer :: ip
4750 integer :: ipiv(16,48)
4751 real(rp) :: tmpv(16)
4752 real(rp) :: a(16,48,48)
4753 real(rp) :: bb(16,4,48)
4754 real(rp) :: tmpb(16,4)
4755 real(rp) :: tmp
4756 !------------------------------
4757
4758 !$omp parallel do collapse(2) private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A)
4759 do ke_xy=1, ne2d
4760 do v2=1, jm
4761 a(:,:,:) = d(:,:,:,v2,ke_xy)
4762 do k=1, 48
4763 tmpv(:) = abs(a(:,k,k))
4764 ipiv(:,k) = k
4765 do i=k+1, 48
4766 do v=1, 16
4767 tmp = abs(a(v,i,k))
4768 if ( tmp > tmpv(v) ) then
4769 ipiv(v,k) = i
4770 tmpv(v) = tmp
4771 end if
4772 end do
4773 end do
4774 do v=1, 16
4775 ip = ipiv(v,k)
4776 if ( ip /= k ) then
4777 do j=1, 48
4778 tmp = a(v,k,j)
4779 a(v,k,j) = a(v,ip,j)
4780 a(v,ip,j) = tmp
4781 end do
4782 end if
4783 end do
4784 !--
4785 tmpv(:) = 1.0_rp / a(:,k,k)
4786 a(:,k,k) = tmpv(:)
4787 do i=k+1, 48
4788 do v=1, 16
4789 a(v,i,k) = a(v,i,k) * tmpv(v)
4790 end do
4791 end do
4792 do j=k+1, 48
4793 do i=k+1, 48
4794 do v=1, 16
4795 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4796 end do
4797 end do
4798 end do
4799 end do
4800 !---------------------------
4801 if ( is_top ) then
4802 do i=1, 48-1
4803 do v=1, 16
4804 ip = ipiv(v,i)
4805 if (ip /= i) then
4806 tmp = b(v,i,v2,ke_xy)
4807 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4808 b(v,ip,v2,ke_xy) = tmp
4809 end if
4810 end do
4811 end do
4812 do i=1, 48
4813 tmpv(:) = b(:,i,v2,ke_xy)
4814 do j=1, i-1
4815 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
4816 end do
4817 b(:,i,v2,ke_xy) = tmpv(:)
4818 end do
4819 do i=48, 1, -1
4820 tmpv(:) = b(:,i,v2,ke_xy)
4821 do j=i+1, 48
4822 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
4823 end do
4824 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
4825 end do
4826 else
4827 do i=1, 48-1
4828 do v=1, 16
4829 ip = ipiv(v,i)
4830 if (ip /= i) then
4831 tmp = b(v,i,v2,ke_xy)
4832 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
4833 b(v,ip,v2,ke_xy) = tmp
4834
4835 tmp = g(v,i,1,v2,ke_xy)
4836 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
4837 g(v,ip,1,v2,ke_xy) = tmp
4838 tmp = g(v,i,2,v2,ke_xy)
4839 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
4840 g(v,ip,2,v2,ke_xy) = tmp
4841 tmp = g(v,i,3,v2,ke_xy)
4842 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
4843 g(v,ip,3,v2,ke_xy) = tmp
4844 end if
4845 end do
4846 end do
4847 do i=1, 48
4848 do v=1,16
4849 bb(v,1,i) = b(v,i,v2,ke_xy)
4850 bb(v,2,i) = g(v,i,1,v2,ke_xy)
4851 bb(v,3,i) = g(v,i,2,v2,ke_xy)
4852 bb(v,4,i) = g(v,i,3,v2,ke_xy)
4853 end do
4854 end do
4855 do i=1, 48
4856 tmpb(:,:) = bb(:,:,i)
4857 do j=1, i-1
4858 do r=1, 4
4859 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
4860 end do
4861 end do
4862 bb(:,:,i) = tmpb(:,:)
4863 end do
4864 do i=48, 1, -1
4865 tmpb(:,:) = bb(:,:,i)
4866 do j=i+1, 48
4867 do r=1, 4
4868 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
4869 end do
4870 end do
4871 do r=1, 4
4872 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
4873 end do
4874 end do
4875 do i=1, 48
4876 b(:,i,v2,ke_xy) = bb(:,1,i)
4877 g(:,i,1,v2,ke_xy) = bb(:,2,i)
4878 g(:,i,2,v2,ke_xy) = bb(:,3,i)
4879 g(:,i,3,v2,ke_xy) = bb(:,4,i)
4880 end do
4881 end if
4882 end do
4883 end do
4884 return
4885 end subroutine solve_nnode16_var3
4886!OCL SERIAL
4887 subroutine solve_uv_im( D, b, G, Nnode_v, im, jm, Ne2D, is_top )
4888 implicit none
4889 integer, intent(in) :: nnode_v
4890 integer, intent(in) :: im, jm
4891 integer, intent(in) :: ne2d
4892 real(rp), intent(in) :: d(im,nnode_v,nnode_v,jm,ne2d)
4893 real(rp), intent(inout) :: b(im,nnode_v,2,jm,ne2d)
4894 real(rp), intent(inout) :: g(im,nnode_v,jm,ne2d)
4895 logical, intent(in) :: is_top
4896
4897 integer :: ke_xy
4898 integer :: v2
4899 integer :: k, i, j, v, r
4900 integer :: ip
4901 integer :: ipiv(im,nnode_v)
4902 real(rp) :: tmpv(im)
4903 real(rp) :: tmpv2(im)
4904 real(rp) :: a(im,nnode_v,nnode_v)
4905 real(rp) :: bb(im,3,nnode_v)
4906 real(rp) :: tmpb(im,3)
4907 real(rp) :: tmp
4908 !------------------------------
4909
4910 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpv2,tmpb,tmp,bb,A) collapse(2)
4911 do ke_xy=1, ne2d
4912 do v2=1, jm
4913 a(:,:,:) = d(:,:,:,v2,ke_xy)
4914 do k=1, nnode_v
4915 tmpv(:) = abs(a(:,k,k))
4916 ipiv(:,k) = k
4917 do i=k+1, nnode_v
4918 do v=1, im
4919 tmp = abs(a(v,i,k))
4920 if ( tmp > tmpv(v) ) then
4921 ipiv(v,k) = i
4922 tmpv(v) = tmp
4923 end if
4924 end do
4925 end do
4926 do v=1, im
4927 ip = ipiv(v,k)
4928 if ( ip /= k ) then
4929 do j=1, nnode_v
4930 tmp = a(v,k,j)
4931 a(v,k,j) = a(v,ip,j)
4932 a(v,ip,j) = tmp
4933 end do
4934 end if
4935 end do
4936 !--
4937 tmpv(:) = 1.0_rp / a(:,k,k)
4938 a(:,k,k) = tmpv(:)
4939 do i=k+1, nnode_v
4940 a(:,i,k) = a(:,i,k) * tmpv(:)
4941 end do
4942 do j=k+1, nnode_v
4943 do i=k+1, nnode_v
4944 do v=1, im
4945 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
4946 end do
4947 end do
4948 end do
4949 end do
4950 !---------------------------
4951 if ( is_top ) then
4952 do i=1, nnode_v-1
4953 do v=1, im
4954 ip = ipiv(v,i)
4955 if (ip /= i) then
4956 tmp = b(v,i,1,v2,ke_xy)
4957 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4958 b(v,ip,1,v2,ke_xy) = tmp
4959 tmp = b(v,i,2,v2,ke_xy)
4960 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4961 b(v,ip,2,v2,ke_xy) = tmp
4962 end if
4963 end do
4964 end do
4965 do i=1, nnode_v
4966 tmpv(:) = b(:,i,1,v2,ke_xy)
4967 tmpv2(:) = b(:,i,2,v2,ke_xy)
4968 do j=1, i-1
4969 tmpv(:) = tmpv(:) - b(:,j,1,v2,ke_xy) * a(:,i,j)
4970 tmpv2(:) = tmpv2(:) - b(:,j,2,v2,ke_xy) * a(:,i,j)
4971 end do
4972 b(:,i,1,v2,ke_xy) = tmpv(:)
4973 b(:,i,2,v2,ke_xy) = tmpv2(:)
4974 end do
4975 do i=nnode_v, 1, -1
4976 tmpv(:) = b(:,i,1,v2,ke_xy)
4977 tmpv2(:) = b(:,i,2,v2,ke_xy)
4978 do j=i+1, nnode_v
4979 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,1,v2,ke_xy)
4980 tmpv2(:) = tmpv2(:) - a(:,i,j) * b(:,j,2,v2,ke_xy)
4981 end do
4982 b(:,i,1,v2,ke_xy) = tmpv(:) * a(:,i,i)
4983 b(:,i,2,v2,ke_xy) = tmpv2(:) * a(:,i,i)
4984 end do
4985 else
4986 do i=1, nnode_v-1
4987 do v=1, im
4988 ip = ipiv(v,i)
4989 if (ip /= i) then
4990 tmp = b(v,i,1,v2,ke_xy)
4991 b(v,i,1,v2,ke_xy) = b(v,ip,1,v2,ke_xy)
4992 b(v,ip,1,v2,ke_xy) = tmp
4993 tmp = b(v,i,2,v2,ke_xy)
4994 b(v,i,2,v2,ke_xy) = b(v,ip,2,v2,ke_xy)
4995 b(v,ip,2,v2,ke_xy) = tmp
4996
4997 tmp = g(v,i,v2,ke_xy)
4998 g(v,i,v2,ke_xy) = g(v,ip,v2,ke_xy)
4999 g(v,ip,v2,ke_xy) = tmp
5000 end if
5001 end do
5002 end do
5003 do i=1, nnode_v
5004 do v=1,im
5005 bb(v,1,i) = b(v,i,1,v2,ke_xy)
5006 bb(v,2,i) = b(v,i,2,v2,ke_xy)
5007 bb(v,3,i) = g(v,i,v2,ke_xy)
5008 end do
5009 end do
5010 do i=1, nnode_v
5011 tmpb(:,:) = bb(:,:,i)
5012 do j=1, i-1
5013 do r=1, 3
5014 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
5015 end do
5016 end do
5017 bb(:,:,i) = tmpb(:,:)
5018 end do
5019 do i=nnode_v, 1, -1
5020 tmpb(:,:) = bb(:,:,i)
5021 do j=i+1, nnode_v
5022 do r=1, 3
5023 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
5024 end do
5025 end do
5026 do r=1, 3
5027 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
5028 end do
5029 end do
5030 do i=1, nnode_v
5031 b(:,i,1,v2,ke_xy) = bb(:,1,i)
5032 b(:,i,2,v2,ke_xy) = bb(:,2,i)
5033 g(:,i,v2,ke_xy) = bb(:,3,i)
5034 end do
5035 end if
5036 end do
5037 end do
5038 return
5039 end subroutine solve_uv_im
5040
5041!OCL SERIAL
5042 subroutine solve_var3_im( D, b, G, Nnode_v, im, jm, Ne2D, is_top )
5043 implicit none
5044 integer, intent(in) :: nnode_v
5045 integer, intent(in) :: im, jm
5046 integer, intent(in) :: ne2d
5047 real(rp), intent(in) :: d(im,3*nnode_v,3*nnode_v,jm,ne2d)
5048 real(rp), intent(inout) :: b(im,3*nnode_v,jm,ne2d)
5049 real(rp), intent(inout) :: g(im,3*nnode_v,3,jm,ne2d)
5050 logical, intent(in) :: is_top
5051
5052 integer :: ke_xy
5053 integer :: v2
5054 integer :: k, i, j, v, r
5055 integer :: ip
5056 integer :: ipiv(im,3*nnode_v)
5057 real(rp) :: tmpv(im)
5058 real(rp) :: a(im,3*nnode_v,3*nnode_v)
5059 real(rp) :: bb(im,4,3*nnode_v)
5060 real(rp) :: tmpb(im,4)
5061 real(rp) :: tmp
5062 integer :: nnode_vx3
5063 !------------------------------
5064
5065 nnode_vx3 = nnode_v * 3
5066
5067 !$omp parallel do private(ke_xy,v2,k,i,j,v,r,ip,ipiv,tmpv,tmpb,tmp,bb,A) collapse(2)
5068 do ke_xy=1, ne2d
5069 do v2=1, jm
5070 a(:,:,:) = d(:,:,:,v2,ke_xy)
5071 do k=1, nnode_vx3
5072 tmpv(:) = abs(a(:,k,k))
5073 ipiv(:,k) = k
5074 do i=k+1, nnode_vx3
5075 do v=1, im
5076 tmp = abs(a(v,i,k))
5077 if ( tmp > tmpv(v) ) then
5078 ipiv(v,k) = i
5079 tmpv(v) = tmp
5080 end if
5081 end do
5082 end do
5083 do v=1, im
5084 ip = ipiv(v,k)
5085 if ( ip /= k ) then
5086 do j=1, nnode_vx3
5087 tmp = a(v,k,j)
5088 a(v,k,j) = a(v,ip,j)
5089 a(v,ip,j) = tmp
5090 end do
5091 end if
5092 end do
5093 !--
5094 tmpv(:) = 1.0_rp / a(:,k,k)
5095 a(:,k,k) = tmpv(:)
5096 do i=k+1, nnode_vx3
5097 a(:,i,k) = a(:,i,k) * tmpv(:)
5098 end do
5099 do j=k+1, nnode_vx3
5100 do i=k+1, nnode_vx3
5101 do v=1, im
5102 a(v,i,j) = a(v,i,j) - a(v,i,k) * a(v,k,j)
5103 end do
5104 end do
5105 end do
5106 end do
5107 !---------------------------
5108 if ( is_top ) then
5109 do i=1, nnode_vx3-1
5110 do v=1, im
5111 ip = ipiv(v,i)
5112 if (ip /= i) then
5113 tmp = b(v,i,v2,ke_xy)
5114 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
5115 b(v,ip,v2,ke_xy) = tmp
5116 end if
5117 end do
5118 end do
5119 do i=1, nnode_vx3
5120 tmpv(:) = b(:,i,v2,ke_xy)
5121 do j=1, i-1
5122 tmpv(:) = tmpv(:) - b(:,j,v2,ke_xy) * a(:,i,j)
5123 end do
5124 b(:,i,v2,ke_xy) = tmpv(:)
5125 end do
5126 do i=nnode_vx3, 1, -1
5127 tmpv(:) = b(:,i,v2,ke_xy)
5128 do j=i+1, nnode_vx3
5129 tmpv(:) = tmpv(:) - a(:,i,j) * b(:,j,v2,ke_xy)
5130 end do
5131 b(:,i,v2,ke_xy) = tmpv(:) * a(:,i,i)
5132 end do
5133 else
5134 do i=1, nnode_vx3-1
5135 do v=1, im
5136 ip = ipiv(v,i)
5137 if (ip /= i) then
5138 tmp = b(v,i,v2,ke_xy)
5139 b(v,i,v2,ke_xy) = b(v,ip,v2,ke_xy)
5140 b(v,ip,v2,ke_xy) = tmp
5141
5142 tmp = g(v,i,1,v2,ke_xy)
5143 g(v,i,1,v2,ke_xy) = g(v,ip,1,v2,ke_xy)
5144 g(v,ip,1,v2,ke_xy) = tmp
5145 tmp = g(v,i,2,v2,ke_xy)
5146 g(v,i,2,v2,ke_xy) = g(v,ip,2,v2,ke_xy)
5147 g(v,ip,2,v2,ke_xy) = tmp
5148 tmp = g(v,i,3,v2,ke_xy)
5149 g(v,i,3,v2,ke_xy) = g(v,ip,3,v2,ke_xy)
5150 g(v,ip,3,v2,ke_xy) = tmp
5151 end if
5152 end do
5153 end do
5154 do i=1, nnode_vx3
5155 do v=1, im
5156 bb(v,1,i) = b(v,i,v2,ke_xy)
5157 bb(v,2,i) = g(v,i,1,v2,ke_xy)
5158 bb(v,3,i) = g(v,i,2,v2,ke_xy)
5159 bb(v,4,i) = g(v,i,3,v2,ke_xy)
5160 end do
5161 end do
5162 do i=1, nnode_vx3
5163 tmpb(:,:) = bb(:,:,i)
5164 do j=1, i-1
5165 do r=1, 4
5166 tmpb(:,r) = tmpb(:,r) - bb(:,r,j) * a(:,i,j)
5167 end do
5168 end do
5169 bb(:,:,i) = tmpb(:,:)
5170 end do
5171 do i=nnode_vx3, 1, -1
5172 tmpb(:,:) = bb(:,:,i)
5173 do j=i+1, nnode_vx3
5174 do r=1, 4
5175 tmpb(:,r) = tmpb(:,r) - a(:,i,j) * bb(:,r,j)
5176 end do
5177 end do
5178 do r=1, 4
5179 bb(:,r,i) = tmpb(:,r) * a(:,i,i)
5180 end do
5181 end do
5182 do i=1, nnode_vx3
5183 b(:,i,v2,ke_xy) = bb(:,1,i)
5184 g(:,i,1,v2,ke_xy) = bb(:,2,i)
5185 g(:,i,2,v2,ke_xy) = bb(:,3,i)
5186 g(:,i,3,v2,ke_xy) = bb(:,4,i)
5187 end do
5188 end if
5189 end do
5190 end do
5191 return
5192 end subroutine solve_var3_im
5193
5194end module scale_atm_dyn_dgm_hevi_common_linalgebra
scale_atm_dyn_dgm_hevi_common_linalgebra
module FElib / Fluid dyn solver / Atmosphere / HEVI / Common
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:17
scale_atm_dyn_dgm_hevi_common_linalgebra::atm_dyn_dgm_hevi_common_linalgebra_solve_var3
subroutine, public atm_dyn_dgm_hevi_common_linalgebra_solve_var3(d, b, g, nnode_v, im, jm, ne2d, is_top)
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:203
scale_atm_dyn_dgm_hevi_common_linalgebra::atm_dyn_dgm_hevi_common_linalgebra_solve_uv
subroutine, public atm_dyn_dgm_hevi_common_linalgebra_solve_uv(d, b, g, nnode_v, im, jm, ne2d, is_top)
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:94
scale_atm_dyn_dgm_hevi_common_linalgebra::atm_dyn_dgm_hevi_common_linalgebra_get_param
subroutine, public atm_dyn_dgm_hevi_common_linalgebra_get_param(im, jm, nnode_h1d)
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:50
Generated by doxygen 1.14.0