CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Sample code for BiCGSTAB - Fortran 90

Sample code for BiCGSTAB - Fortran 90

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
Line 24: Line 24:
!---------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------!
-
 
-
module var
+
        !---------------------------------------------------------------------------------------!
-
implicit none
+
        !      REFERENCE :                                                                     !
-
+
        !      The Improved BiCGStab Method for Large and Sparse Unsymmetric Linear            !
-
integer, parameter   :: n = 4 !-------> size of the matrices
+
        !      Systems on Parallel Distributed Memory Architectures; Yang et. al.              !
-
double precision, parameter :: e = 1e-20 !-------> maximum permiable error
+
        !                                                                                      !
-
double precision, dimension(1:n,1:n) :: A
+
        !      NOTE:                                                                           !
-
double precision, dimension(1:n    ) :: b,x
+
        !      This is the classic version of BICGStab. Stopping criteria of these code        !
 +
        !      needs additional treatment to avoid/detect failure condition!                  !
 +
        !---------------------------------------------------------------------------------------!
-
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
 
-
!---------------------------------------------------------------------------------------!
+
program main
-
! this code is based on the algo given in wikipedia : !
+
        implicit none
-
! http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method !
+
-
!---------------------------------------------------------------------------------------!
+
 +
        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
-
!--------------------------------------main program block---------------------------------------!
+
        integer (kind=4), parameter            :: m=4, n=4
-
program BiCGSTAB
+
        real    (kind=8), dimension(1:m,1:n)    :: A
-
use var
+
        real    (kind=8), dimension(1:m    )    :: x,b
-
implicit none
+
-
integer :: i,j,it
+
!-------------------------------A,b DEFINATION----------------------------------!
-
double precision :: summesion, temp
+
        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--------------------------------!
-
call ab_def !---------> defination of matrics A,b
+
        x = BicgStab(A,b)
-
+
-
call init !---------> Setting initial values of variable
+
-
it = 0 !---------> with 'it', we count the number of iteration
+
        print*, "X is:",x
 +
end program main
-
!-------------------Start of loop-----------------------!
 
-
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)
+
function BicgStab(A,b) result(x)
-
+
        implicit none
-
!------------------------------------------------------! step 04
+
-
do i = 1,n
+
        interface mat_vec_mul
-
summesion = 0.0d0
+
                function mat_vec_mul(a_mvm,b_mvm) result(c_mvm)
-
do j = 1,n
+
                      real    (kind=8), intent(in) :: a_mvm(:,:),b_mvm(:)
-
summesion = summesion + A(i,j)*p(j)
+
                      real    (kind=8)             :: c_mvm(1:size(b_mvm, dim=1))
-
end do
+
                end function mat_vec_mul
-
v(i) = summesion
+
        end interface ! mat_vec_mul
-
end do
+
-
+
-
!------------------------------------------------------! step 05
+
-
summesion = 0.0d0
+
        interface dot_prod
-
do i = 1,n
+
                function dot_prod(a_dp,b_dp) result(c_dp)
-
summesion = summesion + rs(i)*v(i)
+
                      real    (kind=8), intent(in) :: a_dp(:),b_dp(:)
-
end do
+
                      real    (kind=8)            :: c_dp
-
+
                end function dot_prod
-
alpha = rho/summesion
+
        end interface ! dot_prod
-
+
-
!------------------------------------------------------! step 06
+
-
s = r - alpha*v
+
!--------------------------PARAMETER AND VARIABLE-------------------------------!
-
+
        real    (kind=8), intent(in )                  :: A (:,:)
-
!------------------------------------------------------! step 07
+
        real    (kind=8), intent(in )                  :: b ( : )
 +
        real    (kind=8), dimension(1:size(b, dim=1))  :: x
-
do i = 1,n
+
        real    (kind=8), dimension(1:size(b, dim=1))  :: r, rs, v, p, s, t
-
summesion = 0.0d0
+
        real    (kind=8), parameter                    :: e = 1d-20
-
do j = 1,n
+
        real    (kind=8)                                :: rho      , rho_prev
-
summesion = summesion + A(i,j)*s(j)
+
        real    (kind=8)                                :: alpha    , omega  , beta
-
end do
+
        real    (kind=8)                               :: norm_r  , norm_b     
-
t(i) = summesion
+
        real    (kind=8)                               :: summesion, temp
-
end do
+
-
+
-
!-----------------------------------------------------! step 08
+
-
summesion  = 0.0d0
+
        integer                                        :: m,n
-
do i = 1,n
+
        integer                                        :: i,j
-
summesion = summesion + t(i)*s(i)
+
        integer                                        :: it=0,err
-
end do
+
!------------------------END PARAMETER AND VARIABLE-----------------------------!
-
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
+
        m = size(A, dim=1)
-
+
        n = size(A, dim=2)
-
!-----------------------------------------------------! step 10
+
-
+
-
call norm
+
-
+
-
it = it + 1
+
-
+
-
end do !main loop
+
-
!*****************************************************!
+
-
+
-
write(*,*) 'the solution is',x
+
-
write(*,*) 'iteration required = ', it
+
-
+
-
end program BiCGSTAB
+
 +
!---------------------------------------!
 +
        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
 +
!-------------------------------------------------------!
-
!-----------------------------------subroutine init---------------------------------------!
+
!-------------------------------------------------------!
 +
        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
 +
!-------------------------------------------------------!
-
subroutine init
+
        norm_r = dsqrt(dot_prod(r,r))
-
use var
+
        norm_b = dsqrt(dot_prod(b,b)) 
-
implicit none
+
-
integer :: i,j
+
-
double precision :: summesion
+
-
+
-
x = 0.0d0
+
-
!---------------------------------------------! step 01
+
        it = it + 1
-
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 do !END LOOP
-
+
-
end subroutine init
+
-
subroutine norm
 
-
use var
 
-
implicit none
 
-
integer :: i
+
        print*, "Iter :",it
-
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
+
 +
return
 +
end function BicgStab
-
!---------------------------------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
+
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))
-
A(2,1) = 2.0d0
+
        if (size(a_mvm, dim=2) /= size(b_mvm, dim=1)) stop &
-
A(2,2) = 3.0d0
+
        "Error: Improper dimension of matrix/vector in mat_vec_mul."
-
A(2,3) = 1.0d0
+
-
A(2,4) = 1.0d0
+
-
A(3,1) = 3.0d0
+
        c_mvm = matmul(a_mvm,b_mvm)
-
A(3,2) = 4.0d0
+
-
A(3,3) = 1.0d0
+
-
A(3,4) = 1.0d0
+
-
A(4,1) = 3.0d0
+
        return
-
A(4,2) = 4.0d0
+
end function mat_vec_mul
-
A(4,3) = 1.0d0
+
 
-
A(4,4) = 2.0d0
+
 
-
+
 
-
b(1) = 23.0d0
+
 
-
b(2) = 13.0d0
+
 
-
b(3) = 18.0d0
+
function dot_prod(a_dp,b_dp) result(c_dp)
-
b(4) = 23.0d0
+
        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
-
 
-
end subroutine ab_def
 
</pre>
</pre>

Revision as of 22:09, 20 June 2017

	!-----------------------------------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
Retrieved from "https://cfd-online.com/Wiki/Sample_code_for_BiCGSTAB_-_Fortran_90"
My wiki