ALGEBRA: Add cubic and quartic solvers.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
b8fed51178
commit
1930bda0cc
720
src/algebra.F90
720
src/algebra.F90
@ -40,6 +40,8 @@ module algebra
|
||||
! declare public subroutines
|
||||
!
|
||||
public :: quadratic, quadratic_normalized
|
||||
public :: cubic, cubic_normalized
|
||||
public :: quartic
|
||||
|
||||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||
!
|
||||
@ -360,6 +362,724 @@ module algebra
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
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, eps = epsilon(1.0d+00)
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! check if the coefficient a₃ is not zero
|
||||
!
|
||||
if (a(4) /= 0.0d+00) then
|
||||
|
||||
! check if the coefficient a₀ is not zero
|
||||
!
|
||||
if (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 (cc == 0.0d+00) then
|
||||
ec = b3 / a2
|
||||
ec = 2.0d+00 * ec * (1.0d+00 / sqrt(1.0d+00 - 4.0d+00 * ec) &
|
||||
+ 1.0d+00 / ac)
|
||||
else
|
||||
ec = b3 / cc
|
||||
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
|
||||
|
||||
if (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
|
||||
|
||||
else
|
||||
|
||||
! 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 (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
|
||||
|
||||
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 (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 (a(5) /= 0.0d+00) then
|
||||
|
||||
! check if the coefficient a₀ is not zero
|
||||
!
|
||||
if (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
|
||||
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
@ -124,7 +124,7 @@ clean-all: clean-bak clean-data clean-exec clean-logs clean-modules \
|
||||
|
||||
#-------------------------------------------------------------------------------
|
||||
|
||||
algebra.o : algebra.F90
|
||||
algebra.o : algebra.F90 constants.o
|
||||
blocks.o : blocks.F90 error.o timers.o
|
||||
boundaries.o : boundaries.F90 blocks.o coordinates.o equations.o error.o \
|
||||
interpolations.o mpitools.o timers.o
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user