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
 
(4 intermediate revisions not shown)
Line 2: Line 2:
         !---------------------------------------------------------------------------------------!
         !---------------------------------------------------------------------------------------!
-
         !       REFERENCE :                                                                    !
+
         !                     This code is tested with gFortran 5.01                          !
-
        !      The Improved BiCGStab Method for Large and Sparse Unsymmetric Linear            !
+
         !                                   (C) Biswa                                          !
-
        !      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!                  !
+
         !---------------------------------------------------------------------------------------!
         !---------------------------------------------------------------------------------------!
module BiCGStab_mod
module BiCGStab_mod
     implicit none
     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
contains
Line 20: Line 20:
!--------------------------PARAMETER AND VARIABLE-------------------------------!
!--------------------------PARAMETER AND VARIABLE-------------------------------!
-
         real    (kind=8), intent(in )                  :: A (:,:)
+
         real    (kind=wp), intent(in )                  :: A (:,:)
-
         real    (kind=8), intent(in )                  :: b ( : )
+
         real    (kind=wp), intent(in )                  :: b ( : )
-
         real    (kind=8), dimension(1:size(b, dim=1))  :: x
+
         real    (kind=wp), dimension(1:size(b, dim=1))  :: x
-
         real    (kind=8), dimension(1:size(b, dim=1))  :: r, rs, v, p, s, t
+
         real    (kind=wp), dimension(1:size(b, dim=1))  :: r, rs, v, p, s, t
-
         real    (kind=8), parameter                    :: e = 1d-10
+
         real    (kind=wp), parameter                    :: e = 1d-33
-
         real    (kind=8)                                :: rho      , rho_prev
+
         real    (kind=wp)                                :: rho      , rho_prev
-
         real    (kind=8)                                :: alpha    , omega  , beta
+
         real    (kind=wp)                                :: alpha    , omega  , beta
-
         real    (kind=8)                                :: norm_r  , norm_b       
+
         real    (kind=wp)                                :: norm_r  , norm_b       
-
         real    (kind=8)                                :: summesion, temp
+
         real    (kind=wp)                                :: summesion, temp
-
         integer                                         :: it=0,err
+
         integer (kind=sp)                                :: it=0,err
!------------------------END PARAMETER AND VARIABLE-----------------------------!   
!------------------------END PARAMETER AND VARIABLE-----------------------------!   
         if(size(A, dim=1) /= size(A, dim=2)) stop &
         if(size(A, dim=1) /= size(A, dim=2)) stop &
         "Error: Improper dimension of matrix A in BiCGStab."
         "Error: Improper dimension of matrix A in BiCGStab."
 +
 +
 +
 +
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         x  = 0.0d0                                      !-------> INITIAL GUESS
+
         x  = 0.0_wp                                    !-------> INITIAL GUESS
!-------------------------------------------------------!
!-------------------------------------------------------!
         r  = b - matmul(A,x)                            !-------> LINE 1
         r  = b - matmul(A,x)                            !-------> LINE 1
         rs = r                                          !
         rs = r                                          !
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         rho  = 1.0d0; alpha = 1.0d0; omega = 1.0d0    !-------> LINE 2
+
         rho  = 1.0_wp; alpha = 1.0_wp; omega = 1.0_wp  !-------> LINE 2
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         v  = 0.0d0; p  = 0.0d0                          !-------> LINE 3
+
         v  = 0.0_wp; p  = 0.0_wp                        !-------> LINE 3
!                                                      !
!                                                      !
-
         norm_r = dsqrt(dot_product(r,r))               !
+
         norm_r = sqrt(dot_product(r,r))                 !
-
         norm_b = dsqrt(dot_product(b,b))               !
+
         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)                        !
+
 
-
            !-------------------------------------------------------!
+
        do while(norm_r .GT. e*norm_b)                          !-------> START OF LOOP
-
                beta    = (rho/rho_prev) * (alpha/omega)          !-------> LINE 6
+
 
-
            !-------------------------------------------------------!
+
        !-------------------------------------------------------!
-
                p        = r + beta * (p - omega*v)                !-------> LINE 7
+
            rho_prev = rho                                      !-------> LINE 5
-
            !-------------------------------------------------------!
+
            rho      = dot_product(rs,r)                        !
-
                v        = matmul(A,p)                              !-------> LINE 8
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            beta    = (rho/rho_prev) * (alpha/omega)          !-------> LINE 6
-
                alpha    = rho/dot_product(rs,v)                    !-------> LINE 9
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            p        = r + beta * (p - omega*v)                !-------> LINE 7
-
                s        = r - alpha*v                              !-------> LINE 10
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            v        = matmul(A,p)                              !-------> LINE 8
-
                t        = matmul(A,s)                              !-------> LINE 11
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            alpha    = rho/dot_product(rs,v)                    !-------> LINE 9
-
                omega    = dot_product(t,s)/dot_product(t,t)        !-------> LINE 12
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            s        = r - alpha*v                              !-------> LINE 10
-
                x        = x + alpha*p + omega*s                    !-------> LINE 13
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            t        = matmul(A,s)                              !-------> LINE 11
-
                r        = s - omega*t                              !-------> LINE 17
+
        !-------------------------------------------------------!
-
            !-------------------------------------------------------!
+
            omega    = dot_product(t,s)/dot_product(t,t)        !-------> LINE 12
-
                norm_r  = dsqrt(dot_product(r,r))
+
        !-------------------------------------------------------!
-
                norm_b  = dsqrt(dot_product(b,b))  
+
            x        = x + alpha*p + omega*s                    !-------> LINE 13
-
          
+
        !-------------------------------------------------------!
-
                it = it + 1
+
            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
         end do                                                      !-------> END OF LOOP
Line 94: Line 103:
     implicit none
     implicit none
-
     integer (kind=4), parameter            :: m=4, n=4
+
     integer (kind=sp), parameter            :: m=4, n=4
-
     real    (kind=8), dimension(1:m,1:n)    :: A
+
     real    (kind=wp), dimension(1:m,1:n)    :: A
-
     real    (kind=8), dimension(1:m    )    :: x,b
+
     real    (kind=wp), dimension(1:m    )    :: x_calculated,x_actual,b
!-----------------------------A,b DEFINATION--------------------------------!
!-----------------------------A,b DEFINATION--------------------------------!
-
     A(1,:) = (/ 1.0d0, 1.0d0, 1.0d0, 1.0d0/) ![x(1)] = b(1) |
+
     A       (1,:) = [ 1.0_wp, 1.0_wp, 1.0_wp, 1.0_wp] ![x(1)] = b(1) |
-
     A(2,:) = (/ 2.0d0, 3.0d0, 1.0d0, 1.0d0/) ![x(2)] = b(2) |---> [A]{x}={b}
+
     A       (2,:) = [ 2.0_wp, 3.0_wp, 1.0_wp, 1.0_wp] ![x(2)] = b(2) |---> [A]{x}={b}
-
     A(3,:) = (/ 3.0d0, 4.0d0, 1.0d0, 1.0d0/) ![x(3)] = b(3) |
+
     A       (3,:) = [ 3.0_wp, 4.0_wp, 1.0_wp, 1.0_wp] ![x(3)] = b(3) |
-
     A(4,:) = (/ 3.0d0, 4.0d0, 1.0d0, 2.0d0/) ![x(4)] = b(4) |
+
     A       (4,:) = [ 3.0_wp, 4.0_wp, 1.0_wp, 2.0_wp] ![x(4)] = b(4) |
-
     b( : ) = (/ 4.0d0, 7.0d0, 9.0d0,10.0d0/)
+
 
 +
! 
 +
     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------------------------------!
!---------------------------END A,b DEFINATION------------------------------!
-
!----CALL BiCGStab func-----!
+
!----------CALL BiCGStab func-----------!
-
!                           !
+
!                                       !
-
     x = BiCGStab(A,b)       !
+
     x_calculated = BiCGStab(A,b)       !
-
!                           !
+
!                                       !
-
!---------------------------!
+
!---------------------------------------!
-
     print*, x
+
     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
end program main
</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