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 38: Line 38:
         implicit none
         implicit none
-
         interface BicgStab
+
         interface BICGStab
-
                 function BicgStab(A,b) result(x)
+
                 function BICGStab(A,b) result(x)
                         real    (kind=8), intent(in )  :: A(:,:)
                         real    (kind=8), intent(in )  :: A(:,:)
                         real    (kind=8), intent(in )  :: b( : )
                         real    (kind=8), intent(in )  :: b( : )
                         real    (kind=8)                :: x(1:size(b, dim=1))
                         real    (kind=8)                :: x(1:size(b, dim=1))
-
                 end function BicgStab
+
                 end function BICGStab
-
         end interface ! BicgStab
+
         end interface ! BICGStab
         integer (kind=4), parameter            :: m=4, n=4
         integer (kind=4), parameter            :: m=4, n=4
Line 51: Line 51:
!-------------------------------A,b DEFINATION----------------------------------!
!-------------------------------A,b DEFINATION----------------------------------!
-
         A(1,:) = (/1.0d0, 1.0d0, 1.0d0, 1.0d0/) ![x(1)] = b(1) |
+
         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(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(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) |
+
         A(4,:) = (/ 3.0d0, 4.0d0, 1.0d0, 2.0d0/) ![x(4)] = b(4) |
-
         b( : ) = (/5.0d0, 9.0d0,12.0d0,13.0d0/)
+
         b( : ) = (/ 4.0d0, 7.0d0, 9.0d0,10.0d0/)
!-----------------------------END A,b DEFINATION--------------------------------!
!-----------------------------END A,b DEFINATION--------------------------------!
-
         x = BicgStab(A,b)
+
         x = BICGStab(A,b)
-
         print*, "X is:",x
+
         print*, x
end program main
end program main
Line 68: Line 68:
-
function BicgStab(A,b) result(x)
+
function BICGStab(A,b) result(x)
         implicit none
         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-------------------------------!
!--------------------------PARAMETER AND VARIABLE-------------------------------!
Line 91: Line 77:
         real    (kind=8), dimension(1:size(b, dim=1))  :: r, rs, v, p, s, t
         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), parameter                    :: e = 1d-10
         real    (kind=8)                                :: rho      , rho_prev
         real    (kind=8)                                :: rho      , rho_prev
         real    (kind=8)                                :: alpha    , omega  , beta
         real    (kind=8)                                :: alpha    , omega  , beta
Line 97: Line 83:
         real    (kind=8)                                :: summesion, temp
         real    (kind=8)                                :: summesion, temp
-
        integer                                        :: m,n
 
-
        integer                                        :: i,j
 
         integer                                        :: it=0,err
         integer                                        :: it=0,err
!------------------------END PARAMETER AND VARIABLE-----------------------------!   
!------------------------END PARAMETER AND VARIABLE-----------------------------!   
-
         m = size(A, dim=1)
+
         if(size(A, dim=1) /= size(A, dim=2)) stop &
-
        n = size(A, dim=2)
+
        "Error: Improper dimension of matrix A in 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
 
!-------------------------------------------------------!
!-------------------------------------------------------!
-
 
+
        x  = 0.0d0                                      !-------> INITIAL GUESS
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         rho_prev = rho                                  !-------> STEP 5
+
         r  = b - matmul(A,x)                            !-------> LINE 1
-
         rho      = dot_prod(rs,r)                      !
+
         rs = r                                         !
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         beta    = (rho/rho_prev) * (alpha/omega)        !-------> STEP 6
+
         rho   = 1.0d0; alpha = 1.0d0; omega = 1.0d0    !-------> LINE 2
!-------------------------------------------------------!
!-------------------------------------------------------!
-
         p      = r + beta * (p - omega*v)              !-------> STEP 7
+
         v  = 0.0d0; p = 0.0d0                          !-------> LINE 3
-
!-------------------------------------------------------!
+
!                                                       !
-
        v      = mat_vec_mul(A,p)                      !-------> STEP 8
+
         norm_r = dsqrt(dot_product(r,r))               !
-
!-------------------------------------------------------!
+
         norm_b = dsqrt(dot_product(b,b))                !
-
         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))
+
         do while(norm_r .GT. e*norm_b)                                 !-------> START OF LOOP
-
        norm_b = dsqrt(dot_prod(b,b))  
+
-
         it = it + 1
+
                !-------------------------------------------------------!
 +
                        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 LOOP
         end do !END LOOP
-
 
+
         print*, "Iteration Required :", it
-
         print*, "Iter :",it
+
return
return
-
end function BicgStab
+
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
+
</pre>
</pre>

Revision as of 07:35, 21 June 2017

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


        !---------------------------------------------------------------------------------------!
        !       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( : ) = (/ 4.0d0, 7.0d0, 9.0d0,10.0d0/)
!-----------------------------END A,b DEFINATION--------------------------------!

        x = BICGStab(A,b)

        print*, x

end program main





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 LOOP

        print*, "Iteration Required :", it

return
end function BICGStab
Retrieved from "https://cfd-online.com/Wiki/Sample_code_for_BiCGSTAB_-_Fortran_90"
My wiki