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