1364 lines
31 KiB
Fortran
1364 lines
31 KiB
Fortran
!!******************************************************************************
|
|
!!
|
|
!! 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-2023 Grzegorz Kowal <grzegorz@amuncode.org>
|
|
!!
|
|
!! 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 <http://www.gnu.org/licenses/>.
|
|
!!
|
|
!!******************************************************************************
|
|
!!
|
|
!! 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
|