Sample code for BiCGSTAB - Fortran 90
From CFD-Wiki
(Difference between revisions)
(9 intermediate revisions not shown) | |||
Line 1: | Line 1: | ||
<pre> | <pre> | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
!---------------------------------------------------------------------------------------! | !---------------------------------------------------------------------------------------! | ||
- | ! | + | ! 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) | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | function | + | |
implicit none | implicit none | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
!--------------------------PARAMETER AND VARIABLE-------------------------------! | !--------------------------PARAMETER AND VARIABLE-------------------------------! | ||
- | real (kind= | + | real (kind=wp), intent(in ) :: A (:,:) |
- | real (kind= | + | real (kind=wp), intent(in ) :: b ( : ) |
- | real (kind= | + | real (kind=wp), dimension(1:size(b, dim=1)) :: x |
- | real (kind= | + | real (kind=wp), dimension(1:size(b, dim=1)) :: r, rs, v, p, s, t |
- | real (kind= | + | real (kind=wp), parameter :: e = 1d-33 |
- | real (kind= | + | real (kind=wp) :: rho , rho_prev |
- | real (kind= | + | real (kind=wp) :: alpha , omega , beta |
- | real (kind= | + | real (kind=wp) :: norm_r , norm_b |
- | real (kind= | + | real (kind=wp) :: summesion, temp |
- | integer | + | integer (kind=sp) :: it=0,err |
- | + | ||
- | + | ||
!------------------------END PARAMETER AND VARIABLE-----------------------------! | !------------------------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 | + | 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 | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | end | + | |
</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