FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_precond_DIAG_22.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_22
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_22_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(2,2), pw(2)
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(4*np))
49  alu = 0.d0
50 
51  do ii= 1, n
52  alu(4*ii-3) = d(4*ii-3)
53  alu(4*ii-2) = d(4*ii-2)
54  alu(4*ii-1) = d(4*ii-1)
55  alu(4*ii-0) = d(4*ii-0)
56  enddo
57 
58  if (hecmat%cmat%n_val.gt.0) then
59  do k= 1, hecmat%cmat%n_val
60  if (hecmat%cmat%pair(k)%i.ne.hecmat%cmat%pair(k)%j) cycle
61  ii = hecmat%cmat%pair(k)%i
62  alu(4*ii-3) = alu(4*ii-3) + hecmat%cmat%pair(k)%val(1, 1)
63  alu(4*ii-2) = alu(4*ii-2) + hecmat%cmat%pair(k)%val(1, 2)
64  alu(4*ii-1) = alu(4*ii-1) + hecmat%cmat%pair(k)%val(2, 1)
65  alu(4*ii-0) = alu(4*ii-0) + hecmat%cmat%pair(k)%val(2, 2)
66  enddo
67 
68  !call hecmw_cmat_LU( hecMAT )
69 
70  endif
71 
72  !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,SIGMA_DIAG)
73  !$omp do
74  do ii= 1, n
75  alutmp(1,1)= alu(4*ii-3) * sigma_diag
76  alutmp(1,2)= alu(4*ii-2)
77  alutmp(2,1)= alu(4*ii-1)
78  alutmp(2,2)= alu(4*ii-0) * sigma_diag
79 
80  do k= 1, 2
81  alutmp(k,k)= 1.d0/alutmp(k,k)
82  do i= k+1, 2
83  alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
84  do j= k+1, 2
85  pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
86  enddo
87  do j= k+1, 2
88  alutmp(i,j)= pw(j)
89  enddo
90  enddo
91  enddo
92  alu(4*ii-3)= alutmp(1,1)
93  alu(4*ii-2)= alutmp(1,2)
94  alu(4*ii-1)= alutmp(2,1)
95  alu(4*ii-0)= alutmp(2,2)
96  enddo
97  !$omp end do
98  !$omp end parallel
99 
100  initialized = .true.
101  hecmat%Iarray(98) = 0 ! symbolic setup done
102  hecmat%Iarray(97) = 0 ! numerical setup done
103 
104  end subroutine hecmw_precond_diag_22_setup
105 
107  implicit none
108  real(kind=kreal), intent(inout) :: ww(:)
109  integer(kind=kint) :: i
110  real(kind=kreal) :: x1, x2
111 
112  !C
113  !C== Block SCALING
114 
115  !$omp parallel default(none),private(i,X1,X2),shared(N,WW,ALU)
116  !$omp do
117  do i= 1, n
118  x1= ww(2*i-1)
119  x2= ww(2*i-0)
120  x2= x2 - alu(4*i-1)*x1
121  x2= alu(4*i )* x2
122  x1= alu(4*i-3)*( x1 - alu(4*i-2)*x2 )
123  ww(2*i-1)= x1
124  ww(2*i-0)= x2
125  enddo
126  !$omp end do
127  !$omp end parallel
128 
129  end subroutine hecmw_precond_diag_22_apply
130 
132  implicit none
133  if (associated(alu)) deallocate(alu)
134  nullify(alu)
135  initialized = .false.
136  end subroutine hecmw_precond_diag_22_clear
137 
138 end module hecmw_precond_diag_22
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
subroutine, public hecmw_precond_diag_22_clear()
subroutine, public hecmw_precond_diag_22_setup(hecMAT)
subroutine, public hecmw_precond_diag_22_apply(WW)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal