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
m
Line 1: Line 1:
<pre>
<pre>
-
!-----------------------------------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 !
 
-
!---------------------------------------------------------------------------------------!
 
-
 
         !---------------------------------------------------------------------------------------!
         !---------------------------------------------------------------------------------------!
Line 34: Line 11:
         !---------------------------------------------------------------------------------------!
         !---------------------------------------------------------------------------------------!
 +
module BiCGStab_mod
 +
    implicit none
 +
 
 +
contains
-
program main
+
    function BiCGStab(A,b) result(x)
-
        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--------------------------------!
+
-
 
+
-
!--------------CALL BiCGStab function--------------!
+
-
!                                                  !
+
-
                x = BiCGStab(A,b)                  !
+
-
!                                                  !
+
-
!--------------------------------------------------!
+
-
 
+
-
        print*, x
+
-
 
+
-
end program main
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
function BiCGStab(A,b) result(x)
+
         implicit none
         implicit none
Line 107: Line 51:
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         do while(norm_r .GT. e*norm_b)                                 !-------> START OF LOOP
+
         do while(norm_r .GT. e*norm_b)                             !-------> START OF LOOP
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        rho_prev = rho                                 !-------> LINE 5
+
                rho_prev = rho                                     !-------> LINE 5
-
                        rho      = dot_product(rs,r)                   !
+
                rho      = dot_product(rs,r)                       !
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        beta    = (rho/rho_prev) * (alpha/omega)       !-------> LINE 6
+
                beta    = (rho/rho_prev) * (alpha/omega)           !-------> LINE 6
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        p        = r + beta * (p - omega*v)             !-------> LINE 7
+
                p        = r + beta * (p - omega*v)                 !-------> LINE 7
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        v        = matmul(A,p)                         !-------> LINE 8
+
                v        = matmul(A,p)                             !-------> LINE 8
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        alpha    = rho/dot_product(rs,v)               !-------> LINE 9
+
                alpha    = rho/dot_product(rs,v)                   !-------> LINE 9
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        s        = r - alpha*v                         !-------> LINE 10
+
                s        = r - alpha*v                             !-------> LINE 10
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        t        = matmul(A,s)                         !-------> LINE 11
+
                t        = matmul(A,s)                             !-------> LINE 11
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        omega    = dot_product(t,s)/dot_product(t,t)   !-------> LINE 12
+
                omega    = dot_product(t,s)/dot_product(t,t)       !-------> LINE 12
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        x        = x + alpha*p + omega*s               !-------> LINE 13
+
                x        = x + alpha*p + omega*s                   !-------> LINE 13
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
-
                        r        = s - omega*t                         !-------> LINE 17
+
                r        = s - omega*t                             !-------> LINE 17
-
                !-------------------------------------------------------!
+
            !-------------------------------------------------------!
 +
                norm_r  = dsqrt(dot_product(r,r))
 +
                norm_b  = dsqrt(dot_product(b,b)) 
          
          
-
                        norm_r  = dsqrt(dot_product(r,r))
+
                it = it + 1
-
                        norm_b  = dsqrt(dot_product(b,b)) 
+
-
       
+
-
                        it = it + 1
+
-
         end do                                                         !-------> END OF LOOP
+
         end do                                                     !-------> END OF LOOP
         print*, "Iteration Required :", it
         print*, "Iteration Required :", it
return
return
-
end function BiCGStab
+
end function BiCGStab    
 +
 
 +
end module BiCGStab_mod
 +
 
 +
 
 +
program main
 +
    use BiCGStab_mod
 +
    implicit none
 +
 
 +
    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------------------------------!
 +
 
 +
!----CALL BiCGStab func-----!
 +
!                          !
 +
    x = BiCGStab(A,b)      !
 +
!                          !
 +
!---------------------------!
 +
 
 +
    print*, x
 +
 
 +
end program main
</pre>
</pre>

Revision as of 07:33, 25 August 2017


        !---------------------------------------------------------------------------------------!
        !       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!                   !
        !---------------------------------------------------------------------------------------!

module BiCGStab_mod
    implicit none
   
contains

    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 OF LOOP

        print*, "Iteration Required :", it

return
end function BiCGStab     

end module BiCGStab_mod


program main
    use BiCGStab_mod
    implicit none

    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------------------------------!

!----CALL BiCGStab func-----!
!                           !
    x = BiCGStab(A,b)       !
!                           !
!---------------------------!

    print*, x

end program main
Retrieved from "https://cfd-online.com/Wiki/Sample_code_for_BiCGSTAB_-_Fortran_90"
My wiki