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