FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
hecmw_precond_SSOR_66.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_66
9 !C***
10 !C
11 module hecmw_precond_ssor_66
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_66_setup
24  public:: hecmw_precond_ssor_66_apply
25  public:: hecmw_precond_ssor_66_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 contains
45 
46  subroutine hecmw_precond_ssor_66_setup(hecMAT)
47  implicit none
48  type(hecmwst_matrix), intent(in) :: hecmat
49  integer(kind=kint ) :: npl, npu
50  integer(kind=kint ) :: ncolor_in
51  real (kind=kreal) :: sigma_diag
52  real (kind=kreal) :: alutmp(6,6), pw(6)
53  integer(kind=kint ) :: ii, i, j, k
54  integer(kind=kint ) :: nthreads = 1
55  integer(kind=kint ), allocatable :: perm_tmp(:)
56  !real (kind=kreal) :: t0
57 
58  !t0 = hecmw_Wtime()
59  !write(*,*) 'DEBUG: SSOR setup start', hecmw_Wtime()-t0
60 
61 #ifndef _OPENACC
62  !$ nthreads = omp_get_max_threads()
63 #endif
64 
65  n = hecmat%N
66  ! N = hecMAT%NP
67  ncolor_in = hecmw_mat_get_ncolor_in(hecmat)
68  sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
69 
70 #ifdef _OPENACC
71  allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
72  call hecmw_matrix_ordering_rcm(n, hecmat%indexL, hecmat%itemL, &
73  hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
74  !write(*,*) 'DEBUG: RCM ordering done', hecmw_Wtime()-t0
75  call hecmw_matrix_ordering_mc(n, hecmat%indexL, hecmat%itemL, &
76  hecmat%indexU, hecmat%itemU, perm_tmp, &
77  ncolor_in, ncolor, colorindex, perm, iperm)
78  !write(*,*) 'DEBUG: MC ordering done', hecmw_Wtime()-t0
79  deallocate(perm_tmp)
80 
81  !call write_debug_info
82 
83  npl = hecmat%indexL(n)
84  npu = hecmat%indexU(n)
85  allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
86  call hecmw_matrix_reorder_profile(n, perm, iperm, &
87  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
88  indexl, indexu, iteml, itemu)
89  !write(*,*) 'DEBUG: reordering profile done', hecmw_Wtime()-t0
90 
91  call check_ordering
92 
93  allocate(d(36*n), al(36*npl), au(36*npu))
94  call hecmw_matrix_reorder_values(n, 6, perm, iperm, &
95  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
96  hecmat%AL, hecmat%AU, hecmat%D, &
97  indexl, indexu, iteml, itemu, al, au, d)
98  !write(*,*) 'DEBUG: reordering values done', hecmw_Wtime()-t0
99 
100  call hecmw_matrix_reorder_renum_item(n, perm, indexl, iteml)
101  call hecmw_matrix_reorder_renum_item(n, perm, indexu, itemu)
102 #else
103  if (nthreads == 1) then
104  ncolor = 1
105  allocate(colorindex(0:1), perm(n), iperm(n))
106  colorindex(0) = 0
107  colorindex(1) = n
108  do i=1,n
109  perm(i) = i
110  iperm(i) = i
111  end do
112 
113  d => hecmat%D
114  al => hecmat%AL
115  au => hecmat%AU
116  indexl => hecmat%indexL
117  indexu => hecmat%indexU
118  iteml => hecmat%itemL
119  itemu => hecmat%itemU
120  else
121  allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
122  call hecmw_matrix_ordering_rcm(n, hecmat%indexL, hecmat%itemL, &
123  hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
124  !write(*,*) 'DEBUG: RCM ordering done', hecmw_Wtime()-t0
125  call hecmw_matrix_ordering_mc(n, hecmat%indexL, hecmat%itemL, &
126  hecmat%indexU, hecmat%itemU, perm_tmp, &
127  ncolor_in, ncolor, colorindex, perm, iperm)
128  !write(*,*) 'DEBUG: MC ordering done', hecmw_Wtime()-t0
129  deallocate(perm_tmp)
130 
131  !call write_debug_info
132 
133  npl = hecmat%indexL(n)
134  npu = hecmat%indexU(n)
135  allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
136  call hecmw_matrix_reorder_profile(n, perm, iperm, &
137  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
138  indexl, indexu, iteml, itemu)
139  !write(*,*) 'DEBUG: reordering profile done', hecmw_Wtime()-t0
140 
141  call check_ordering
142 
143  allocate(d(36*n), al(36*npl), au(36*npu))
144  call hecmw_matrix_reorder_values(n, 6, perm, iperm, &
145  hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
146  hecmat%AL, hecmat%AU, hecmat%D, &
147  indexl, indexu, iteml, itemu, al, au, d)
148  !write(*,*) 'DEBUG: reordering values done', hecmw_Wtime()-t0
149 
150  call hecmw_matrix_reorder_renum_item(n, perm, indexl, iteml)
151  call hecmw_matrix_reorder_renum_item(n, perm, indexu, itemu)
152  end if
153 #endif
154 
155  allocate(alu(36*n))
156  alu = 0.d0
157 
158  do ii= 1, 36*n
159  alu(ii) = d(ii)
160  enddo
161 
162 #ifdef _OPENACC
163  !$acc kernels
164  !$acc loop independent private(ALUtmp,PW)
165 #else
166  !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,SIGMA_DIAG)
167  !$omp do
168 #endif
169  do ii= 1, n
170  alutmp(1,1)= alu(36*ii-35) * sigma_diag
171  alutmp(1,2)= alu(36*ii-34)
172  alutmp(1,3)= alu(36*ii-33)
173  alutmp(1,4)= alu(36*ii-32)
174  alutmp(1,5)= alu(36*ii-31)
175  alutmp(1,6)= alu(36*ii-30)
176 
177  alutmp(2,1)= alu(36*ii-29)
178  alutmp(2,2)= alu(36*ii-28) * sigma_diag
179  alutmp(2,3)= alu(36*ii-27)
180  alutmp(2,4)= alu(36*ii-26)
181  alutmp(2,5)= alu(36*ii-25)
182  alutmp(2,6)= alu(36*ii-24)
183 
184  alutmp(3,1)= alu(36*ii-23)
185  alutmp(3,2)= alu(36*ii-22)
186  alutmp(3,3)= alu(36*ii-21) * sigma_diag
187  alutmp(3,4)= alu(36*ii-20)
188  alutmp(3,5)= alu(36*ii-19)
189  alutmp(3,6)= alu(36*ii-18)
190 
191  alutmp(4,1)= alu(36*ii-17)
192  alutmp(4,2)= alu(36*ii-16)
193  alutmp(4,3)= alu(36*ii-15)
194  alutmp(4,4)= alu(36*ii-14) * sigma_diag
195  alutmp(4,5)= alu(36*ii-13)
196  alutmp(4,6)= alu(36*ii-12)
197 
198  alutmp(5,1)= alu(36*ii-11)
199  alutmp(5,2)= alu(36*ii-10)
200  alutmp(5,3)= alu(36*ii-9 )
201  alutmp(5,4)= alu(36*ii-8 )
202  alutmp(5,5)= alu(36*ii-7 ) * sigma_diag
203  alutmp(5,6)= alu(36*ii-6 )
204 
205  alutmp(6,1)= alu(36*ii-5 )
206  alutmp(6,2)= alu(36*ii-4 )
207  alutmp(6,3)= alu(36*ii-3 )
208  alutmp(6,4)= alu(36*ii-2 )
209  alutmp(6,5)= alu(36*ii-1 )
210  alutmp(6,6)= alu(36*ii ) * sigma_diag
211 
212  do k= 1, 6
213  alutmp(k,k)= 1.d0/alutmp(k,k)
214  do i= k+1, 6
215  alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
216  do j= k+1, 6
217  pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
218  enddo
219  do j= k+1, 6
220  alutmp(i,j)= pw(j)
221  enddo
222  enddo
223  enddo
224 
225  alu(36*ii-35)= alutmp(1,1)
226  alu(36*ii-34)= alutmp(1,2)
227  alu(36*ii-33)= alutmp(1,3)
228  alu(36*ii-32)= alutmp(1,4)
229  alu(36*ii-31)= alutmp(1,5)
230  alu(36*ii-30)= alutmp(1,6)
231  alu(36*ii-29)= alutmp(2,1)
232  alu(36*ii-28)= alutmp(2,2)
233  alu(36*ii-27)= alutmp(2,3)
234  alu(36*ii-26)= alutmp(2,4)
235  alu(36*ii-25)= alutmp(2,5)
236  alu(36*ii-24)= alutmp(2,6)
237  alu(36*ii-23)= alutmp(3,1)
238  alu(36*ii-22)= alutmp(3,2)
239  alu(36*ii-21)= alutmp(3,3)
240  alu(36*ii-20)= alutmp(3,4)
241  alu(36*ii-19)= alutmp(3,5)
242  alu(36*ii-18)= alutmp(3,6)
243  alu(36*ii-17)= alutmp(4,1)
244  alu(36*ii-16)= alutmp(4,2)
245  alu(36*ii-15)= alutmp(4,3)
246  alu(36*ii-14)= alutmp(4,4)
247  alu(36*ii-13)= alutmp(4,5)
248  alu(36*ii-12)= alutmp(4,6)
249  alu(36*ii-11)= alutmp(5,1)
250  alu(36*ii-10)= alutmp(5,2)
251  alu(36*ii-9 )= alutmp(5,3)
252  alu(36*ii-8 )= alutmp(5,4)
253  alu(36*ii-7 )= alutmp(5,5)
254  alu(36*ii-6 )= alutmp(5,6)
255  alu(36*ii-5 )= alutmp(6,1)
256  alu(36*ii-4 )= alutmp(6,2)
257  alu(36*ii-3 )= alutmp(6,3)
258  alu(36*ii-2 )= alutmp(6,4)
259  alu(36*ii-1 )= alutmp(6,5)
260  alu(36*ii )= alutmp(6,6)
261  enddo
262 #ifdef _OPENACC
263  !$acc end kernels
264 #else
265  !$omp end do
266  !$omp end parallel
267 #endif
268 
269  isfirst = .true.
270 
271  !write(*,*) 'DEBUG: SSOR setup done', hecmw_Wtime()-t0
272 
273  end subroutine hecmw_precond_ssor_66_setup
274 
275  subroutine hecmw_precond_ssor_66_apply(ZP)
276  use hecmw_tuning_fx
277  implicit none
278  real(kind=kreal), intent(inout) :: zp(:)
279  integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
280  real(kind=kreal) :: x1, x2, x3, x4, x5, x6
281  real(kind=kreal) :: sw1, sw2, sw3, sw4, sw5, sw6
282 
283  ! added for tuning >>>
284  integer(kind=kint), parameter :: numofblockperthread = 100
285  integer(kind=kint), save :: numofthread = 1, numofblock
286  integer(kind=kint), save, allocatable :: ictoblockindex(:)
287  integer(kind=kint), save, allocatable :: blockindextocolorindex(:)
288  integer(kind=kint), save :: sectorcachesize0, sectorcachesize1
289  integer(kind=kint) :: blockindex, elementcount, numofelement, ii
290  real(kind=kreal) :: numofelementperblock
291  integer(kind=kint) :: my_rank
292 
293 #ifndef _OPENACC
294  if (isfirst) then
295  !$ numOfThread = omp_get_max_threads()
296  numofblock = numofthread * numofblockperthread
297  if (allocated(ictoblockindex)) deallocate(ictoblockindex)
298  if (allocated(blockindextocolorindex)) deallocate(blockindextocolorindex)
299  allocate (ictoblockindex(0:ncolor), &
300  blockindextocolorindex(0:numofblock + ncolor))
301  numofelement = n + indexl(n) + indexu(n)
302  numofelementperblock = dble(numofelement) / numofblock
303  blockindex = 0
304  ictoblockindex = -1
305  ictoblockindex(0) = 0
306  blockindextocolorindex = -1
307  blockindextocolorindex(0) = 0
308  my_rank = hecmw_comm_get_rank()
309  ! write(9000+my_rank,*) &
310  ! '# numOfElementPerBlock =', numOfElementPerBlock
311  ! write(9000+my_rank,*) &
312  ! '# ic, blockIndex, colorIndex, elementCount'
313  do ic = 1, ncolor
314  elementcount = 0
315  ii = 1
316  do i = colorindex(ic-1)+1, colorindex(ic)
317  elementcount = elementcount + 1
318  elementcount = elementcount + (indexl(i) - indexl(i-1))
319  elementcount = elementcount + (indexu(i) - indexu(i-1))
320  if (elementcount > ii * numofelementperblock &
321  .or. i == colorindex(ic)) then
322  ii = ii + 1
323  blockindex = blockindex + 1
324  blockindextocolorindex(blockindex) = i
325  ! write(9000+my_rank,*) ic, blockIndex, &
326  ! blockIndexToColorIndex(blockIndex), elementCount
327  endif
328  enddo
329  ictoblockindex(ic) = blockindex
330  enddo
331  numofblock = blockindex
332 
333  call hecmw_tuning_fx_calc_sector_cache( n, 3, &
334  sectorcachesize0, sectorcachesize1 )
335 
336  isfirst = .false.
337  endif
338 #endif
339  ! <<< added for tuning
340 
341 #ifndef _OPENACC
342  !call start_collection("loopInPrecond66")
343 
344  !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
345  !OCL CACHE_SUBSECTOR_ASSIGN(ZP)
346 
347  !$omp parallel default(none) &
348  !$omp&shared(NColor,indexL,itemL,indexU,itemU,AL,AU,D,ALU,perm,&
349  !$omp& ZP,icToBlockIndex,blockIndexToColorIndex) &
350  !$omp&private(SW1,SW2,SW3,SW4,SW5,SW6,X1,X2,X3,X4,X5,X6,ic,i,iold,isL,ieL,isU,ieU,j,k,blockIndex)
351 #endif
352 
353  !C-- FORWARD
354  do ic=1,ncolor
355 #ifdef _OPENACC
356  !$acc kernels
357  !$acc loop independent
358  do i = colorindex(ic-1)+1, colorindex(ic)
359 #else
360  !$omp do schedule (static, 1)
361  do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
362  do i = blockindextocolorindex(blockindex-1)+1, &
363  blockindextocolorindex(blockindex)
364 #endif
365  ! do i = startPos(threadNum, ic), endPos(threadNum, ic)
366  iold = perm(i)
367  sw1= zp(6*iold-5)
368  sw2= zp(6*iold-4)
369  sw3= zp(6*iold-3)
370  sw4= zp(6*iold-2)
371  sw5= zp(6*iold-1)
372  sw6= zp(6*iold )
373  isl= indexl(i-1)+1
374  iel= indexl(i)
375  do j= isl, iel
376  !k= perm(itemL(j))
377  k= iteml(j)
378  x1= zp(6*k-5)
379  x2= zp(6*k-4)
380  x3= zp(6*k-3)
381  x4= zp(6*k-2)
382  x5= zp(6*k-1)
383  x6= zp(6*k )
384  sw1= sw1 -al(36*j-35)*x1 -al(36*j-34)*x2 -al(36*j-33)*x3 -al(36*j-32)*x4 -al(36*j-31)*x5 -al(36*j-30)*x6
385  sw2= sw2 -al(36*j-29)*x1 -al(36*j-28)*x2 -al(36*j-27)*x3 -al(36*j-26)*x4 -al(36*j-25)*x5 -al(36*j-24)*x6
386  sw3= sw3 -al(36*j-23)*x1 -al(36*j-22)*x2 -al(36*j-21)*x3 -al(36*j-20)*x4 -al(36*j-19)*x5 -al(36*j-18)*x6
387  sw4= sw4 -al(36*j-17)*x1 -al(36*j-16)*x2 -al(36*j-15)*x3 -al(36*j-14)*x4 -al(36*j-13)*x5 -al(36*j-12)*x6
388  sw5= sw5 -al(36*j-11)*x1 -al(36*j-10)*x2 -al(36*j-9 )*x3 -al(36*j-8 )*x4 -al(36*j-7 )*x5 -al(36*j-6 )*x6
389  sw6= sw6 -al(36*j-5 )*x1 -al(36*j-4 )*x2 -al(36*j-3 )*x3 -al(36*j-2 )*x4 -al(36*j-1 )*x5 -al(36*j )*x6
390  enddo ! j
391 
392  x1= sw1
393  x2= sw2
394  x3= sw3
395  x4= sw4
396  x5= sw5
397  x6= sw6
398  x2= x2 -alu(36*i-29)*x1
399  x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
400  x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-5)*x3
401  x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9)*x3 -alu(36*i-8)*x4
402  x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3)*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
403  x6= alu(36*i )* x6
404  x5= alu(36*i-7 )*( x5 -alu(36*i-6)*x6 )
405  x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
406  x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
407  x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
408  x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
409  zp(6*iold-5)= x1
410  zp(6*iold-4)= x2
411  zp(6*iold-3)= x3
412  zp(6*iold-2)= x4
413  zp(6*iold-1)= x5
414  zp(6*iold )= x6
415 #ifdef _OPENACC
416  enddo
417  !$acc end kernels
418 #else
419  enddo ! i
420  enddo ! blockIndex
421  !$omp end do
422 #endif
423  enddo ! ic
424 
425  !C-- BACKWARD
426  do ic=ncolor, 1, -1
427 #ifdef _OPENACC
428  !$acc kernels
429  !$acc loop independent
430  do i = colorindex(ic-1)+1, colorindex(ic)
431 #else
432  !$omp do schedule (static, 1)
433  do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
434  do i = blockindextocolorindex(blockindex), &
435  blockindextocolorindex(blockindex-1)+1, -1
436 #endif
437  ! do blockIndex = icToBlockIndex(ic-1)+1, icToBlockIndex(ic)
438  ! do i = blockIndexToColorIndex(blockIndex-1)+1, &
439  ! blockIndexToColorIndex(blockIndex)
440  ! do i = endPos(threadNum, ic), startPos(threadNum, ic), -1
441  sw1= 0.d0
442  sw2= 0.d0
443  sw3= 0.d0
444  sw4= 0.d0
445  sw5= 0.d0
446  sw6= 0.d0
447  isu= indexu(i-1) + 1
448  ieu= indexu(i)
449  do j= ieu, isu, -1
450  !k= perm(itemU(j))
451  k= itemu(j)
452  x1= zp(6*k-5)
453  x2= zp(6*k-4)
454  x3= zp(6*k-3)
455  x4= zp(6*k-2)
456  x5= zp(6*k-1)
457  x6= zp(6*k )
458  sw1= sw1 +au(36*j-35)*x1 +au(36*j-34)*x2 +au(36*j-33)*x3 +au(36*j-32)*x4 +au(36*j-31)*x5 +au(36*j-30)*x6
459  sw2= sw2 +au(36*j-29)*x1 +au(36*j-28)*x2 +au(36*j-27)*x3 +au(36*j-26)*x4 +au(36*j-25)*x5 +au(36*j-24)*x6
460  sw3= sw3 +au(36*j-23)*x1 +au(36*j-22)*x2 +au(36*j-21)*x3 +au(36*j-20)*x4 +au(36*j-19)*x5 +au(36*j-18)*x6
461  sw4= sw4 +au(36*j-17)*x1 +au(36*j-16)*x2 +au(36*j-15)*x3 +au(36*j-14)*x4 +au(36*j-13)*x5 +au(36*j-12)*x6
462  sw5= sw5 +au(36*j-11)*x1 +au(36*j-10)*x2 +au(36*j-9 )*x3 +au(36*j-8 )*x4 +au(36*j-7 )*x5 +au(36*j-6 )*x6
463  sw6= sw6 +au(36*j-5 )*x1 +au(36*j-4 )*x2 +au(36*j-3 )*x3 +au(36*j-2 )*x4 +au(36*j-1 )*x5 +au(36*j )*x6
464  enddo ! j
465 
466  x1= sw1
467  x2= sw2
468  x3= sw3
469  x4= sw4
470  x5= sw5
471  x6= sw6
472  x2= x2 -alu(36*i-29)*x1
473  x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
474  x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-5)*x3
475  x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9)*x3 -alu(36*i-8)*x4
476  x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3)*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
477  x6= alu(36*i )* x6
478  x5= alu(36*i-7 )*( x5 -alu(36*i-6)*x6 )
479  x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
480  x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
481  x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
482  x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
483  iold = perm(i)
484  zp(6*iold-5)= zp(6*iold-5) -x1
485  zp(6*iold-4)= zp(6*iold-4) -x2
486  zp(6*iold-3)= zp(6*iold-3) -x3
487  zp(6*iold-2)= zp(6*iold-2) -x4
488  zp(6*iold-1)= zp(6*iold-1) -x5
489  zp(6*iold )= zp(6*iold ) -x6
490 #ifdef _OPENACC
491  enddo
492  !$acc end kernels
493 #else
494  enddo ! i
495  enddo ! blockIndex
496  !$omp end do
497 #endif
498  enddo ! ic
499 #ifndef _OPENACC
500  !$omp end parallel
501 
502  !OCL END_CACHE_SUBSECTOR
503  !OCL END_CACHE_SECTOR_SIZE
504 
505  !call stop_collection("loopInPrecond66")
506 #endif
507 
508  end subroutine hecmw_precond_ssor_66_apply
509 
510  subroutine hecmw_precond_ssor_66_clear(hecMAT)
511  implicit none
512  type(hecmwst_matrix), intent(inout) :: hecmat
513  integer(kind=kint ) :: nthreads = 1
514 #ifndef _OPENACC
515  !$ nthreads = omp_get_max_threads()
516 #endif
517  if (associated(colorindex)) deallocate(colorindex)
518  if (associated(perm)) deallocate(perm)
519  if (associated(iperm)) deallocate(iperm)
520  if (associated(alu)) deallocate(alu)
521  if (nthreads >= 1) then
522  if (associated(d)) deallocate(d)
523  if (associated(al)) deallocate(al)
524  if (associated(au)) deallocate(au)
525  if (associated(indexl)) deallocate(indexl)
526  if (associated(indexu)) deallocate(indexu)
527  if (associated(iteml)) deallocate(iteml)
528  if (associated(itemu)) deallocate(itemu)
529  end if
530  nullify(colorindex)
531  nullify(perm)
532  nullify(iperm)
533  nullify(alu)
534  nullify(d)
535  nullify(al)
536  nullify(au)
537  nullify(indexl)
538  nullify(indexu)
539  nullify(iteml)
540  nullify(itemu)
541  end subroutine hecmw_precond_ssor_66_clear
542 
543  subroutine write_debug_info
544  implicit none
545  integer(kind=kint) :: my_rank, ic, in
546  my_rank = hecmw_comm_get_rank()
547  !--------------------> debug: shizawa
548  if (my_rank.eq.0) then
549  write(*,*) 'DEBUG: Output fort.19000+myrank and fort.29000+myrank for coloring information'
550  endif
551  write(19000+my_rank,'(a)') '#NCOLORTot'
552  write(19000+my_rank,*) ncolor
553  write(19000+my_rank,'(a)') '#ic COLORindex(ic-1)+1 COLORindex(ic)'
554  do ic=1,ncolor
555  write(19000+my_rank,*) ic, colorindex(ic-1)+1,colorindex(ic)
556  enddo ! ic
557  write(29000+my_rank,'(a)') '#n_node'
558  write(29000+my_rank,*) n
559  write(29000+my_rank,'(a)') '#in OLDtoNEW(in) NEWtoOLD(in)'
560  do in=1,n
561  write(29000+my_rank,*) in, iperm(in), perm(in)
562  if (perm(iperm(in)) .ne. in) then
563  write(29000+my_rank,*) '** WARNING **: NEWtoOLD and OLDtoNEW: ',in
564  endif
565  enddo
566  end subroutine write_debug_info
567 
568  subroutine check_ordering
569  implicit none
570  integer(kind=kint) :: ic, i, j, k
571  integer(kind=kint), allocatable :: iicolor(:)
572  ! check color dependence of neighbouring nodes
573  if (ncolor.gt.1) then
574  allocate(iicolor(n))
575  do ic=1,ncolor
576  do i= colorindex(ic-1)+1, colorindex(ic)
577  iicolor(i) = ic
578  enddo ! i
579  enddo ! ic
580  ! FORWARD: L-part
581  do ic=1,ncolor
582  do i= colorindex(ic-1)+1, colorindex(ic)
583  do j= indexl(i-1)+1, indexl(i)
584  k= iteml(j)
585  if (iicolor(i).eq.iicolor(k)) then
586  write(*,*) .eq.'** ERROR **: L-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
587  endif
588  enddo ! j
589  enddo ! i
590  enddo ! ic
591  ! BACKWARD: U-part
592  do ic=ncolor, 1, -1
593  do i= colorindex(ic), colorindex(ic-1)+1, -1
594  do j= indexu(i-1)+1, indexu(i)
595  k= itemu(j)
596  if (iicolor(i).eq.iicolor(k)) then
597  write(*,*) .eq.'** ERROR **: U-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
598  endif
599  enddo ! j
600  enddo ! i
601  enddo ! ic
602  deallocate(iicolor)
603  endif ! if (NColor.gt.1)
604  !--------------------< debug: shizawa
605  end subroutine check_ordering
606 
607 end module hecmw_precond_ssor_66
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:700
hecmw_matrix_misc::hecmw_mat_get_ncolor_in
integer(kind=kint) function, public hecmw_mat_get_ncolor_in(hecMAT)
Definition: hecmw_matrix_misc.F90:537
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:71
hecmw_matrix_reorder::hecmw_matrix_reorder_renum_item
subroutine, public hecmw_matrix_reorder_renum_item(N, perm, indexXp, itemXp)
Definition: hecmw_matrix_reorder.f90:144
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:22
hecmw_precond_ssor_66
Definition: hecmw_precond_SSOR_66.F90:11
hecmw_precond_ssor_66::hecmw_precond_ssor_66_setup
subroutine, public hecmw_precond_ssor_66_setup(hecMAT)
Definition: hecmw_precond_SSOR_66.F90:47
hecmw_precond_ssor_66::hecmw_precond_ssor_66_clear
subroutine, public hecmw_precond_ssor_66_clear(hecMAT)
Definition: hecmw_precond_SSOR_66.F90:511
hecmw_precond_ssor_66::hecmw_precond_ssor_66_apply
subroutine, public hecmw_precond_ssor_66_apply(ZP)
Definition: hecmw_precond_SSOR_66.F90:276
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:590
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:17
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:449