FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
hecmw_precond_SSOR_44.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
6 !C
7 !C***
8 !C*** module hecmw_precond_SSOR_44
9 !C***
10 !C
12  use hecmw_util
17 #ifndef _OPENACC
18  !$ use omp_lib
19 #endif
20 
21  private
22 
26 
27  integer(kind=kint) :: N
28  real(kind=kreal), pointer :: d(:) => null()
29  real(kind=kreal), pointer :: al(:) => null()
30  real(kind=kreal), pointer :: au(:) => null()
31  integer(kind=kint), pointer :: indexL(:) => null()
32  integer(kind=kint), pointer :: indexU(:) => null()
33  integer(kind=kint), pointer :: itemL(:) => null()
34  integer(kind=kint), pointer :: itemU(:) => null()
35  real(kind=kreal), pointer :: alu(:) => null()
36 
37  integer(kind=kint) :: NColor
38  integer(kind=kint), pointer :: COLORindex(:) => null()
39  integer(kind=kint), pointer :: perm(:) => null()
40  integer(kind=kint), pointer :: iperm(:) => null()
41 
42  logical, save :: isFirst = .true.
43 
44  logical, save :: INITIALIZED = .false.
45 
46  ! relaxation parameter of
47  ! M = 1/(2-OMEGA) (D/OMEGA + L) (D/OMEGA)^-1 (D/OMEGA + U)
48  ! OMEGA = 1 is the classic (D+L) D^-1 (D+U). _apply is handed only ZP, so setup
49  ! leaves the value here for it.
50  real(kind=kreal), save :: omega = 1.d0
51 
52  ! for tuning
53  integer(kind=kint), parameter :: numOfBlockPerThread = 100
54  integer(kind=kint), save :: numOfThread = 1, numofblock
55  integer(kind=kint), save, allocatable :: icToBlockIndex(:)
56  integer(kind=kint), save, allocatable :: blockIndexToColorIndex(:)
57  integer(kind=kint), save :: sectorCacheSize0, sectorCacheSize1
58 
59  integer(kind=kint), parameter :: SSOR_DEBUG = 0
60 
61 contains
62 
63  subroutine hecmw_precond_ssor_44_setup(hecMAT)
64  implicit none
65  type(hecmwst_matrix), intent(inout) :: hecmat
66  integer(kind=kint ) :: npl, npu
67  integer(kind=kint ) :: ncolor_in
68  real (kind=kreal) :: omega_inv
69  real (kind=kreal) :: alutmp(4,4), pw(4)
70  integer(kind=kint ) :: ii, i, j, k
71  integer(kind=kint ) :: nthreads = 1
72  integer(kind=kint ), allocatable :: perm_tmp(:)
73  real (kind=kreal) :: t0
74 
75  if (ssor_debug >= 1) then
76  t0 = hecmw_wtime()
77  write(*,*) 'DEBUG: SSOR setup start', hecmw_wtime()-t0
78  endif
79 
80  if (initialized) then
81  if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
82  call hecmw_precond_ssor_44_clear(hecmat)
83  else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
84  call hecmw_precond_ssor_44_clear(hecmat) ! TEMPORARY
85  else
86  return
87  endif
88  endif
89 
90 #ifndef _OPENACC
91  !$ nthreads = omp_get_max_threads()
92 #endif
93 
94  n = hecmat%N
95  ! N = hecMAT%NP
96  ncolor_in = hecmw_mat_get_ncolor_in(hecmat)
97  omega = hecmw_mat_get_omega(hecmat)
98  omega_inv = 1.d0 / omega
99 
100 #ifdef _OPENACC
101  allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
102  call hecmw_matrix_ordering_rcm(n, hecmat%indexL, hecmat%itemL, &
103  hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
104  if (ssor_debug >= 1) write(*,*) 'DEBUG: RCM ordering done', hecmw_wtime()-t0
105  call hecmw_matrix_ordering_mc(n, hecmat%indexL, hecmat%itemL, &
106  hecmat%indexU, hecmat%itemU, perm_tmp, &
107  ncolor_in, ncolor, colorindex, perm, iperm)
108  if (ssor_debug >= 1) write(*,*) 'DEBUG: MC ordering done', hecmw_wtime()-t0
109  deallocate(perm_tmp)
110 
111 #else
112  if (nthreads == 1) then
113  ncolor = 1
114  allocate(colorindex(0:1), perm(n), iperm(n))
115  colorindex(0) = 0
116  colorindex(1) = n
117  do i=1,n
118  perm(i) = i
119  iperm(i) = i
120  end do
121  else
122  allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
123  call hecmw_matrix_ordering_rcm(n, hecmat%indexL, hecmat%itemL, &
124  hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
125  if (ssor_debug >= 1) write(*,*) 'DEBUG: RCM ordering done', hecmw_wtime()-t0
126  call hecmw_matrix_ordering_mc(n, hecmat%indexL, hecmat%itemL, &
127  hecmat%indexU, hecmat%itemU, perm_tmp, &
128  ncolor_in, ncolor, colorindex, perm, iperm)
129  if (ssor_debug >= 1) write(*,*) 'DEBUG: MC ordering done', hecmw_wtime()-t0
130  deallocate(perm_tmp)
131 
132  endif
133 #endif
134 
135  npl = 0
136  do i=1,n
137  do j=hecmat%indexU(i-1)+1,hecmat%indexU(i)
138  if( hecmat%itemU(j) > n ) exit
139  npl = npl + 1
140  enddo
141  enddo
142  npl = max(hecmat%indexL(n),npl)
143  npu = hecmat%indexU(n)
144  allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
145  call hecmw_matrix_reorder_profile(n, perm, iperm, &
146  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
147  indexl, indexu, iteml, itemu)
148  if (ssor_debug >= 1) write(*,*) 'DEBUG: reordering profile done', hecmw_wtime()-t0
149 
150 
151  allocate(d(16*n), al(16*npl), au(16*npu))
152  call hecmw_matrix_reorder_values(n, 4, perm, iperm, &
153  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
154  hecmat%AL, hecmat%AU, hecmat%D, &
155  indexl, indexu, iteml, itemu, al, au, d)
156  if (ssor_debug >= 1) write(*,*) 'DEBUG: reordering values done', hecmw_wtime()-t0
157 
158  call hecmw_matrix_reorder_renum_item(n, perm, indexl, iteml)
159  call hecmw_matrix_reorder_renum_item(n, perm, indexu, itemu)
160 
161  allocate(alu(16*n))
162  alu = 0.d0
163 
164  do ii= 1, 16*n
165  alu(ii) = d(ii)
166  enddo
167 
168 #ifdef _OPENACC
169  !$acc kernels
170  !$acc loop independent private(ALUtmp,PW)
171 #else
172  !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,OMEGA_INV)
173  !$omp do
174 #endif
175  do ii= 1, n
176  alutmp(1,1)= alu(16*ii-15) * omega_inv
177  alutmp(1,2)= alu(16*ii-14) * omega_inv
178  alutmp(1,3)= alu(16*ii-13) * omega_inv
179  alutmp(1,4)= alu(16*ii-12) * omega_inv
180  alutmp(2,1)= alu(16*ii-11) * omega_inv
181  alutmp(2,2)= alu(16*ii-10) * omega_inv
182  alutmp(2,3)= alu(16*ii- 9) * omega_inv
183  alutmp(2,4)= alu(16*ii- 8) * omega_inv
184  alutmp(3,1)= alu(16*ii- 7) * omega_inv
185  alutmp(3,2)= alu(16*ii- 6) * omega_inv
186  alutmp(3,3)= alu(16*ii- 5) * omega_inv
187  alutmp(3,4)= alu(16*ii- 4) * omega_inv
188  alutmp(4,1)= alu(16*ii- 3) * omega_inv
189  alutmp(4,2)= alu(16*ii- 2) * omega_inv
190  alutmp(4,3)= alu(16*ii- 1) * omega_inv
191  alutmp(4,4)= alu(16*ii ) * omega_inv
192 
193  do k= 1, 4
194  alutmp(k,k)= 1.d0/alutmp(k,k)
195  do i= k+1, 4
196  alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
197  do j= k+1, 4
198  pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
199  enddo
200  do j= k+1, 4
201  alutmp(i,j)= pw(j)
202  enddo
203  enddo
204  enddo
205  alu(16*ii-15)= alutmp(1,1)
206  alu(16*ii-14)= alutmp(1,2)
207  alu(16*ii-13)= alutmp(1,3)
208  alu(16*ii-12)= alutmp(1,4)
209  alu(16*ii-11)= alutmp(2,1)
210  alu(16*ii-10)= alutmp(2,2)
211  alu(16*ii- 9)= alutmp(2,3)
212  alu(16*ii- 8)= alutmp(2,4)
213  alu(16*ii- 7)= alutmp(3,1)
214  alu(16*ii- 6)= alutmp(3,2)
215  alu(16*ii- 5)= alutmp(3,3)
216  alu(16*ii- 4)= alutmp(3,4)
217  alu(16*ii- 3)= alutmp(4,1)
218  alu(16*ii- 2)= alutmp(4,2)
219  alu(16*ii- 1)= alutmp(4,3)
220  alu(16*ii )= alutmp(4,4)
221 
222  enddo
223 #ifdef _OPENACC
224  !$acc end kernels
225 #else
226  !$omp end do
227  !$omp end parallel
228 #endif
229 
230  isfirst = .true.
231 
232  initialized = .true.
233  hecmat%Iarray(98) = 0 ! symbolic setup done
234  hecmat%Iarray(97) = 0 ! numerical setup done
235 
236  if (ssor_debug >= 1) write(*,*) 'DEBUG: SSOR setup done', hecmw_wtime()-t0
237 
238  end subroutine hecmw_precond_ssor_44_setup
239 
240  subroutine setup_tuning_parameters
241  use hecmw_tuning_fx
242  implicit none
243  integer(kind=kint) :: blockindex, elementcount, numofelement, ii
244  real(kind=kreal) :: numofelementperblock
245  integer(kind=kint) :: ic, i
246  if (ssor_debug >= 1) write(*,*) 'DEBUG: setting up tuning parameters for SSOR'
247 #ifndef _OPENACC
248  !$ numOfThread = omp_get_max_threads()
249 #endif
250 
251  numofblock = numofthread * numofblockperthread
252  if (allocated(ictoblockindex)) deallocate(ictoblockindex)
253  if (allocated(blockindextocolorindex)) deallocate(blockindextocolorindex)
254  allocate (ictoblockindex(0:ncolor), &
255  blockindextocolorindex(0:numofblock + ncolor))
256  numofelement = n + indexl(n) + indexu(n)
257  numofelementperblock = dble(numofelement) / numofblock
258  blockindex = 0
259  ictoblockindex = -1
260  ictoblockindex(0) = 0
261  blockindextocolorindex = -1
262  blockindextocolorindex(0) = 0
263  do ic = 1, ncolor
264  elementcount = 0
265  ii = 1
266  do i = colorindex(ic-1)+1, colorindex(ic)
267  elementcount = elementcount + 1
268  elementcount = elementcount + (indexl(i) - indexl(i-1))
269  elementcount = elementcount + (indexu(i) - indexu(i-1))
270  if (elementcount > ii * numofelementperblock &
271  .or. i == colorindex(ic)) then
272  ii = ii + 1
273  blockindex = blockindex + 1
274  blockindextocolorindex(blockindex) = i
275  endif
276  enddo
277  ictoblockindex(ic) = blockindex
278  enddo
279  numofblock = blockindex
280 
282  sectorcachesize0, sectorcachesize1 )
283  end subroutine setup_tuning_parameters
284 
286  implicit none
287  real(kind=kreal), intent(inout) :: zp(:)
288  integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
289  real(kind=kreal) :: sw1, sw2, sw3, sw4, x1, x2, x3, x4
290 
291  real(kind=kreal) :: omega_fac
292 
293  ! added for tuning >>>
294  integer(kind=kint) :: blockindex
295 
296 #ifndef _OPENACC
297  if (isfirst) then
298  call setup_tuning_parameters
299  isfirst = .false.
300  endif
301 #endif
302  ! <<< added for tuning
303 
304  omega_fac = 2.d0 - omega
305 
306 #ifndef _OPENACC
307  !call start_collection("loopInPrecond44")
308 
309  !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
310  !OCL CACHE_SUBSECTOR_ASSIGN(ZP)
311 
312  !$omp parallel default(none) &
313  !$omp&shared(NColor,indexL,itemL,indexU,itemU,AL,AU,D,ALU,perm,&
314  !$omp& ZP,icToBlockIndex,blockIndexToColorIndex,OMEGA_FAC) &
315  !$omp&private(SW1,SW2,SW3,SW4,X1,X2,X3,X4,ic,i,iold,isL,ieL,isU,ieU,j,k,blockIndex)
316 #endif
317 
318  !C-- FORWARD
319  do ic=1,ncolor
320 #ifdef _OPENACC
321  !$acc kernels
322  !$acc loop independent
323  do i = colorindex(ic-1)+1, colorindex(ic)
324 #else
325  !$omp do schedule (static, 1)
326  do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
327  do i = blockindextocolorindex(blockindex-1)+1, &
328  blockindextocolorindex(blockindex)
329 #endif
330  ! do i = startPos(threadNum, ic), endPos(threadNum, ic)
331  iold = perm(i)
332  sw1= omega_fac * zp(4*iold-3)
333  sw2= omega_fac * zp(4*iold-2)
334  sw3= omega_fac * zp(4*iold-1)
335  sw4= omega_fac * zp(4*iold )
336  isl= indexl(i-1)+1
337  iel= indexl(i)
338  do j= isl, iel
339  !k= perm(itemL(j))
340  k= iteml(j)
341  x1= zp(4*k-3)
342  x2= zp(4*k-2)
343  x3= zp(4*k-1)
344  x4= zp(4*k )
345  sw1= sw1 - al(16*j-15)*x1 - al(16*j-14)*x2 - al(16*j-13)*x3 - al(16*j-12)*x4
346  sw2= sw2 - al(16*j-11)*x1 - al(16*j-10)*x2 - al(16*j- 9)*x3 - al(16*j- 8)*x4
347  sw3= sw3 - al(16*j- 7)*x1 - al(16*j- 6)*x2 - al(16*j- 5)*x3 - al(16*j- 4)*x4
348  sw4= sw4 - al(16*j- 3)*x1 - al(16*j- 2)*x2 - al(16*j- 1)*x3 - al(16*j- 0)*x4
349  enddo ! j
350 
351  x1= sw1
352  x2= sw2
353  x3= sw3
354  x4= sw4
355  x2= x2 - alu(16*i-11)*x1
356  x3= x3 - alu(16*i- 7)*x1 - alu(16*i- 6)*x2
357  x4= x4 - alu(16*i- 3)*x1 - alu(16*i- 2)*x2 - alu(16*i- 1)*x3
358 
359  x4= alu(16*i )* x4
360  x3= alu(16*i- 5)*( x3 - alu(16*i- 4)*x4 )
361  x2= alu(16*i-10)*( x2 - alu(16*i- 8)*x4 - alu(16*i- 9)*x3 )
362  x1= alu(16*i-15)*( x1 - alu(16*i-12)*x4 - alu(16*i-13)*x3 - alu(16*i-14)*x2 )
363 
364  zp(4*iold-3)= x1
365  zp(4*iold-2)= x2
366  zp(4*iold-1)= x3
367  zp(4*iold )= x4
368 #ifdef _OPENACC
369  enddo
370  !$acc end kernels
371 #else
372  enddo ! i
373  enddo ! blockIndex
374  !$omp end do
375 #endif
376  enddo ! ic
377 
378  !C-- BACKWARD
379  do ic=ncolor, 1, -1
380 #ifdef _OPENACC
381  !$acc kernels
382  !$acc loop independent
383  do i = colorindex(ic-1)+1, colorindex(ic)
384 #else
385  !$omp do schedule (static, 1)
386  do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
387  do i = blockindextocolorindex(blockindex), &
388  blockindextocolorindex(blockindex-1)+1, -1
389 #endif
390  ! do blockIndex = icToBlockIndex(ic-1)+1, icToBlockIndex(ic)
391  ! do i = blockIndexToColorIndex(blockIndex-1)+1, &
392  ! blockIndexToColorIndex(blockIndex)
393  ! do i = endPos(threadNum, ic), startPos(threadNum, ic), -1
394  sw1= 0.d0
395  sw2= 0.d0
396  sw3= 0.d0
397  sw4= 0.d0
398  isu= indexu(i-1) + 1
399  ieu= indexu(i)
400  do j= ieu, isu, -1
401  !k= perm(itemU(j))
402  k= itemu(j)
403  x1= zp(4*k-3)
404  x2= zp(4*k-2)
405  x3= zp(4*k-1)
406  x4= zp(4*k )
407  sw1= sw1 + au(16*j-15)*x1 + au(16*j-14)*x2 + au(16*j-13)*x3 + au(16*j-12)*x4
408  sw2= sw2 + au(16*j-11)*x1 + au(16*j-10)*x2 + au(16*j- 9)*x3 + au(16*j- 8)*x4
409  sw3= sw3 + au(16*j- 7)*x1 + au(16*j- 6)*x2 + au(16*j- 5)*x3 + au(16*j- 4)*x4
410  sw4= sw4 + au(16*j- 3)*x1 + au(16*j- 2)*x2 + au(16*j- 1)*x3 + au(16*j- 0)*x4
411 
412  enddo ! j
413 
414  x1= sw1
415  x2= sw2
416  x3= sw3
417  x4= sw4
418  x2= x2 - alu(16*i-11)*x1
419  x3= x3 - alu(16*i- 7)*x1 - alu(16*i- 6)*x2
420  x4= x4 - alu(16*i- 3)*x1 - alu(16*i- 2)*x2 - alu(16*i- 1)*x3
421  x4= alu(16*i )* x4
422  x3= alu(16*i- 5)*( x3 - alu(16*i- 4)*x4 )
423  x2= alu(16*i-10)*( x2 - alu(16*i- 8)*x4 - alu(16*i- 9)*x3 )
424  x1= alu(16*i-15)*( x1 - alu(16*i-12)*x4 - alu(16*i-13)*x3 - alu(16*i-14)*x2 )
425  iold = perm(i)
426  zp(4*iold-3)= zp(4*iold-3) - x1
427  zp(4*iold-2)= zp(4*iold-2) - x2
428  zp(4*iold-1)= zp(4*iold-1) - x3
429  zp(4*iold )= zp(4*iold ) - x4
430 #ifdef _OPENACC
431  enddo
432  !$acc end kernels
433 #else
434  enddo ! i
435  enddo ! blockIndex
436  !$omp end do
437 #endif
438  enddo ! ic
439 #ifndef _OPENACC
440  !$omp end parallel
441 
442  !OCL END_CACHE_SUBSECTOR
443  !OCL END_CACHE_SECTOR_SIZE
444 
445  !call stop_collection("loopInPrecond44")
446 #endif
447 
448  end subroutine hecmw_precond_ssor_44_apply
449 
450  subroutine hecmw_precond_ssor_44_clear(hecMAT)
451  implicit none
452  type(hecmwst_matrix), intent(inout) :: hecmat
453  if (associated(colorindex)) deallocate(colorindex)
454  if (associated(perm)) deallocate(perm)
455  if (associated(iperm)) deallocate(iperm)
456  if (associated(alu)) deallocate(alu)
457  if (associated(d)) deallocate(d)
458  if (associated(al)) deallocate(al)
459  if (associated(au)) deallocate(au)
460  if (associated(indexl)) deallocate(indexl)
461  if (associated(indexu)) deallocate(indexu)
462  if (associated(iteml)) deallocate(iteml)
463  if (associated(itemu)) deallocate(itemu)
464  nullify(colorindex)
465  nullify(perm)
466  nullify(iperm)
467  nullify(alu)
468  nullify(d)
469  nullify(al)
470  nullify(au)
471  nullify(indexl)
472  nullify(indexu)
473  nullify(iteml)
474  nullify(itemu)
475  initialized = .false.
476  end subroutine hecmw_precond_ssor_44_clear
477 
478 end module hecmw_precond_ssor_44
real(kind=kreal) function, public hecmw_mat_get_omega(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_ncolor_in(hecMAT)
subroutine, public hecmw_matrix_reorder_values(N, NDOF, perm, iperm, indexL, indexU, itemL, itemU, AL, AU, D, indexLp, indexUp, itemLp, itemUp, ALp, AUp, Dp)
subroutine, public hecmw_matrix_reorder_renum_item(N, perm, indexXp, itemXp)
subroutine, public hecmw_matrix_reorder_profile(N, perm, iperm, indexL, indexU, itemL, itemU, indexLp, indexUp, itemLp, itemUp)
subroutine, public hecmw_precond_ssor_44_clear(hecMAT)
subroutine, public hecmw_precond_ssor_44_apply(ZP)
subroutine, public hecmw_precond_ssor_44_setup(hecMAT)
subroutine, public hecmw_tuning_fx_calc_sector_cache(N, NDOF, sectorCacheSize0, sectorCacheSize1)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
subroutine, public hecmw_matrix_ordering_rcm(N, indexL, itemL, indexU, itemU, perm, iperm)
subroutine, public hecmw_matrix_ordering_mc(N, indexL, itemL, indexU, itemU, perm_cur, ncolor_in, ncolor_out, COLORindex, perm, iperm)