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
11module hecmw_precond_DIAG_33
12  use  hecmw_util
13
14  private
15
16  public:: hecmw_precond_DIAG_33_setup
17  public:: hecmw_precond_DIAG_33_apply
18  public:: hecmw_precond_DIAG_33_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_33_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(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
40      call hecmw_precond_DIAG_33_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(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
125  subroutine hecmw_precond_DIAG_33_apply(WW)
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
154  subroutine hecmw_precond_DIAG_33_clear()
155    implicit none
156    if (associated(ALU)) deallocate(ALU)
157    nullify(ALU)
158    INITIALIZED = .false.
159  end subroutine hecmw_precond_DIAG_33_clear
160
161end module     hecmw_precond_DIAG_33
162