FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_precond_DIAG_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_DIAG_66
9 !C***
10 !C
11 module hecmw_precond_diag_66
12  use hecmw_util
13 
14  private
15 
16  public:: hecmw_precond_diag_66_setup
17  public:: hecmw_precond_diag_66_apply
18  public:: hecmw_precond_diag_66_clear
19 
20  integer(kind=kint) :: N
21  real(kind=kreal), pointer :: alu(:) => null()
22 
23 contains
24 
25  subroutine hecmw_precond_diag_66_setup(hecMAT)
26  use hecmw_matrix_misc
27  implicit none
28  type(hecmwst_matrix), intent(in) :: hecmat
29  integer(kind=kint ) :: np
30  real (kind=kreal) :: sigma_diag
31  real(kind=kreal), pointer:: d(:)
32 
33  real (kind=kreal):: alutmp(6,6), pw(6)
34  integer(kind=kint ):: ii, i, j, k
35 
36  n = hecmat%N
37  np = hecmat%NP
38  d => hecmat%D
39  sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
40 
41  allocate(alu(36*np))
42  alu = 0.d0
43 
44  do ii= 1, n
45  alu(36*ii-35) = d(36*ii-35)
46  alu(36*ii-34) = d(36*ii-34)
47  alu(36*ii-33) = d(36*ii-33)
48  alu(36*ii-32) = d(36*ii-32)
49  alu(36*ii-31) = d(36*ii-31)
50  alu(36*ii-30) = d(36*ii-30)
51  alu(36*ii-29) = d(36*ii-29)
52  alu(36*ii-28) = d(36*ii-28)
53  alu(36*ii-27) = d(36*ii-27)
54  alu(36*ii-26) = d(36*ii-26)
55  alu(36*ii-25) = d(36*ii-25)
56  alu(36*ii-24) = d(36*ii-24)
57  alu(36*ii-23) = d(36*ii-23)
58  alu(36*ii-22) = d(36*ii-22)
59  alu(36*ii-21) = d(36*ii-21)
60  alu(36*ii-20) = d(36*ii-20)
61  alu(36*ii-19) = d(36*ii-19)
62  alu(36*ii-18) = d(36*ii-18)
63  alu(36*ii-17) = d(36*ii-17)
64  alu(36*ii-16) = d(36*ii-16)
65  alu(36*ii-15) = d(36*ii-15)
66  alu(36*ii-14) = d(36*ii-14)
67  alu(36*ii-13) = d(36*ii-13)
68  alu(36*ii-12) = d(36*ii-12)
69  alu(36*ii-11) = d(36*ii-11)
70  alu(36*ii-10) = d(36*ii-10)
71  alu(36*ii-9 ) = d(36*ii-9 )
72  alu(36*ii-8 ) = d(36*ii-8 )
73  alu(36*ii-7 ) = d(36*ii-7 )
74  alu(36*ii-6 ) = d(36*ii-6 )
75  alu(36*ii-5 ) = d(36*ii-5 )
76  alu(36*ii-4 ) = d(36*ii-4 )
77  alu(36*ii-3 ) = d(36*ii-3 )
78  alu(36*ii-2 ) = d(36*ii-2 )
79  alu(36*ii-1 ) = d(36*ii-1 )
80  alu(36*ii ) = d(36*ii )
81  enddo
82 
83  !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,SIGMA_DIAG)
84  !$omp do
85  do ii= 1, n
86  alutmp(1,1)= alu(36*ii-35) * sigma_diag
87  alutmp(1,2)= alu(36*ii-34)
88  alutmp(1,3)= alu(36*ii-33)
89  alutmp(1,4)= alu(36*ii-32)
90  alutmp(1,5)= alu(36*ii-31)
91  alutmp(1,6)= alu(36*ii-30)
92 
93  alutmp(2,1)= alu(36*ii-29)
94  alutmp(2,2)= alu(36*ii-28) * sigma_diag
95  alutmp(2,3)= alu(36*ii-27)
96  alutmp(2,4)= alu(36*ii-26)
97  alutmp(2,5)= alu(36*ii-25)
98  alutmp(2,6)= alu(36*ii-24)
99 
100  alutmp(3,1)= alu(36*ii-23)
101  alutmp(3,2)= alu(36*ii-22)
102  alutmp(3,3)= alu(36*ii-21) * sigma_diag
103  alutmp(3,4)= alu(36*ii-20)
104  alutmp(3,5)= alu(36*ii-19)
105  alutmp(3,6)= alu(36*ii-18)
106 
107  alutmp(4,1)= alu(36*ii-17)
108  alutmp(4,2)= alu(36*ii-16)
109  alutmp(4,3)= alu(36*ii-15)
110  alutmp(4,4)= alu(36*ii-14) * sigma_diag
111  alutmp(4,5)= alu(36*ii-13)
112  alutmp(4,6)= alu(36*ii-12)
113 
114  alutmp(5,1)= alu(36*ii-11)
115  alutmp(5,2)= alu(36*ii-10)
116  alutmp(5,3)= alu(36*ii-9 )
117  alutmp(5,4)= alu(36*ii-8 )
118  alutmp(5,5)= alu(36*ii-7 ) * sigma_diag
119  alutmp(5,6)= alu(36*ii-6 )
120 
121  alutmp(6,1)= alu(36*ii-5 )
122  alutmp(6,2)= alu(36*ii-4 )
123  alutmp(6,3)= alu(36*ii-3 )
124  alutmp(6,4)= alu(36*ii-2 )
125  alutmp(6,5)= alu(36*ii-1 )
126  alutmp(6,6)= alu(36*ii ) * sigma_diag
127 
128  do k= 1, 6
129  alutmp(k,k)= 1.d0/alutmp(k,k)
130  do i= k+1, 6
131  alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
132  do j= k+1, 6
133  pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
134  enddo
135  do j= k+1, 6
136  alutmp(i,j)= pw(j)
137  enddo
138  enddo
139  enddo
140 
141  alu(36*ii-35)= alutmp(1,1)
142  alu(36*ii-34)= alutmp(1,2)
143  alu(36*ii-33)= alutmp(1,3)
144  alu(36*ii-32)= alutmp(1,4)
145  alu(36*ii-31)= alutmp(1,5)
146  alu(36*ii-30)= alutmp(1,6)
147  alu(36*ii-29)= alutmp(2,1)
148  alu(36*ii-28)= alutmp(2,2)
149  alu(36*ii-27)= alutmp(2,3)
150  alu(36*ii-26)= alutmp(2,4)
151  alu(36*ii-25)= alutmp(2,5)
152  alu(36*ii-24)= alutmp(2,6)
153  alu(36*ii-23)= alutmp(3,1)
154  alu(36*ii-22)= alutmp(3,2)
155  alu(36*ii-21)= alutmp(3,3)
156  alu(36*ii-20)= alutmp(3,4)
157  alu(36*ii-19)= alutmp(3,5)
158  alu(36*ii-18)= alutmp(3,6)
159  alu(36*ii-17)= alutmp(4,1)
160  alu(36*ii-16)= alutmp(4,2)
161  alu(36*ii-15)= alutmp(4,3)
162  alu(36*ii-14)= alutmp(4,4)
163  alu(36*ii-13)= alutmp(4,5)
164  alu(36*ii-12)= alutmp(4,6)
165  alu(36*ii-11)= alutmp(5,1)
166  alu(36*ii-10)= alutmp(5,2)
167  alu(36*ii-9 )= alutmp(5,3)
168  alu(36*ii-8 )= alutmp(5,4)
169  alu(36*ii-7 )= alutmp(5,5)
170  alu(36*ii-6 )= alutmp(5,6)
171  alu(36*ii-5 )= alutmp(6,1)
172  alu(36*ii-4 )= alutmp(6,2)
173  alu(36*ii-3 )= alutmp(6,3)
174  alu(36*ii-2 )= alutmp(6,4)
175  alu(36*ii-1 )= alutmp(6,5)
176  alu(36*ii )= alutmp(6,6)
177  enddo
178  !$omp end do
179  !$omp end parallel
180 
181  end subroutine hecmw_precond_diag_66_setup
182 
183  subroutine hecmw_precond_diag_66_apply(WW)
184  implicit none
185  real(kind=kreal), intent(inout) :: ww(:)
186  integer(kind=kint) :: i
187  real(kind=kreal) :: x1, x2, x3, x4, x5, x6
188 
189  !C
190  !C== Block SCALING
191  !$omp parallel default(none),private(i,X1,X2,X3,X4,X5,X6),shared(N,WW,ALU)
192  !$omp do
193  do i= 1, n
194  x1= ww(6*i-5)
195  x2= ww(6*i-4)
196  x3= ww(6*i-3)
197  x4= ww(6*i-2)
198  x5= ww(6*i-1)
199  x6= ww(6*i )
200  x2= x2 -alu(36*i-29)*x1
201  x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
202  x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-15)*x3
203  x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9 )*x3 -alu(36*i-8)*x4
204  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
205  x6= alu(36*i )* x6
206  x5= alu(36*i-7 )*( x5 -alu(36*i-6 )*x6 )
207  x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
208  x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
209  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)
210  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)
211  ww(6*i-5)= x1
212  ww(6*i-4)= x2
213  ww(6*i-3)= x3
214  ww(6*i-2)= x4
215  ww(6*i-1)= x5
216  ww(6*i )= x6
217  enddo
218  !$omp end do
219  !$omp end parallel
220 
221  end subroutine hecmw_precond_diag_66_apply
222 
223  subroutine hecmw_precond_diag_66_clear()
224  implicit none
225  if (associated(alu)) deallocate(alu)
226  nullify(alu)
227  end subroutine hecmw_precond_diag_66_clear
228 
229 end module hecmw_precond_diag_66
hecmw_precond_diag_66::hecmw_precond_diag_66_clear
subroutine, public hecmw_precond_diag_66_clear()
Definition: hecmw_precond_DIAG_66.f90:224
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_matrix_misc::hecmw_mat_get_sigma_diag
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
Definition: hecmw_matrix_misc.f90:694
hecmw_precond_diag_66
Definition: hecmw_precond_DIAG_66.f90:11
hecmw_precond_diag_66::hecmw_precond_diag_66_setup
subroutine, public hecmw_precond_diag_66_setup(hecMAT)
Definition: hecmw_precond_DIAG_66.f90:26
hecmw_matrix_misc
Definition: hecmw_matrix_misc.f90:6
hecmw_precond_diag_66::hecmw_precond_diag_66_apply
subroutine, public hecmw_precond_diag_66_apply(WW)
Definition: hecmw_precond_DIAG_66.f90:184
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444