From CFD-Wiki
!-----------------------------------MIT LICENCE-----------------------------------------!
! !
! Copyright (c) 2015, Biswajit Ghosh !
! http://www.biswa.me/cfd !
! 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 !
! THE SOFTWARE. !
!---------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------!
! NOTE : THIS CODE IS DESIGNED ONLY FOR LEARNING PURPOSE. !
! IMPLEMENTATION IN PRACTICAL PROBLEM MAY RESULT CONVERGENCE ISSUE !
!---------------------------------------------------------------------------------------!
module var
implicit none
integer, parameter :: n = 4 !-------> size of the matrices
double precision, parameter :: e = 1e-20 !-------> maximum permiable error
double precision, dimension(1:n,1:n) :: A
double precision, dimension(1:n ) :: b,x
double precision, dimension(1:n ) :: r, rs, v, p, s, t
double precision :: rho, rho_prev, alpha, omega, beta
double precision :: norm_r, norm_b
end module var
!---------------------------------------------------------------------------------------!
! this code is based on the algo given in wikipedia : !
! http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method !
!---------------------------------------------------------------------------------------!
!--------------------------------------main program block---------------------------------------!
program BiCGSTAB
use var
implicit none
integer :: i,j,it
double precision :: summesion, temp
call ab_def !---------> defination of matrics A,b
call init !---------> Setting initial values of variable
it = 0 !---------> with 'it', we count the number of iteration
!*******************************************************!
do while(norm_r .GT. e*norm_b) !-----> convergence criteria
!-------------------------------------------------------! step 01
rho_prev = rho !---------> rho_prev stores the
rho = 0.0d0 ! value of rho at n-1 th step
do i = 1,n
rho = rho + rs(i)*r(i)
end do
!------------------------------------------------------! step 02
beta = (rho/rho_prev) * (alpha/omega)
!------------------------------------------------------! step 03
p = r + beta * (p - omega*v)
!------------------------------------------------------! step 04
do i = 1,n
summesion = 0.0d0
do j = 1,n
summesion = summesion + A(i,j)*p(j)
end do
v(i) = summesion
end do
!------------------------------------------------------! step 05
summesion = 0.0d0
do i = 1,n
summesion = summesion + rs(i)*v(i)
end do
alpha = rho/summesion
!------------------------------------------------------! step 06
s = r - alpha*v
!------------------------------------------------------! step 07
do i = 1,n
summesion = 0.0d0
do j = 1,n
summesion = summesion + A(i,j)*s(j)
end do
t(i) = summesion
end do
!-----------------------------------------------------! step 08
summesion = 0.0d0
do i = 1,n
summesion = summesion + t(i)*s(i)
end do
temp = summesion
summesion = 0.0d0
do i = 1,n
summesion = summesion + t(i)*t(i)
end do
omega = temp/summesion
!-----------------------------------------------------! step 09
x = x + alpha*p + omega*s
!-----------------------------------------------------! step 11
r = s - omega*t
!-----------------------------------------------------! step 10
call norm
it = it + 1
end do !main loop
!*****************************************************!
write(*,*) 'the solution is',x
write(*,*) 'iteration required = ', it
end program BiCGSTAB
!-----------------------------------------------------------------------------------------!
!-----------------------------------subroutine init---------------------------------------!
subroutine init
use var
implicit none
integer :: i,j
double precision :: summesion
x = 0.0d0
!---------------------------------------------! step 01
do i = 1,n
summesion = 0.0d0
do j = 1,n
summesion = summesion + A(i,j)*x(j)
end do
r(i) = b(i) - summesion
end do
!---------------------------------------------! step 02
rs = r
!---------------------------------------------! step 03
rho = 1.0d0
alpha = 1.0d0
omega = 1.0d0
!---------------------------------------------! step 04
v = 0.0d0
p = 0.0d0
!---------------------------------------------!
call norm
end subroutine init
subroutine norm
use var
implicit none
integer :: i
double precision :: summesion
!-------------------------------------------
summesion = 0.0d0
do i = 1,n
summesion = summesion + r(i)*r(i)
end do
norm_r = dsqrt(summesion)
!-------------------------------------------
summesion = 0.0d0
do i = 1,n
summesion = summesion + b(i)*b(i)
end do
norm_b = dsqrt(summesion)
end subroutine norm
!---------------------------------subroutine ab_def---------------------------------------!
subroutine ab_def
use var
implicit none
!---------------------first row
A(1,1) = 1.0d0
A(1,2) = 1.0d0
A(1,3) = 1.0d0
A(1,4) = 1.0d0
!---------------------end first row
A(2,1) = 2.0d0
A(2,2) = 3.0d0
A(2,3) = 1.0d0
A(2,4) = 1.0d0
A(3,1) = 3.0d0
A(3,2) = 4.0d0
A(3,3) = 1.0d0
A(3,4) = 1.0d0
A(4,1) = 3.0d0
A(4,2) = 4.0d0
A(4,3) = 1.0d0
A(4,4) = 2.0d0
b(1) = 23.0d0
b(2) = 13.0d0
b(3) = 18.0d0
b(4) = 23.0d0
end subroutine ab_def
