From CFD-Wiki
!---------------------------------------------------------------------------------------!
! This code is tested with gFortran 5.01 !
! (C) Biswa !
!---------------------------------------------------------------------------------------!
module BiCGStab_mod
implicit none
integer (kind=2 ), parameter :: sp = kind(1.000)
integer (kind=2 ), parameter :: dp = kind(1.0d0)
integer (kind=2 ), parameter :: qp = kind(1.0q0)
integer (kind=2 ), parameter :: wp = dp ! can be changed to "qp" to get quad precision.
contains
function BiCGStab(A,b) result(x)
implicit none
!--------------------------PARAMETER AND VARIABLE-------------------------------!
real (kind=wp), intent(in ) :: A (:,:)
real (kind=wp), intent(in ) :: b ( : )
real (kind=wp), dimension(1:size(b, dim=1)) :: x
real (kind=wp), dimension(1:size(b, dim=1)) :: r, rs, v, p, s, t
real (kind=wp), parameter :: e = 1d-33
real (kind=wp) :: rho , rho_prev
real (kind=wp) :: alpha , omega , beta
real (kind=wp) :: norm_r , norm_b
real (kind=wp) :: summesion, temp
integer (kind=sp) :: it=0,err
!------------------------END PARAMETER AND VARIABLE-----------------------------!
if(size(A, dim=1) /= size(A, dim=2)) stop &
"Error: Improper dimension of matrix A in BiCGStab."
!-------------------------------------------------------!
x = 0.0_wp !-------> INITIAL GUESS
!-------------------------------------------------------!
r = b - matmul(A,x) !-------> LINE 1
rs = r !
!-------------------------------------------------------!
rho = 1.0_wp; alpha = 1.0_wp; omega = 1.0_wp !-------> LINE 2
!-------------------------------------------------------!
v = 0.0_wp; p = 0.0_wp !-------> LINE 3
! !
norm_r = sqrt(dot_product(r,r)) !
norm_b = sqrt(dot_product(b,b)) !
!-------------------------------------------------------!
do while(norm_r .GT. e*norm_b) !-------> START OF LOOP
!-------------------------------------------------------!
rho_prev = rho !-------> LINE 5
rho = dot_product(rs,r) !
!-------------------------------------------------------!
beta = (rho/rho_prev) * (alpha/omega) !-------> LINE 6
!-------------------------------------------------------!
p = r + beta * (p - omega*v) !-------> LINE 7
!-------------------------------------------------------!
v = matmul(A,p) !-------> LINE 8
!-------------------------------------------------------!
alpha = rho/dot_product(rs,v) !-------> LINE 9
!-------------------------------------------------------!
s = r - alpha*v !-------> LINE 10
!-------------------------------------------------------!
t = matmul(A,s) !-------> LINE 11
!-------------------------------------------------------!
omega = dot_product(t,s)/dot_product(t,t) !-------> LINE 12
!-------------------------------------------------------!
x = x + alpha*p + omega*s !-------> LINE 13
!-------------------------------------------------------!
r = s - omega*t !-------> LINE 17
!-------------------------------------------------------!
norm_r = sqrt(dot_product(r,r)) !
norm_b = sqrt(dot_product(b,b)) !
!-------------------------------------------------------!
it = it + 1 !
!-------------------------------------------------------!
end do !-------> END OF LOOP
print*, "Iteration Required :", it
return
end function BiCGStab
end module BiCGStab_mod
program main
use BiCGStab_mod
implicit none
integer (kind=sp), parameter :: m=4, n=4
real (kind=wp), dimension(1:m,1:n) :: A
real (kind=wp), dimension(1:m ) :: x_calculated,x_actual,b
!-----------------------------A,b DEFINATION--------------------------------!
A (1,:) = [ 1.0_wp, 1.0_wp, 1.0_wp, 1.0_wp] ![x(1)] = b(1) |
A (2,:) = [ 2.0_wp, 3.0_wp, 1.0_wp, 1.0_wp] ![x(2)] = b(2) |---> [A]{x}={b}
A (3,:) = [ 3.0_wp, 4.0_wp, 1.0_wp, 1.0_wp] ![x(3)] = b(3) |
A (4,:) = [ 3.0_wp, 4.0_wp, 1.0_wp, 2.0_wp] ![x(4)] = b(4) |
!
x_actual( : ) = [ 5.0_wp, 8.0_wp,10.0_wp,12.0_wp]
! settting b = Ax
b = matmul(A,x_actual)
!---------------------------END A,b DEFINATION------------------------------!
!----------CALL BiCGStab func-----------!
! !
x_calculated = BiCGStab(A,b) !
! !
!---------------------------------------!
print*, "analytical solution =", x_actual
print*, "computed solution =", x_calculated
print*, "Error in computed solution =", abs(x_actual-x_calculated)/x_actual
end program main
