From 25dbf151e1e27216ca2f6f94edf3603e58451710 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 29 Sep 2020 12:04:59 -0300 Subject: [PATCH] ALGEBRA: Rewrite slightly the tridiagonal solver. Signed-off-by: Grzegorz Kowal --- sources/algebra.F90 | 66 +++++++++++--------------------------- sources/interpolations.F90 | 45 ++++++++++++-------------- 2 files changed, 40 insertions(+), 71 deletions(-) diff --git a/sources/algebra.F90 b/sources/algebra.F90 index 517f8e4..7251e91 100644 --- a/sources/algebra.F90 +++ b/sources/algebra.F90 @@ -1093,83 +1093,55 @@ module algebra ! subroutine TRIDIAG: ! ------------------ ! -! Subroutine solves the system of tridiagonal linear equations. +! Subroutine solves the system of tridiagonal linear equations using +! Thomas algorithm. ! ! Arguments: ! ! n - the size of the matrix; -! l(:) - the vector of the lower off-diagonal values; -! d(:) - the vector of the diagonal values; -! u(:) - the vector of the upper off-diagonal values; -! r(:) - the vector of the right side values; +! a(:) - the vector of the lower off-diagonal values; +! b(:) - the vector of the diagonal values; +! c(:) - the vector of the upper off-diagonal values; +! d(:) - the vector of the right side values; ! x(:) - the solution vector; ! ! References: ! -! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., -! "Numerical Recipes in Fortran", -! Cambridge University Press, Cambridge, 1992 +! [1] https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm ! !=============================================================================== ! - subroutine tridiag(n, l, d, u, r, x, iret) + subroutine tridiag(n, a, b, c, d, x) -! import external procedures -! - use iso_fortran_env, only : error_unit - -! local variables are not implicit by default -! implicit none -! input/output arguments -! integer , intent(in) :: n - real(kind=8), dimension(n), intent(in) :: l, d, u, r + real(kind=8), dimension(n), intent(in) :: a, b, c, d real(kind=8), dimension(n), intent(out) :: x - integer , intent(out) :: iret -! local variables -! - integer :: i, j + integer :: i real(kind=8) :: t - real(kind=8), dimension(n) :: g - -! local parameters -! - character(len=*), parameter :: loc = 'ALGEBRA::tridiag()' + real(kind=8), dimension(n) :: cp, dp ! !------------------------------------------------------------------------------- ! ! decomposition and forward substitution ! - t = d(1) - x(1) = r(1) / t + cp(1) = c(1) / b(1) + dp(1) = d(1) / b(1) do i = 2, n - j = i - 1 - g(i) = u(j) / t - t = d(i) - l(i) * g(i) - if (abs(t) > 0.0d+00) then - x(i) = (r(i) - l(i) * x(j)) / t - else - write(error_unit,"('[',a,']: ',a)") trim(loc) & - , "Could not find solution!" - iret = 1 - return - end if + t = b(i) - a(i) * cp(i-1) + cp(i) = c(i) / t + dp(i) = (d(i) - a(i) * dp(i-1)) / t end do -! backsubstitution +! backward substitution ! + x(n) = dp(n) do i = n - 1, 1, -1 - j = i + 1 - x(i) = x(i) - g(j) * x(j) + x(i) = dp(i) - cp(i) * x(i+1) end do -! set return value to success -! - iret = 0 - !------------------------------------------------------------------------------- ! end subroutine tridiag diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 338e2a7..cf8b6e6 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -2631,7 +2631,6 @@ module interpolations ! local variables ! integer :: n, i, im1, ip1, im2, ip2 - integer :: iret real(kind=8) :: bl, bc, br, tt real(kind=8) :: wl, wc, wr, ww real(kind=8) :: ql, qc, qr @@ -2948,7 +2947,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) ! substitute the left-side values ! @@ -2956,7 +2955,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret) + call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n)) ! substitute the right-side values ! @@ -3027,7 +3026,6 @@ module interpolations ! local variables ! integer :: n, i, im1, ip1, im2, ip2 - integer :: iret real(kind=8) :: bl, bc, br, tt real(kind=8) :: wl, wc, wr, ww real(kind=8) :: ql, qc, qr @@ -3347,7 +3345,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) ! substitute the left-side values ! @@ -3355,7 +3353,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret) + call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n)) ! substitute the right-side values ! @@ -3426,7 +3424,6 @@ module interpolations ! local variables ! integer :: n, i, im1, ip1, im2, ip2 - integer :: iret real(kind=8) :: bl, bc, br real(kind=8) :: wl, wc, wr, ww real(kind=8) :: df, lq, l3, zt @@ -3766,7 +3763,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) ! substitute the left-side values ! @@ -3774,7 +3771,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret) + call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n)) ! substitute the right-side values ! @@ -4649,7 +4646,7 @@ module interpolations ! local variables ! - integer :: n, i, iret + integer :: n, i ! local arrays for derivatives ! @@ -4718,7 +4715,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -4755,7 +4752,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -4829,7 +4826,7 @@ module interpolations ! local variables ! - integer :: n, i, iret + integer :: n, i ! local arrays for derivatives ! @@ -4899,7 +4896,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -4936,7 +4933,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -5003,7 +5000,7 @@ module interpolations ! local variables ! - integer :: n, i, iret + integer :: n, i ! local arrays for derivatives ! @@ -5078,7 +5075,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -5117,7 +5114,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -5185,7 +5182,7 @@ module interpolations ! local variables ! - integer :: n, i, iret + integer :: n, i ! local arrays for derivatives ! @@ -5261,7 +5258,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -5300,7 +5297,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -5353,7 +5350,7 @@ module interpolations real(kind=8), dimension(:), intent(in) :: fc real(kind=8), dimension(:), intent(out) :: fl, fr - integer :: n, i, iret + integer :: n, i real(kind=8), dimension(size(fc)) :: fi real(kind=8), dimension(size(fc)) :: a, b, c @@ -5417,7 +5414,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiter ! @@ -5454,7 +5451,7 @@ module interpolations ! solve the tridiagonal system of equations ! - call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret) + call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) ! apply the monotonicity preserving limiter !