!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2024 Grzegorz Kowal !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation, either version 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: ALGEBRA !! !! This module provides all sort of algebra subroutines. !! !!****************************************************************************** ! module algebra ! include external procedures ! use iso_fortran_env, only : real32, real64, real128 ! module variables are not implicit by default ! implicit none ! maximum real kind ! integer, parameter :: max_prec = max(real32, real64, real128) ! by default everything is private ! private ! declare public subroutines ! public :: max_prec public :: quadratic, quadratic_normalized public :: cubic, cubic_normalized public :: quartic public :: tridiag, pentadiag, invert !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== ! ! function QUADRATIC: ! ------------------ ! ! Function finds real roots of a quadratic equation: ! ! f(x) = a₂ x² + a₁ x + a₀ = 0 ! ! Remark: ! ! The coefficients are renormalized in order to avoid overflow and ! unnecessary underflow. The catastrophic cancellation is avoided as well. ! ! Arguments: ! ! a₂, a₁, a₀ - the quadratic equation coefficients; ! x - the root array; if there are two roots, x(1) corresponds to ! the one with '-' and x(2) to the one with '+'; ! ! Return value: ! ! The number of roots found. ! !=============================================================================== ! function quadratic(a, x) result(nr) ! local variables are not implicit by default ! implicit none ! input/output arguments ! real(kind=8), dimension(3), intent(in) :: a real(kind=8), dimension(2), intent(out) :: x integer :: nr ! local variables ! real(kind=8), dimension(3) :: b real(kind=8) :: bh, dl, dr, tm ! !------------------------------------------------------------------------------- ! ! in order to avoid overflow or needless underflow, divide all coefficients by ! the maximum among them, but only if the maximum value is really large ! b(1:3) = a(1:3) / maxval(abs(a(1:3))) ! check if the coefficient a₂ is not zero ! if (abs(b(3)) > 0.0d+00) then ! check if the coefficient a₀ is not zero ! if (abs(b(1)) > 0.0d+00) then ! coefficients a₂ and a₀ are nonzero, so we solve the quadratic equation ! ! a₂ x² + a₁ x + a₀ = 0 ! ! calculate half of a₁ ! bh = 0.5d+00 * b(2) ! calculate Δ = a₁² - 4 a₂ a₀ ! dl = bh * bh - b(3) * b(1) ! check if Δ is larger then zero ! if (dl > 0.0d+00) then ! calculate quare root of Δ ! dr = sqrt(dl) ! Δ > 0, so the quadratic has two real roots ! if (b(2) > 0.0d+00) then tm = - bh - dr x(1) = tm / b(3) x(2) = b(1) / tm else if (b(2) < 0.0d+00) then tm = - bh + dr x(1) = b(1) / tm x(2) = tm / b(3) else x(2) = dr / b(3) x(1) = - x(2) end if ! update the number of roots ! nr = 2 else if (dl < 0.0d+00) then ! Δ < 0 ! Δ < 0, so the quadratic does not have any real roots ! x(1) = 0.0d+00 x(2) = 0.0d+00 ! update the number of roots ! nr = 0 else ! Δ = 0 ! Δ = 0, so the quadratic has two identical real roots ! x(1) = - 0.5d+00 * b(2) / b(3) x(2) = x(1) ! update the number of roots ! nr = 2 end if ! Δ < 0 else ! a₀ = 0 ! since a₀ is zero, we have one zero root and the second root is calculated from ! the linear formula ! ! (a₂ x + a₁) x = 0 ! tm = - b(2) / b(3) if (tm >= 0.0d+00) then x(1) = 0.0d+00 x(2) = tm else x(1) = tm x(2) = 0.0d+00 end if ! update the number of roots ! nr = 2 end if ! a₀ = 0 else ! a₂ = 0 ! since coefficient a₂ is zero, the quadratic equation is reduced to a linear ! one; one root is undetermined in this case ! ! a₁ x + a₀ = 0 ! x(2) = 0.0d+00 ! find the unique root ! if (abs(b(2)) > 0.0d+00) then ! the coefficient a₁ is not zero, so the quadratic has only one root ! x(1) = - b(1) / b(2) ! update the number of roots ! nr = 1 else ! a₁ = 0 ! the coefficient a₁ is zero, so the quadratic has no roots ! ! a₀ = 0 ! x(1) = 0.0d+00 ! update the number of roots ! nr = 0 end if ! a₁ = 0 end if ! a₂ = 0 !------------------------------------------------------------------------------- ! end function quadratic ! !=============================================================================== ! ! function QUADRATIC_NORMALIZED: ! ----------------------------- ! ! Function finds real roots of the normalized quadratic equation: ! ! f(x) = x² + a₁ x + a₀ = 0 ! ! Remark: ! ! The coefficients are renormalized in order to avoid overflow and ! unnecessary underflow. The catastrophic cancellation is avoided as well. ! ! Arguments: ! ! a₁, a₀ - the quadratic equation coefficients; ! x - the root array; ! ! Return value: ! ! The number of roots found. ! !=============================================================================== ! function quadratic_normalized(a, x) result(nr) ! local variables are not implicit by default ! implicit none ! input/output arguments ! real(kind=8), dimension(2), intent(in) :: a real(kind=8), dimension(2), intent(out) :: x integer :: nr ! local variables ! real(kind=8) :: bh, dl, dr, tm ! !------------------------------------------------------------------------------- ! ! check if the coefficient a₀ is not zero ! if (abs(a(1)) > 0.0d+00) then ! coefficients a₀ is nonzero, so we solve the quadratic equation ! ! x² + a₁ x + a₀ = 0 ! ! calculate half of a₁ ! bh = 0.5d+00 * a(2) ! calculate Δ = ¼ a₁² - a₀ ! dl = bh * bh - a(1) ! check if Δ is larger then zero ! if (dl > 0.0d+00) then ! calculate quare root of Δ ! dr = sqrt(dl) ! Δ > 0, so the quadratic has two real roots, x₁ = - ½ a₁ ± √Δ ! if (a(2) > 0.0d+00) then tm = bh + dr x(1) = - a(1) / tm x(2) = - tm else if (a(2) < 0.0d+00) then tm = bh - dr x(1) = - tm x(2) = - a(1) / tm else x(1) = dr x(2) = - dr end if ! update the number of roots ! nr = 2 else if (dl < 0.0d+00) then ! Δ < 0, so the quadratic does not have any real root ! x(1) = 0.0d+00 x(2) = 0.0d+00 ! update the number of roots ! nr = 0 else ! Δ = 0, so the quadratic has two identical real roots ! x(:) = - bh ! update the number of roots ! nr = 2 end if else ! since a₀ is zero, we have one zero root and the second root is calculated from ! the linear formula ! ! (x + a₁) x = 0 ! x(1) = 0.0d+00 x(2) = - a(2) ! update the number of roots ! nr = 2 end if ! a₀ = 0 !------------------------------------------------------------------------------- ! end function quadratic_normalized ! !=============================================================================== ! ! function CUBIC: ! -------------- ! ! Function finds real roots of a cubic equation: ! ! f(x) = a₃ x³ + a₂ x² + a₁ x + a₀ = 0 ! ! Remark: ! ! The coefficients are renormalized in order to avoid overflow and ! unnecessary underflow. ! ! Arguments: ! ! a(:) - the cubic equation coefficients vector (in increasing moments); ! x(:) - the root vector; ! ! Return value: ! ! The number of roots found. ! !=============================================================================== ! function cubic(a, x) result(nr) ! include external constants ! use constants, only : pi2 ! local variables are not implicit by default ! implicit none ! input/output arguments ! real(kind=8), dimension(4), intent(in) :: a real(kind=8), dimension(3), intent(out) :: x integer :: nr ! local variables ! real(kind=8) :: ac, bc, qc, cc, dc, ec, rc, th, cs, a2, b3 ! local parameters ! real(kind=8), parameter :: c0 = 1.0d+00 / 3.0d+00 ! !------------------------------------------------------------------------------- ! ! check if the coefficient a₃ is not zero ! if (abs(a(4)) > 0.0d+00) then ! check if the coefficient a₀ is not zero ! if (abs(a(1)) > 0.0d+00) then ! both coefficients a₃ and a₀ are not zero, so we have to solve full cubic ! equation ! ! a₃ x³ + a₂ x² + a₁ x + a₀ = 0 ! ! compute coefficient A = 2 a₂³ - 9 a₃ a₂ a₁ + 27 a₃² a₀ ! ac = 2.0d+00 * a(3)**3 - 9.0d+00 * a(4) * a(3) * a(2) & + 27.0d+00 * a(4)**2 * a(1) ! compute coefficient B = a₂² - 3 a₃ a₁ ! bc = a(3)**2 - 3.0d+00 * a(4) * a(2) ! calculate A² and B³ ! a2 = ac**2 b3 = bc**3 ! compute coefficient Q² = A² - 4 B³ ! qc = a2 - 4.0d+00 * b3 ! check the condition for any root existance ! if (qc > 0.0d+00) then ! the case Q > 0, in which only one root is real ! ! calculate C³ = ½ (Q + A) ! cc = 0.5d+00 * (sqrt(qc) + ac) dc = sign(1.0d+00, cc) * abs(cc)**c0 if (abs(cc) > 0.0d+00) then ec = b3 / cc else ec = b3 / a2 ec = 2.0d+00 * ec * (1.0d+00 / sqrt(1.0d+00 - 4.0d+00 * ec) & + 1.0d+00 / ac) end if ec = sign(1.0d+00, ec) * abs(ec)**c0 ! calculate the unique real root ! x(1) = - c0 * (a(3) + dc + ec) / a(4) x(2) = 0.0d+00 x(3) = 0.0d+00 ! update the number of roots ! nr = 1 else if (qc < 0.0d+00) then ! Q < 0, so there are three distinct and real roots ! ! compute Q and R ! th = 3.0d+00 * a(4) qc = bc / th**2 rc = 0.5d+00 * ac / th**3 ! compute cosine ! cs = sign(0.5d+00, ac * bc) * sqrt(ac**2 / abs(bc)**3) ! check condition ! if (cs <= 1.0d+00) then ! there are three real roots; prepare coefficients to calculate them ! th = acos(cs) ac = - 2.0d+00 * sqrt(qc) bc = - c0 * (a(3) / a(4)) ! calculate roots ! x(1) = ac * cos(c0 * th ) + bc x(2) = ac * cos(c0 * (th + pi2)) + bc x(3) = ac * cos(c0 * (th - pi2)) + bc ! update the number of roots ! nr = 3 else ! there is only one real root; prepare coefficients to calculate it ! ac = - sign(1.0d+00, rc) * (abs(rc) + sqrt(rc**2 - qc**3))**c0 if (abs(ac) > 0.0d+00) then bc = qc / ac else bc = 0.0d+00 end if ! calculate the root ! x(1) = (ac + bc) - c0 * (a(3) / a(4)) ! reset remaining roots ! x(2) = 0.0d+00 x(3) = 0.0d+00 ! update the number of roots ! nr = 1 end if else if (abs(bc) > 0.0d+00) then ! the case Q = 0 and B ≠ 0, in which all three roots are real and one root ! is double ! x(1) = 0.5d+00 * (9.0d+00 * a(4) * a(1) - a(3) * a(2)) / bc x(2) = x(1) x(3) = (- 9.0d+00 * a(4) * a(1) + 4.0d+00 * a(3) * a(2) & - a(3)**3 / a(4)) / bc ! update the number of roots ! nr = 3 else ! the case Q = 0 and B = 0, in which all three roots are real and equal ! x(:) = - c0 * a(3) / a(4) ! update the number of roots ! nr = 3 end if end if else ! since the coefficient a₀ is zero, there is one zero root, and the remaining ! roots are found from a quadratic equation ! ! (a₃ x² + a₂ x + a₁) x = 0 ! x(3) = 0.0d+00 ! call function quadratic() to find the remaining roots ! nr = quadratic(a(2:4), x(1:2)) ! increase the number of roots ! nr = nr + 1 end if else ! since the coefficient a₃ is zero, the equation reduces to a quadratic one, ! and one root is undetermined ! ! a₂ x² + a₁ x + a₀ = 0 ! x(3) = 0.0d+00 ! call function quadratic() to find the remaining roots ! nr = quadratic(a(1:3), x(1:2)) end if ! sort roots in the increasing order ! call sort(nr, x(:)) !------------------------------------------------------------------------------- ! end function cubic ! !=============================================================================== ! ! function CUBIC_NORMALIZED: ! ------------------------- ! ! Function finds real roots of the normalized cubic equation: ! ! f(x) = x³ + a₂ x² + a₁ x + a₀ = 0 ! ! Arguments: ! ! a(:) - the cubic equation coefficients vector (in increasing moments); ! x(:) - the root vector; ! ! Return value: ! ! The number of roots found. ! !=============================================================================== ! function cubic_normalized(a, x) result(nr) ! include external constants ! use constants, only : pi2 ! local variables are not implicit by default ! implicit none ! input/output arguments ! real(kind=8), dimension(3), intent(in) :: a real(kind=8), dimension(3), intent(out) :: x integer :: nr ! local variables ! real(kind=8) :: ac, bc, qc, a3, b2, aa, bb, sq, cs, ph ! local parameters ! real(kind=8), parameter :: c0 = 1.0d+00 / 3.0d+00, c1 = 1.0d+00 / 2.7d+01 real(kind=8), parameter :: c2 = 1.0d+00 / 4.0d+00 ! !------------------------------------------------------------------------------- ! ! check if the coefficient a₀ is not zero ! if (abs(a(1)) > 0.0d+00) then ! both coefficients a₃ and a₀ are not zero, so we have to solve full cubic ! equation ! ! x³ + a₂ x² + a₁ x + a₀ = 0 ! ! compute coefficient A = a₂² - 3 a₁ ! ac = c0 * (a(3)**2 - 3.0d+00 * a(2)) ! compute coefficient B = (2 a₂³ - 9 a₂ a₁ + 27 a₀) / 27 ! bc = c1 * (2.0d+00 * a(3)**3 - 9.0d+00 * a(3) * a(2) + 2.7d+01 * a(1)) ! calculate powers of A and B ! a3 = c1 * ac**3 b2 = c2 * bc**2 ! compute coefficient Q² = B² / 4 + A³ / 27 ! qc = b2 - a3 ! check the condition for any root existance ! if (qc > 0.0d+00) then ! the case Q > 0, in which only one root is real ! ! calculate A and B ! sq = sqrt(qc) ph = - 0.5d+00 * bc + sq aa = sign(1.0d+00, ph) * abs(ph)**c0 ph = - 0.5d+00 * bc - sq bb = sign(1.0d+00, ph) * abs(ph)**c0 ! calculate the unique real root ! x(1) = aa + bb - (c0 * a(3)) x(2) = 0.0d+00 x(3) = 0.0d+00 ! update the number of roots ! nr = 1 else if (qc < 0.0d+00) then ! calculate angle φ ! if (bc > 0.0d+00) then cs = - sqrt(b2 / a3) else if (bc < 0.0d+00) then cs = sqrt(b2 / a3) else ! B = 0 cs = 0.0d+00 end if ! there are three real roots; prepare coefficients to calculate them ! ph = acos(cs) aa = 2.0d+00 * sqrt(c0 * ac) ! calculate all three roots ! x(1) = aa * cos(c0 * (ph - pi2)) x(2) = aa * cos(c0 * ph ) x(3) = aa * cos(c0 * (ph + pi2)) ! return to the original variables ! x(:) = x(:) - (c0 * a(3)) ! update the number of roots ! nr = 3 else ! Q = 0 ! the case Q = 0, for which all three roots are real and two roots are equal ! if (bc > 0.0d+00) then aa = sqrt(c0 * ac) x(1) = - 2.0d+00 * aa x(2) = aa x(3) = aa else if (bc < 0.0d+00) then aa = sqrt(c0 * ac) x(1) = 2.0d+00 * aa x(2) = - aa x(3) = - aa else x(:) = 0.0d+00 end if ! return to the original variables ! x(:) = x(:) - (c0 * a(3)) ! update the number of roots ! nr = 3 end if ! Q = 0 else ! a₀ = 0 ! since the coefficient a₀ is zero, there is one zero root, and the remaining ! roots are found by solving the quadratic equation ! ! (x² + a₂ x + a₁) x = 0 ! x(3) = 0.0d+00 ! call function quadratic() to find the remaining roots ! nr = quadratic_normalized(a(2:3), x(1:2)) ! increase the number of roots ! nr = nr + 1 end if ! a₀ = 0 ! sort roots in the increasing order ! call sort(nr, x(:)) !------------------------------------------------------------------------------- ! end function cubic_normalized ! !=============================================================================== ! ! subroutine QUARTIC: ! ------------------ ! ! Subroutine finds real roots of a quartic equation: ! ! f(x) = a₄ * x⁴ + a₃ * x³ + a₂ * x² + a₁ * x + a₀ = 0 ! ! Arguments: ! ! a(:) - the quartic equation coefficients vector (in increasing moments); ! x(:) - the root vector; ! ! References: ! ! Hacke, Jr. J. E., "A Simple Solution of the General Quartic", ! The American Mathematical Monthly, 1941, vol. 48, no. 5, pp. 327-328 ! !=============================================================================== ! function quartic(a, x) result(nr) ! local variables are not implicit by default ! implicit none ! input/output arguments ! real(kind=8), dimension(5), intent(in) :: a real(kind=8), dimension(4), intent(out) :: x integer :: nr ! local variables ! integer :: n, m real(kind=8) :: pc, qc, rc real(kind=8) :: zm, k2, l2, kc, lc ! local arrays ! real(kind=8), dimension(4) :: b real(kind=8), dimension(2) :: y real(kind=8), dimension(3) :: c, z ! !------------------------------------------------------------------------------- ! ! check if the coefficient a₄ is not zero ! if (abs(a(5)) > 0.0d+00) then ! check if the coefficient a₀ is not zero ! if (abs(a(1)) > 0.0d+00) then ! both coefficients a₄ and a₀ are not zero, so we have to solve full quartic ! equation ! ! a₄ * x⁴ + a₃ * x³ + a₂ * x² + a₁ * x + a₀ = 0 ! ! calculate coefficients of P, Q, and R ! pc = (- 3.0d+00 * a(4)**2 + 8.0d+00 * a(5) * a(3)) / (8.0d+00 * a(5)**2) qc = (a(4)**3 - 4.0d+00 * a(5) * a(4) * a(3) & + 8.0d+00 * a(5)**2 * a(2)) / (8.0d+00 * a(5)**3) rc = (- 3.0d+00 * a(4)**4 + 16.0d+00 * a(5) * a(4)**2 * a(3) & - 64.0d+00 * a(5)**2 * a(4) * a(2) + 256.0d+00 * a(5)**3 * a(1)) & / (256.0d+00 * a(5)**4) ! prepare coefficients for the cubic equation ! b(4) = 8.0d+00 b(3) = - 4.0d+00 * pc b(2) = - 8.0d+00 * rc b(1) = 4.0d+00 * pc * rc - qc**2 ! find roots of the cubic equation ! nr = cubic(b(1:4), z(1:3)) ! check if the cubic solver returned any roots ! if (nr > 0) then ! take the maximum root ! zm = z(1) ! calculate coefficients to find quartic roots ! k2 = 2.0d+00 * zm - pc l2 = zm * zm - rc ! check if both coefficients are larger then zero ! if (k2 > 0.0d+00 .or. l2 > 0.0d+00) then ! prepare coefficients K and L ! if (k2 > l2) then kc = sqrt(k2) lc = 0.5d+00 * qc / kc else if (l2 > k2) then lc = sqrt(l2) kc = 0.5d+00 * qc / lc else kc = sqrt(k2) lc = sqrt(l2) end if ! prepare coefficients for the first quadratic equation ! c(3) = 1.0d+00 c(2) = - kc c(1) = zm + lc ! find the first pair of roots ! m = quadratic(c(1:3), y(1:2)) ! update the roots ! do n = 1, m x(n) = y(n) - 0.25d+00 * a(4) / a(5) end do ! update the number of roots ! nr = m ! prepare coefficients for the second quadratic equation ! c(3) = 1.0d+00 c(2) = kc c(1) = zm - lc ! find the second pair of roots ! m = quadratic(c(1:3), y(1:2)) ! update the roots ! do n = 1, m x(nr + n) = y(n) - 0.25d+00 * a(4) / a(5) end do ! update the number of roots ! nr = nr + m ! reset the remainings roots ! do n = nr + 1, 4 x(n) = 0.0d+00 end do else ! both K² and L² are negative, which signifies that the quartic equation has ! no real roots; reset them ! x(:) = 0.0d+00 ! update the number of roots ! nr = 0 end if else ! weird, no cubic roots... it might mean that this is not real quartic equation ! or something is wrong; reset the roots ! x(:) = 0.0d+00 end if else ! since the coefficient a₀ is zero, there is one zero root, and the remaining ! roots are found from a cubic equation ! ! (a₄ * x³ + a₃ * x² + a₂ * x + a₁) * x = 0 ! x(4) = 0.0d+00 ! call function cubic() to find the remaining roots ! nr = cubic(a(2:5), x(1:3)) ! increase the number of roots ! nr = nr + 1 end if else ! since the coefficient a₄ is zero, the equation reduces to a cubic one, and ! one root is undetermined ! ! a₃ * x³ + a₂ * x² + a₁ * x + a₀ = 0 ! x(4) = 0.0d+00 ! call function cubic() to find the remaining roots ! nr = cubic(a(1:4), x(1:3)) end if ! sort roots in the increasing order ! call sort(nr, x(:)) !------------------------------------------------------------------------------- ! end function quartic ! !=============================================================================== ! ! subroutine SORT: ! --------------- ! ! Subroutine sorts the input vector in the ascending order by straight ! insertion. ! ! Arguments: ! ! n - the number of elements to sort; ! x(:) - the vector which needs to be sorted; ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine sort(n, x) ! local variables are not implicit by default ! implicit none ! input/output arguments ! integer , intent(in) :: n real(kind=8), dimension(:), intent(inout) :: x ! local variables ! integer :: k, l real(kind=8) :: y ! !------------------------------------------------------------------------------- ! ! return if only one element ! if (n < 2) return ! iterate over all elements ! do k = 2, n ! pick the element ! y = x(k) ! go back to find the right place ! do l = k - 1, 1, -1 if (x(l) <= y) goto 10 x(l+1) = x(l) end do ! reset the index ! l = 0 ! insert the element ! 10 x(l+1) = y end do ! k = 2, n !------------------------------------------------------------------------------- ! end subroutine sort ! !=============================================================================== ! ! subroutine TRIDIAG: ! ------------------ ! ! Subroutine solves the system of tridiagonal linear equations using ! Thomas algorithm. ! ! Arguments: ! ! n - the size of the matrix; ! 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] https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm ! !=============================================================================== ! subroutine tridiag(n, a, b, c, d, x) implicit none integer , intent(in) :: n real(kind=8), dimension(n), intent(in) :: a, b, c, d real(kind=8), dimension(n), intent(out) :: x integer :: i real(kind=8) :: t real(kind=8), dimension(n) :: cp, dp ! !------------------------------------------------------------------------------- ! ! decomposition and forward substitution ! cp(1) = c(1) / b(1) dp(1) = d(1) / b(1) do i = 2, n 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 ! backward substitution ! x(n) = dp(n) do i = n - 1, 1, -1 x(i) = dp(i) - cp(i) * x(i+1) end do !------------------------------------------------------------------------------- ! end subroutine tridiag ! !=============================================================================== ! ! subroutine PENTADIAG: ! -------------------- ! ! Subroutine solves the system of pentadiagonal linear equations using ! PTRANS-I algorithm. ! ! Arguments: ! ! n - the size of the matrix; ! e(:) - the vector of the outer lower off-diagonal values; ! c(:) - the vector of the inner lower off-diagonal values; ! d(:) - the vector of the diagonal values; ! a(:) - the vector of the inner upper off-diagonal values; ! b(:) - the vector of the outer upper off-diagonal values; ! y(:) - the vector of the right side values; ! x(:) - the solution vector; ! ! References: ! ! [1] Askar, S. S., and Karawia, A. A., ! "On Solving Pentadiagonal Linear Systems via Transformations", ! Mathematical Problems in Engineering, ! 2015, Volume 2015, Article ID 232456, 9 pages ! http://dx.doi.org/10.1155/2015/232456 ! !=============================================================================== ! subroutine pentadiag(n, e, c, d, a, b, y, x) implicit none integer , intent(in) :: n real(kind=8), dimension(n), intent(in) :: e, c, d, a, b, y real(kind=8), dimension(n), intent(out) :: x integer :: i real(kind=8), dimension(n) :: al, bt, gm, mu, z ! !------------------------------------------------------------------------------- ! mu(1) = d(1) al(1) = a(1) / mu(1) bt(1) = b(1) / mu(1) z (1) = y(1) / mu(1) gm(2) = c(2) mu(2) = d(2) - al(1) * gm(2) al(2) = (a(2) - bt(1) * gm(2)) / mu(2) bt(2) = b(2) / mu(2) z (2) = (y(2) - z(1) * gm(2)) / mu(2) do i = 3, n-2 gm(i) = c(i) - al(i-2) * e(i) mu(i) = d(i) - bt(i-2) * e(i) - al(i-1) * gm(i) al(i) = (a(i) - bt(i-1) * gm(i)) / mu(i) bt(i) = b(i) / mu(i) z (i) = (y(i) - z(i-2) * e(i) - z(i-1) * gm(i)) / mu(i) end do gm(n-1) = c(n-1) - al(n-3) * e(n-1) mu(n-1) = d(n-1) - bt(n-3) * e(n-1) - al(n-2) * gm(n-1) al(n-1) = (a(n-1) - bt(n-2) * gm(n-1)) / mu(n-1) gm(n) = c(n) - al(n-2) * e(n) mu(n) = d(n) - bt(n-3) * e(n) - al(n-1) * gm(n) z(n-1) = (y(n-1) - z(n-2) * e(n-1) - z(n-2) * gm(n-1)) / mu(n-1) z(n ) = (y(n) - z(n-1) * e(n) - z(n-1) * gm(n)) / mu(n) x(n ) = z(n) x(n-1) = z(n-1) - al(n-1) * x(n) do i = n - 2, 1, -1 x(i) = z(i) - al(i) * x(i+1) - bt(i) * x(i+2) end do !------------------------------------------------------------------------------- ! end subroutine pentadiag ! !=============================================================================== ! ! subroutine INVERT: ! ----------------- ! ! Subroutine inverts the squared matrix. ! ! Arguments: ! ! n - the size of the matrix; ! m(:,:) - the matrix elements; ! r(:,:) - the inverted matrix elements; ! f - the solution flag; ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine invert(n, m, r, f) ! local variables are not implicit by default ! implicit none ! input/output arguments ! integer , intent(in) :: n real(kind=max_prec), dimension(n,n), intent(in) :: m real(kind=max_prec), dimension(n,n), intent(out) :: r logical , intent(out) :: f ! local variables ! logical :: flag = .true. integer :: i, j, k real(kind=max_prec) :: t real(kind=max_prec), dimension(n,2*n) :: g ! !------------------------------------------------------------------------------- ! ! augment input matrix with an identity matrix ! do i = 1, n do j = 1, 2 * n if (j <= n) then g(i,j) = m(i,j) else if ((i+n) == j) then g(i,j) = 1.0_max_prec else g(i,j) = 0.0_max_prec end if end do end do ! reduce augmented matrix to upper traingular form ! do k = 1, n - 1 if (.not. (abs(g(k,k)) > 0.0_max_prec)) then flag = .false. do i = k + 1, n if (abs(g(i,k)) > 0.0_max_prec) then do j = 1, 2 * n g(k,j) = g(k,j) + g(i,j) end do flag = .true. exit end if if (flag .eqv. .false.) then r(:,:) = 0.0_max_prec f = .false. return end if end do end if do j = k + 1, n t = g(j,k) / g(k,k) do i = k, 2 * n g(j,i) = g(j,i) - t * g(k,i) end do end do end do ! test for invertibility ! do i = 1, n if (.not. (abs(g(i,i)) > 0.0_max_prec)) then r(:,:) = 0.0_max_prec f = .false. return end if end do ! make diagonal elements as 1 ! do i = 1, n t = g(i,i) do j = i , 2 * n g(i,j) = g(i,j) / t end do end do ! reduced right side half of augmented matrix to identity matrix ! do k = n - 1, 1, -1 do i = 1, k t = g(i,k+1) do j = k, 2 * n g(i,j) = g(i,j) - g(k+1,j) * t end do end do end do ! store answer ! do i = 1, n do j = 1, n r(i,j) = g(i,j+n) end do end do f = .true. !------------------------------------------------------------------------------- ! end subroutine invert !=============================================================================== ! end module algebra