Sample code for BiCGSTAB - Fortran 90
From CFD-Wiki
(Difference between revisions)
m |
|||
Line 38: | Line 38: | ||
implicit none | implicit none | ||
- | interface | + | interface BiCGStab |
- | function | + | 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 | + | end function BiCGStab |
- | end interface ! | + | end interface ! BiCGStab |
integer (kind=4), parameter :: m=4, n=4 | integer (kind=4), parameter :: m=4, n=4 | ||
Line 58: | Line 58: | ||
!-----------------------------END A,b DEFINATION--------------------------------! | !-----------------------------END A,b DEFINATION--------------------------------! | ||
- | + | !--------------CALL BiCGStab function--------------! | |
+ | ! ! | ||
+ | x = BiCGStab(A,b) ! | ||
+ | ! ! | ||
+ | !--------------------------------------------------! | ||
print*, x | print*, x | ||
Line 68: | Line 72: | ||
- | function | + | function BiCGStab(A,b) result(x) |
implicit none | implicit none | ||
Line 87: | Line 91: | ||
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 | + | "Error: Improper dimension of matrix A in BiCGStab." |
!-------------------------------------------------------! | !-------------------------------------------------------! | ||
Line 138: | Line 142: | ||
return | return | ||
- | end function | + | end function BiCGStab |
</pre> | </pre> |
Revision as of 15:14, 23 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--------------------------------! !--------------CALL BiCGStab function--------------! ! ! 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