From CFD-Wiki
!-----------------------------------MIT LICENCE-----------------------------------------!
! !
! Copyright (c) 2015, Biswajit Ghosh !
! mebiswajithit@gmail.com !
! All rights reserved. !
! !
! Permission is hereby granted, free of charge, to any person obtaining a copy !
! of this software and associated documentation files (the "Software"), to deal !
! in the Software without restriction, including without limitation the rights !
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell !
! copies of the Software, and to permit persons to whom the Software is !
! furnished to do so, subject to the following conditions: !
! !
! The above copyright notice and this permission notice shall be included in !
! all copies or substantial portions of the Software. !
! !
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR !
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, !
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE !
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER !
! LIABILITY, WHE THER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, !
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN !
!---------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------!
! REFERENCE : !
! The Improved BiCGStab Method for Large and Sparse Unsymmetric Linear !
! Systems on Parallel Distributed Memory Architectures; Yang et. al. !
! !
! NOTE: !
! This is the classic version of BICGStab. Stopping criteria of these code !
! needs additional treatment to avoid/detect failure condition! !
!---------------------------------------------------------------------------------------!
program main
implicit none
interface BicgStab
function BicgStab(A,b) result(x)
real (kind=8), intent(in ) :: A(:,:)
real (kind=8), intent(in ) :: b( : )
real (kind=8) :: x(1:size(b, dim=1))
end function BicgStab
end interface ! BicgStab
integer (kind=4), parameter :: m=4, n=4
real (kind=8), dimension(1:m,1:n) :: A
real (kind=8), dimension(1:m ) :: x,b
!-------------------------------A,b DEFINATION----------------------------------!
A(1,:) = (/1.0d0, 1.0d0, 1.0d0, 1.0d0/) ![x(1)] = b(1) |
A(2,:) = (/2.0d0, 3.0d0, 1.0d0, 1.0d0/) ![x(2)] = b(2) |---> [A]{x}={b}
A(3,:) = (/3.0d0, 4.0d0, 1.0d0, 1.0d0/) ![x(3)] = b(3) |
A(4,:) = (/3.0d0, 4.0d0, 1.0d0, 2.0d0/) ![x(4)] = b(4) |
b( : ) = (/5.0d0, 9.0d0,12.0d0,13.0d0/)
!-----------------------------END A,b DEFINATION--------------------------------!
x = BicgStab(A,b)
print*, "X is:",x
end program main
function BicgStab(A,b) result(x)
implicit none
interface mat_vec_mul
function mat_vec_mul(a_mvm,b_mvm) result(c_mvm)
real (kind=8), intent(in) :: a_mvm(:,:),b_mvm(:)
real (kind=8) :: c_mvm(1:size(b_mvm, dim=1))
end function mat_vec_mul
end interface ! mat_vec_mul
interface dot_prod
function dot_prod(a_dp,b_dp) result(c_dp)
real (kind=8), intent(in) :: a_dp(:),b_dp(:)
real (kind=8) :: c_dp
end function dot_prod
end interface ! dot_prod
!--------------------------PARAMETER AND VARIABLE-------------------------------!
real (kind=8), intent(in ) :: A (:,:)
real (kind=8), intent(in ) :: b ( : )
real (kind=8), dimension(1:size(b, dim=1)) :: x
real (kind=8), dimension(1:size(b, dim=1)) :: r, rs, v, p, s, t
real (kind=8), parameter :: e = 1d-20
real (kind=8) :: rho , rho_prev
real (kind=8) :: alpha , omega , beta
real (kind=8) :: norm_r , norm_b
real (kind=8) :: summesion, temp
integer :: m,n
integer :: i,j
integer :: it=0,err
!------------------------END PARAMETER AND VARIABLE-----------------------------!
m = size(A, dim=1)
n = size(A, dim=2)
!---------------------------------------!
x = 0.0d0 !-------> GUESS X
!---------------------------------------!
r = b - mat_vec_mul(A,x) !-------> STEP 1
rs = r !
!---------------------------------------!
rho = 1.0d0 !
alpha = 1.0d0 !-------> STEP 2
omega = 1.0d0 !
!---------------------------------------!
v = 0.0d0 !-------> STEP 3
p = 0.0d0 !
! !
norm_r = dsqrt(dot_prod(r,r)) !
norm_b = dsqrt(dot_prod(b,b)) !
!---------------------------------------!
!-------------------START OF LOOP-----------------------!
do while(norm_r .GT. e*norm_b) !-------> STEP 4
!-------------------------------------------------------!
!-------------------------------------------------------!
rho_prev = rho !-------> STEP 5
rho = dot_prod(rs,r) !
!-------------------------------------------------------!
beta = (rho/rho_prev) * (alpha/omega) !-------> STEP 6
!-------------------------------------------------------!
p = r + beta * (p - omega*v) !-------> STEP 7
!-------------------------------------------------------!
v = mat_vec_mul(A,p) !-------> STEP 8
!-------------------------------------------------------!
alpha = rho/dot_prod(rs,v) !-------> STEP 9
!-------------------------------------------------------!
s = r - alpha*v !-------> STEP 10
!-------------------------------------------------------!
t = mat_vec_mul(A,s) !-------> STEP 11
!-------------------------------------------------------!
omega = dot_prod(t,s)/dot_prod(t,t) !-------> STEP 12
!-------------------------------------------------------!
x = x + alpha*p + omega*s !-------> STEP 13
!-------------------------------------------------------!
r = s - omega*t !-------> STEP 17
!-------------------------------------------------------!
norm_r = dsqrt(dot_prod(r,r))
norm_b = dsqrt(dot_prod(b,b))
it = it + 1
end do !END LOOP
print*, "Iter :",it
return
end function BicgStab
function mat_vec_mul(a_mvm,b_mvm) result(c_mvm)
implicit none
real (kind=8), intent(in) :: a_mvm(:,:),b_mvm(:)
real (kind=8) :: c_mvm(1:size(b_mvm, dim=1))
if (size(a_mvm, dim=2) /= size(b_mvm, dim=1)) stop &
"Error: Improper dimension of matrix/vector in mat_vec_mul."
c_mvm = matmul(a_mvm,b_mvm)
return
end function mat_vec_mul
function dot_prod(a_dp,b_dp) result(c_dp)
implicit none
real (kind=8), intent(in) :: a_dp(:),b_dp(:)
real (kind=8) :: c_dp
if (size(a_dp, dim=1) /= size(b_dp, dim=1)) stop &
"Error: Improper dimension of matrices in dot_prod."
c_dp = dot_product(a_dp,b_dp)
return
end function dot_prod
