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
 
(16 intermediate revisions not shown)
Line 1: Line 1:
<pre>
<pre>
-
!-----------------------------------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.                         !
+
        !                     This code is tested with gFortran 5.01                          !
-
! IMPLEMENTATION IN PRACTICAL PROBLEM MAY RESULT CONVERGENCE ISSUE !
+
        !                                   (C) Biswa                                          !
-
!---------------------------------------------------------------------------------------!
+
        !---------------------------------------------------------------------------------------!
-
+
-
module var
+
module BiCGStab_mod
-
implicit none
+
    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
+
    integer (kind=2 ), parameter    :: sp = kind(1.000)
-
double precision   :: rho, rho_prev, alpha, omega, beta
+
    integer (kind=2 ), parameter    :: dp = kind(1.0d0)
-
double precision :: norm_r, norm_b
+
    integer (kind=2 ), parameter    :: qp = kind(1.0q0)
-
end module var
+
    integer (kind=2 ), parameter    :: wp = dp        ! can be changed to "qp" to get quad precision.
-
!---------------------------------------------------------------------------------------!
+
contains
-
! this code is based on the algo given in wikipedia : !
+
-
! http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method !
+
-
!---------------------------------------------------------------------------------------!
+
 +
    function BiCGStab(A,b) result(x)
 +
        implicit none
-
!--------------------------------------main program block---------------------------------------!
+
!--------------------------PARAMETER AND VARIABLE-------------------------------!
-
program BiCGSTAB
+
        real    (kind=wp), intent(in )                  :: A (:,:)
-
use var
+
        real    (kind=wp), intent(in )                  :: b ( : )
-
implicit none
+
        real    (kind=wp), dimension(1:size(b, dim=1))  :: x
-
integer :: i,j,it
+
        real    (kind=wp), dimension(1:size(b, dim=1))  :: r, rs, v, p, s, t
-
double precision :: summesion, temp
+
        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
-
call ab_def !---------> defination of matrics A,b
+
        integer (kind=sp)                                :: it=0,err
-
+
!------------------------END PARAMETER AND VARIABLE-----------------------------! 
-
call init !---------> Setting initial values of variable
+
-
it = 0 !---------> with 'it', we count the number of iteration
+
        if(size(A, dim=1) /= size(A, dim=2)) stop &
 +
        "Error: Improper dimension of matrix A in BiCGStab."
-
!*******************************************************!
 
-
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)
+
!-------------------------------------------------------!
-
+
        x  = 0.0_wp                                    !-------> INITIAL GUESS
-
!------------------------------------------------------! step 03
+
!-------------------------------------------------------!
 +
        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))                !
 +
!-------------------------------------------------------!
-
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
+
        do while(norm_r .GT. e*norm_b)                         !-------> START OF LOOP
-
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
+
            rho_prev = rho                                      !-------> LINE 5
-
summesion = summesion + t(i)*s(i)
+
            rho      = dot_product(rs,r)                        !
-
end do
+
        !-------------------------------------------------------!
-
temp = summesion
+
            beta    = (rho/rho_prev) * (alpha/omega)           !-------> LINE 6
-
+
        !-------------------------------------------------------!
-
summesion  = 0.0d0
+
            p        = r + beta * (p - omega*v)                !-------> LINE 7
-
do i = 1,n
+
        !-------------------------------------------------------!
-
summesion = summesion + t(i)*t(i)
+
            v        = matmul(A,p)                              !-------> LINE 8
-
end do
+
        !-------------------------------------------------------!
-
+
            alpha    = rho/dot_product(rs,v)                   !-------> LINE 9
-
omega = temp/summesion
+
        !-------------------------------------------------------!
-
+
            s        = r - alpha*v                              !-------> LINE 10
-
!-----------------------------------------------------! step 09
+
        !-------------------------------------------------------!
-
+
            t       = matmul(A,s)                             !-------> LINE 11
-
x = x + alpha*p + omega*s
+
        !-------------------------------------------------------!
-
+
            omega   = dot_product(t,s)/dot_product(t,t)        !-------> LINE 12
-
!-----------------------------------------------------! step 11
+
        !-------------------------------------------------------!
 +
            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                                        !
 +
        !-------------------------------------------------------!
-
r = s - omega*t
+
        end do                                                      !-------> END OF LOOP
-
+
-
!-----------------------------------------------------! step 10
+
-
+
-
call norm
+
-
+
-
it = it + 1
+
-
+
-
end do !main loop
+
-
!*****************************************************!
+
-
+
-
write(*,*) 'the solution is',x
+
-
write(*,*) 'iteration required = ', it
+
-
+
-
end program BiCGSTAB
+
 +
        print*, "Iteration Required :", it
-
!-----------------------------------------------------------------------------------------!
+
return
 +
end function BiCGStab   
-
!-----------------------------------subroutine init---------------------------------------!
+
end module BiCGStab_mod
-
subroutine init
 
-
use var
 
-
implicit none
 
-
integer :: i,j
 
-
double precision :: summesion
 
-
 
-
x = 0.0d0
 
-
!---------------------------------------------! step 01
+
program main
-
do i = 1,n
+
    use BiCGStab_mod
-
summesion = 0.0d0
+
    implicit none
-
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
+
    integer (kind=sp), parameter            :: m=4, n=4
-
+
    real    (kind=wp), dimension(1:m,1:n)    :: A
-
end subroutine init
+
    real    (kind=wp), dimension(1:m    )    :: x_calculated,x_actual,b
-
subroutine norm
+
!-----------------------------A,b DEFINATION--------------------------------!
-
use var
+
    A      (1,:) = [ 1.0_wp, 1.0_wp, 1.0_wp, 1.0_wp] ![x(1)] = b(1) |
-
implicit none
+
    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) |
-
integer :: i
+
!  
-
double precision :: summesion
+
    x_actual( : ) = [ 5.0_wp, 8.0_wp,10.0_wp,12.0_wp]
-
!-------------------------------------------!
+
-
+
-
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
+
 +
!  settting b = Ax
 +
    b = matmul(A,x_actual)
 +
!---------------------------END A,b DEFINATION------------------------------!
-
!---------------------------------subroutine ab_def---------------------------------------!
+
!----------CALL BiCGStab func-----------!
 +
!                                      !
 +
    x_calculated = BiCGStab(A,b)        !
 +
!                                      !
 +
!---------------------------------------!
-
subroutine ab_def
+
    print*, "analytical solution =", x_actual
-
use var
+
    print*, "computed  solution =", x_calculated
-
implicit none
+
    print*, "Error in computed solution =", abs(x_actual-x_calculated)/x_actual
-
!---------------------first row
+
end program main
-
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
 
</pre>
</pre>

Latest revision as of 07:14, 9 May 2018


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