FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_precond_DIAG_33.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_33
9 !C***
10 !C
12  use hecmw_util
13 
14  private
15 
19 
20  integer(kind=kint) :: N
21  real(kind=kreal), pointer :: alu(:) => null()
22 
23  logical, save :: INITIALIZED = .false.
24 
25 contains
26 
27  subroutine hecmw_precond_diag_33_setup(hecMAT)
29  implicit none
30  type(hecmwst_matrix), intent(inout) :: hecmat
31  integer(kind=kint ) :: np
32  real (kind=kreal) :: sigma_diag
33  real(kind=kreal), pointer:: d(:)
34 
35  real (kind=kreal):: alutmp(3,3), pw(3)
36  integer(kind=kint ):: ii, i, j, k
37 
38  if (initialized) then
39  if (hecmat%Iarray(98) == 0 .and. hecmat%Iarray(97) == 0) return
41  endif
42 
43  n = hecmat%N
44  np = hecmat%NP
45  d => hecmat%D
46  sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
47 
48  allocate(alu(9*np))
49  alu = 0.d0
50 
51  do ii= 1, n
52  alu(9*ii-8) = d(9*ii-8)
53  alu(9*ii-7) = d(9*ii-7)
54  alu(9*ii-6) = d(9*ii-6)
55  alu(9*ii-5) = d(9*ii-5)
56  alu(9*ii-4) = d(9*ii-4)
57  alu(9*ii-3) = d(9*ii-3)
58  alu(9*ii-2) = d(9*ii-2)
59  alu(9*ii-1) = d(9*ii-1)
60  alu(9*ii ) = d(9*ii )
61  enddo
62 
63  if (hecmat%cmat%n_val.gt.0) then
64  do k= 1, hecmat%cmat%n_val
65  if (hecmat%cmat%pair(k)%i.ne.hecmat%cmat%pair(k)%j) cycle
66  ii = hecmat%cmat%pair(k)%i
67  alu(9*ii-8) = alu(9*ii-8) + hecmat%cmat%pair(k)%val(1, 1)
68  alu(9*ii-7) = alu(9*ii-7) + hecmat%cmat%pair(k)%val(1, 2)
69  alu(9*ii-6) = alu(9*ii-6) + hecmat%cmat%pair(k)%val(1, 3)
70  alu(9*ii-5) = alu(9*ii-5) + hecmat%cmat%pair(k)%val(2, 1)
71  alu(9*ii-4) = alu(9*ii-4) + hecmat%cmat%pair(k)%val(2, 2)
72  alu(9*ii-3) = alu(9*ii-3) + hecmat%cmat%pair(k)%val(2, 3)
73  alu(9*ii-2) = alu(9*ii-2) + hecmat%cmat%pair(k)%val(3, 1)
74  alu(9*ii-1) = alu(9*ii-1) + hecmat%cmat%pair(k)%val(3, 2)
75  alu(9*ii ) = alu(9*ii ) + hecmat%cmat%pair(k)%val(3, 3)
76  enddo
77 
78  !call hecmw_cmat_LU( hecMAT )
79 
80  endif
81 
82  !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,SIGMA_DIAG)
83  !$omp do
84  do ii= 1, n
85  alutmp(1,1)= alu(9*ii-8) * sigma_diag
86  alutmp(1,2)= alu(9*ii-7)
87  alutmp(1,3)= alu(9*ii-6)
88  alutmp(2,1)= alu(9*ii-5)
89  alutmp(2,2)= alu(9*ii-4) * sigma_diag
90  alutmp(2,3)= alu(9*ii-3)
91  alutmp(3,1)= alu(9*ii-2)
92  alutmp(3,2)= alu(9*ii-1)
93  alutmp(3,3)= alu(9*ii ) * sigma_diag
94  do k= 1, 3
95  alutmp(k,k)= 1.d0/alutmp(k,k)
96  do i= k+1, 3
97  alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
98  do j= k+1, 3
99  pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
100  enddo
101  do j= k+1, 3
102  alutmp(i,j)= pw(j)
103  enddo
104  enddo
105  enddo
106  alu(9*ii-8)= alutmp(1,1)
107  alu(9*ii-7)= alutmp(1,2)
108  alu(9*ii-6)= alutmp(1,3)
109  alu(9*ii-5)= alutmp(2,1)
110  alu(9*ii-4)= alutmp(2,2)
111  alu(9*ii-3)= alutmp(2,3)
112  alu(9*ii-2)= alutmp(3,1)
113  alu(9*ii-1)= alutmp(3,2)
114  alu(9*ii )= alutmp(3,3)
115  enddo
116  !$omp end do
117  !$omp end parallel
118 
119  initialized = .true.
120  hecmat%Iarray(98) = 0 ! symbolic setup done
121  hecmat%Iarray(97) = 0 ! numerical setup done
122 
123  end subroutine hecmw_precond_diag_33_setup
124 
126  implicit none
127  real(kind=kreal), intent(inout) :: ww(:)
128  integer(kind=kint) :: i
129  real(kind=kreal) :: x1, x2, x3
130 
131  !C
132  !C== Block SCALING
133 
134  !$omp parallel default(none),private(i,X1,X2,X3),shared(N,WW,ALU)
135  !$omp do
136  do i= 1, n
137  x1= ww(3*i-2)
138  x2= ww(3*i-1)
139  x3= ww(3*i )
140  x2= x2 - alu(9*i-5)*x1
141  x3= x3 - alu(9*i-2)*x1 - alu(9*i-1)*x2
142  x3= alu(9*i )* x3
143  x2= alu(9*i-4)*( x2 - alu(9*i-3)*x3 )
144  x1= alu(9*i-8)*( x1 - alu(9*i-6)*x3 - alu(9*i-7)*x2)
145  ww(3*i-2)= x1
146  ww(3*i-1)= x2
147  ww(3*i )= x3
148  enddo
149  !$omp end do
150  !$omp end parallel
151 
152  end subroutine hecmw_precond_diag_33_apply
153 
155  implicit none
156  if (associated(alu)) deallocate(alu)
157  nullify(alu)
158  initialized = .false.
159  end subroutine hecmw_precond_diag_33_clear
160 
161 end module hecmw_precond_diag_33
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
subroutine, public hecmw_precond_diag_33_apply(WW)
subroutine, public hecmw_precond_diag_33_clear()
subroutine, public hecmw_precond_diag_33_setup(hecMAT)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal