From CFD-Wiki
!-----------------------------------MIT LICENCE-----------------------------------------!
! !
! Copyright (c) 2015, Biswajit Ghosh !
! 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( : ) = (/ 4.0d0, 7.0d0, 9.0d0,10.0d0/)
!-----------------------------END A,b DEFINATION--------------------------------!
x = BICGStab(A,b)
print*, x
end program main
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 LOOP
print*, "Iteration Required :", it
return
end function BICGStab
