Sample code for BiCGSTAB - Fortran 90
From CFD-Wiki
(Difference between revisions)
Line 24: | Line 24: | ||
!---------------------------------------------------------------------------------------! | !---------------------------------------------------------------------------------------! | ||
- | |||
- | + | !---------------------------------------------------------------------------------------! | |
- | + | ! 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( : ) = (/5.0d0, 9.0d0,12.0d0,13.0d0/) | ||
+ | !-----------------------------END A,b DEFINATION--------------------------------! | ||
- | + | x = BicgStab(A,b) | |
- | + | ||
- | + | ||
- | + | print*, "X is:",x | |
+ | end program main | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | + | function BicgStab(A,b) result(x) | |
- | + | 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-------------------------------! | |
- | + | 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-20 | |
- | + | 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 :: m,n | |
- | + | integer :: i,j | |
- | + | integer :: it=0,err | |
- | + | !------------------------END PARAMETER AND VARIABLE-----------------------------! | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | !-----------------------------------------------------! | + | |
- | + | m = size(A, dim=1) | |
- | + | n = size(A, dim=2) | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
+ | !---------------------------------------! | ||
+ | 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 | ||
+ | !-------------------------------------------------------! | ||
- | !----------------------------------- | + | !-------------------------------------------------------! |
+ | rho_prev = rho !-------> STEP 5 | ||
+ | rho = dot_prod(rs,r) ! | ||
+ | !-------------------------------------------------------! | ||
+ | beta = (rho/rho_prev) * (alpha/omega) !-------> STEP 6 | ||
+ | !-------------------------------------------------------! | ||
+ | p = r + beta * (p - omega*v) !-------> STEP 7 | ||
+ | !-------------------------------------------------------! | ||
+ | v = mat_vec_mul(A,p) !-------> STEP 8 | ||
+ | !-------------------------------------------------------! | ||
+ | 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)) | |
- | + | norm_b = dsqrt(dot_prod(b,b)) | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | it = it + 1 | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | end do !END LOOP | |
- | + | ||
- | end | + | |
- | |||
- | |||
- | |||
- | + | print*, "Iter :",it | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
+ | return | ||
+ | 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 22:09, 20 June 2017
!-----------------------------------MIT LICENCE-----------------------------------------! ! ! ! Copyright (c) 2015, Biswajit Ghosh ! ! mebiswajithit@gmail.com ! ! 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( : ) = (/5.0d0, 9.0d0,12.0d0,13.0d0/) !-----------------------------END A,b DEFINATION--------------------------------! x = BicgStab(A,b) print*, "X is:",x end program main function BicgStab(A,b) result(x) 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-------------------------------! 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-20 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 :: m,n integer :: i,j integer :: it=0,err !------------------------END PARAMETER AND VARIABLE-----------------------------! m = size(A, dim=1) n = size(A, dim=2) !---------------------------------------! 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 !-------------------------------------------------------! !-------------------------------------------------------! rho_prev = rho !-------> STEP 5 rho = dot_prod(rs,r) ! !-------------------------------------------------------! beta = (rho/rho_prev) * (alpha/omega) !-------> STEP 6 !-------------------------------------------------------! p = r + beta * (p - omega*v) !-------> STEP 7 !-------------------------------------------------------! v = mat_vec_mul(A,p) !-------> STEP 8 !-------------------------------------------------------! 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)) norm_b = dsqrt(dot_prod(b,b)) it = it + 1 end do !END LOOP print*, "Iter :",it return 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