From CFD-Wiki
!---------------------------------------------------------------------------------------!
! 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! !
! This code is tested with gFortran 5.01 !
!---------------------------------------------------------------------------------------!
module BiCGStab_mod
implicit none
contains
function BiCGStab(A,b) result(x)
implicit none
!--------------------------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-10
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 :: 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.0d0 !-------> INITIAL GUESS
!-------------------------------------------------------!
r = b - matmul(A,x) !-------> LINE 1
rs = r !
!-------------------------------------------------------!
rho = 1.0d0; alpha = 1.0d0; omega = 1.0d0 !-------> LINE 2
!-------------------------------------------------------!
v = 0.0d0; p = 0.0d0 !-------> LINE 3
! !
norm_r = dsqrt(dot_product(r,r)) !
norm_b = dsqrt(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 = dsqrt(dot_product(r,r))
norm_b = dsqrt(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=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( : ) = (/ 4.0d0, 7.0d0, 9.0d0,10.0d0/)
!---------------------------END A,b DEFINATION------------------------------!
!----CALL BiCGStab func-----!
! !
x = BiCGStab(A,b) !
! !
!---------------------------!
print*, x
end program main
