Skip to content

Commit

Permalink
Merge pull request #91 from ipqa-research/test_volume_solvers
Browse files Browse the repository at this point in the history 
test for exact root solvers
  • Loading branch information
@fedebenelli
fedebenelli authored Jul 18, 2024
2 parents 9613363 + 78a9be3 commit 77fa993
Show file tree
Hide file tree
Showing 2 changed files with 102 additions and 12 deletions.
9 changes: 3 additions & 6 deletions src/math/linalg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,21 +75,18 @@ subroutine cubic_roots(parameters, real_roots, complex_roots, flag)
complex_roots = nan

if (abs(disc) < 1e-15) then
!print *, "disc: ", disc
flag = 0
real_roots(1) = 3*q/p
real_roots(2) = -3*q/(2*p)
real_roots(3) = real_roots(2)
elseif (disc < 0) then
!print *, "disc: ", disc
flag = -1
theta = acos(0.5_pr * 3 * q / p * sqrt(-3/p))
real_roots(1) = 2 * sqrt(-p/3) * cos(theta/3)
real_roots(2) = 2 * sqrt(-p/3) * cos((theta + 2*pi)/3)
real_roots(3) = 2 * sqrt(-p/3) * cos((theta + 4*pi)/3)
call sort(real_roots)
elseif (disc > 0) then
!print *, "disc: ", disc
flag = 1
u = ((-q + sqrt(disc))/2)
v = ((-q - sqrt(disc))/2)
Expand Down Expand Up @@ -124,8 +121,8 @@ subroutine cubic_roots_rosendo(parameters, real_roots, complex_roots, flag)
theta = acos(R / sqrt(Q**3))

real_roots(1) = -2 * sqrt(Q) * cos(theta / 3.0_pr) - d1 / 3.0_pr
real_roots(2) = -2 * sqrt(Q) * cos((theta + 2 * acos(-1.0_pr)) / 3.0_pr) - d1 / 3.0_pr
real_roots(3) = -2 * sqrt(Q) * cos((theta - 2 * acos(-1.0_pr)) / 3.0_pr) - d1 / 3.0_pr
real_roots(2) = -2 * sqrt(Q) * cos((theta + 2 * pi) / 3.0_pr) - d1 / 3.0_pr
real_roots(3) = -2 * sqrt(Q) * cos((theta - 2 * pi) / 3.0_pr) - d1 / 3.0_pr

! Correction??
! do i=1,100
Expand All @@ -139,7 +136,7 @@ subroutine cubic_roots_rosendo(parameters, real_roots, complex_roots, flag)
else
A = - sign((abs(R) + sqrt(R**2 - Q**3))**(1.0_pr/3.0_pr), R)

if (-log10(A) > 15.0_pr) then
if (abs(A) < 1e-6) then
A = 0.0_pr
B = 0.0_pr
else
Expand Down
105 changes: 99 additions & 6 deletions test/test_math/test_math.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ subroutine collect_suite(testsuite)
new_unittest("Test dx_to_dn", test_dx_to_dn), &
new_unittest("Test sq_error", test_sq_error), &
new_unittest("Test cardano", test_cardano_method), &
new_unittest("Test rosendo", test_rosendo_method), &
new_unittest("Test newton", test_newton_method) &
]
end subroutine collect_suite
Expand Down Expand Up @@ -68,8 +69,100 @@ subroutine test_cardano_method(error)
call check(error, &
maxval(abs(rr - [-0.454216, 0.8601986, 25.594016])) < 1e-5_pr &
)

! =======================================================================
! 3 Real roots
! x1 = 1, x2 = 3, x3 = 4
! -----------------------------------------------------------------------
p = [1._pr, -8._pr, 19._pr, -12._pr]

call cubic_roots(p, rr, cr, flag)

call check(error, maxval(rr - [1.0_pr, 3.0_pr, 4.0_pr]) < 1e-7)
call check(error, flag == -1)

! =======================================================================
! 1 real root different, 2 equal real roots
! x1 = 1, x2 = x3 = 4
! -----------------------------------------------------------------------
p = [1._pr, -9._pr, 24._pr, -16._pr]

call cubic_roots(p, rr, cr, flag)

call check(error, maxval(rr - [1.0_pr, 4.0_pr, 4.0_pr]) < 1e-7)
call check(error, flag == 0)

! =======================================================================
! 1 real root, 2 complex
! x1 = 1, x2 = i, x3 = -i
! -----------------------------------------------------------------------
p = [1._pr, -1._pr, 1._pr, -1._pr]

call cubic_roots(p, rr, cr, flag)

call check(error, (rr(1) - 1.0_pr) < 1e-7)
call check(error, flag == 1)
end subroutine test_cardano_method

subroutine test_rosendo_method(error)
use yaeos__math_linalg, only: pr, cubic_roots_rosendo

type(error_type), allocatable, intent(out) :: error

real(pr) :: p(4)
real(pr) :: a, b
real(pr) :: roots(3)
complex(pr) :: c_roots(3)
integer :: flag

! =======================================================================
! 3 Real roots
! x1 = 1, x2 = 3, x3 = 4
! -----------------------------------------------------------------------
p = [1._pr, -8._pr, 19._pr, -12._pr]

call cubic_roots_rosendo(p, roots, c_roots, flag)

call check(error, maxval(roots - [1.0_pr, 3.0_pr, 4.0_pr]) < 1e-7)
call check(error, flag == -1)

! =======================================================================
! 1 real root different, 2 equal real roots
! x1 = 1, x2 = x3 = 4
! -----------------------------------------------------------------------
p = [1._pr, -9._pr, 24._pr, -16._pr]

call cubic_roots_rosendo(p, roots, c_roots, flag)

call check(error, maxval(roots - [1.0_pr, 4.0_pr, 4.0_pr]) < 1e-7)
call check(error, flag == -1)

! =======================================================================
! 1 real root, 2 complex
! x1 = 1, x2 = i, x3 = -i
! -----------------------------------------------------------------------
p = [1._pr, -1._pr, 1._pr, -1._pr]

call cubic_roots_rosendo(p, roots, c_roots, flag)

call check(error, maxval(roots - [1.0_pr, 1.0_pr, 1.0_pr]) < 1e-7)
call check(error, flag == 1)

! =======================================================================
! 1 real root and two small complex roots
! x1 = 1, x2 = a+bi, x3 = a-bi
! -----------------------------------------------------------------------
a = 1._pr
b = 1e-6_pr

p = [1._pr, -(2*a + 1), (a**2 + b**2 + 2 * a), -(a**2 + b**2)]

call cubic_roots_rosendo(p, roots, c_roots, flag)

call check(error, maxval(roots - [1.0_pr, 1.0_pr, 1.0_pr]) < 1e-7)
call check(error, flag == 1)
end subroutine test_rosendo_method

subroutine test_newton_method(error)
use yaeos__math, only: pr, newton
type(error_type), allocatable, intent(out) :: error
Expand All @@ -82,11 +175,11 @@ subroutine test_newton_method(error)
call check(error, abs(x - sqrt(2._pr)) < tol)
contains
subroutine foo(xx, f, df)
real(pr), intent(in) :: xx
real(pr), intent(out) :: f
real(pr), intent(out) :: df
f = xx**2 - 2
df = 2*xx
end subroutine foo
real(pr), intent(in) :: xx
real(pr), intent(out) :: f
real(pr), intent(out) :: df
f = xx**2 - 2
df = 2*xx
end subroutine foo
end subroutine test_newton_method
end module test_math

0 comments on commit 77fa993

Please sign in to comment.