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