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_44
9!C***
10!C
11module hecmw_precond_DIAG_44
12  use  hecmw_util
13
14  private
15
16  public:: hecmw_precond_DIAG_44_setup
17  public:: hecmw_precond_DIAG_44_apply
18  public:: hecmw_precond_DIAG_44_clear
19
20  integer(kind=kint) :: N
21  real(kind=kreal), pointer :: ALU(:) => null()
22
23  logical, save :: INITIALIZED = .false.
24
25contains
26
27  subroutine hecmw_precond_DIAG_44_setup(hecMAT)
28    use hecmw_matrix_misc
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(4,4), PW(4)
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
40      call hecmw_precond_DIAG_44_clear()
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(16*NP))
49    ALU = 0.d0
50
51    do ii= 1, N
52      ALU(16*ii-15) = D(16*ii-15)
53      ALU(16*ii-14) = D(16*ii-14)
54      ALU(16*ii-13) = D(16*ii-13)
55      ALU(16*ii-12) = D(16*ii-12)
56      ALU(16*ii-11) = D(16*ii-11)
57      ALU(16*ii-10) = D(16*ii-10)
58      ALU(16*ii-9 ) = D(16*ii-9 )
59      ALU(16*ii-8 ) = D(16*ii-8 )
60      ALU(16*ii-7 ) = D(16*ii-7 )
61      ALU(16*ii-6 ) = D(16*ii-6 )
62      ALU(16*ii-5 ) = D(16*ii-5 )
63      ALU(16*ii-4 ) = D(16*ii-4 )
64      ALU(16*ii-3 ) = D(16*ii-3 )
65      ALU(16*ii-2 ) = D(16*ii-2 )
66      ALU(16*ii-1 ) = D(16*ii-1 )
67      ALU(16*ii   ) = D(16*ii   )
68    enddo
69
70    do ii= 1, N
71      ALUtmp(1,1)= ALU(16*ii-15) * SIGMA_DIAG
72      ALUtmp(1,2)= ALU(16*ii-14)
73      ALUtmp(1,3)= ALU(16*ii-13)
74      ALUtmp(1,4)= ALU(16*ii-12)
75
76      ALUtmp(2,1)= ALU(16*ii-11)
77      ALUtmp(2,2)= ALU(16*ii-10) * SIGMA_DIAG
78      ALUtmp(2,3)= ALU(16*ii- 9)
79      ALUtmp(2,4)= ALU(16*ii- 8)
80
81      ALUtmp(3,1)= ALU(16*ii- 7)
82      ALUtmp(3,2)= ALU(16*ii- 6)
83      ALUtmp(3,3)= ALU(16*ii- 5) * SIGMA_DIAG
84      ALUtmp(3,4)= ALU(16*ii- 4)
85
86      ALUtmp(4,1)= ALU(16*ii- 3)
87      ALUtmp(4,2)= ALU(16*ii- 2)
88      ALUtmp(4,3)= ALU(16*ii- 1)
89      ALUtmp(4,4)= ALU(16*ii   ) * SIGMA_DIAG
90
91      do k= 1, 4
92        ALUtmp(k,k)= 1.d0/ALUtmp(k,k)
93        do i= k+1, 4
94          ALUtmp(i,k)= ALUtmp(i,k) * ALUtmp(k,k)
95          do j= k+1, 4
96            PW(j)= ALUtmp(i,j) - ALUtmp(i,k)*ALUtmp(k,j)
97          enddo
98          do j= k+1, 4
99            ALUtmp(i,j)= PW(j)
100          enddo
101        enddo
102      enddo
103
104      ALU(16*ii-15)= ALUtmp(1,1)
105      ALU(16*ii-14)= ALUtmp(1,2)
106      ALU(16*ii-13)= ALUtmp(1,3)
107      ALU(16*ii-12)= ALUtmp(1,4)
108      ALU(16*ii-11)= ALUtmp(2,1)
109      ALU(16*ii-10)= ALUtmp(2,2)
110      ALU(16*ii- 9)= ALUtmp(2,3)
111      ALU(16*ii- 8)= ALUtmp(2,4)
112      ALU(16*ii- 7)= ALUtmp(3,1)
113      ALU(16*ii- 6)= ALUtmp(3,2)
114      ALU(16*ii- 5)= ALUtmp(3,3)
115      ALU(16*ii- 4)= ALUtmp(3,4)
116      ALU(16*ii- 3)= ALUtmp(4,1)
117      ALU(16*ii- 2)= ALUtmp(4,2)
118      ALU(16*ii- 1)= ALUtmp(4,3)
119      ALU(16*ii   )= ALUtmp(4,4)
120    enddo
121
122    INITIALIZED = .true.
123    hecMAT%Iarray(98) = 0 ! symbolic setup done
124    hecMAT%Iarray(97) = 0 ! numerical setup done
125  end subroutine hecmw_precond_DIAG_44_setup
126
127  subroutine hecmw_precond_DIAG_44_apply(WW)
128    implicit none
129    real(kind=kreal), intent(inout) :: WW(:)
130    integer(kind=kint) :: i
131    real(kind=kreal) :: X1, X2, X3, X4
132
133    !C
134    !C== Block SCALING
135    do i= 1, N
136      X1= WW(4*i-3)
137      X2= WW(4*i-2)
138      X3= WW(4*i-1)
139      X4= WW(4*i  )
140      X2= X2 -ALU(16*i-11)*X1
141      X3= X3 -ALU(16*i- 7)*X1 -ALU(16*i- 6)*X2
142      X4= X4 -ALU(16*i- 3)*X1 -ALU(16*i- 2)*X2 -ALU(16*i- 1)*X3
143      X4= ALU(16*i   )* X4
144      X3= ALU(16*i- 5)*(X3 -ALU(16*i- 4)*X4)
145      X2= ALU(16*i-10)*(X2 -ALU(16*i- 8)*X4 -ALU(16*i- 9)*X3)
146      X1= ALU(16*i-15)*(X1 -ALU(16*i-12)*X4 -ALU(16*i-13)*X3 -ALU(16*i-14)*X2)
147      WW(4*i-3)= X1
148      WW(4*i-2)= X2
149      WW(4*i-1)= X3
150      WW(4*i  )= X4
151    enddo
152  end subroutine hecmw_precond_DIAG_44_apply
153
154  subroutine hecmw_precond_DIAG_44_clear()
155    implicit none
156    if (associated(ALU)) deallocate(ALU)
157    nullify(ALU)
158    INITIALIZED = .false.
159  end subroutine hecmw_precond_DIAG_44_clear
160
161end module     hecmw_precond_DIAG_44
162