Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Compute Gauss-Turán quadratures #252

Draft
wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Compute Gauss-Turán quadratures #252

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8026cde
Set up PR
SouthEndMusic Jul 22, 2024
b60c6a4
Use TaylorDiff
SouthEndMusic Jul 22, 2024
7e0b9ab
Fix using TaylorDiff with DoubleFloats
SouthEndMusic Jul 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 61 additions & 0 deletions GaussTuranExampleTemp/example.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
using Optim
using DoubleFloats
using TaylorDiff
using PreallocationTools
using Integrals

# This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201
include("linear_algebra_fix.jl")

FT = Double64
n = 5
s = 1
a = FT(0.0)
b = FT(1.0)

"""
Example functions ϕⱼ from [1]
"""
function f(x::T, j::Int)::T where {T <: Number}
pow = if j % 2 == 0
j / 2 - 1
else
(j - 1) / 2 - 1 / 3
end
x^pow
end

I = Integrals.GaussTuran(
f,
a,
b,
n,
s;
optimization_options = Optim.Options(;
x_abstol = FT(1e-250),
g_tol = FT(1e-250),
show_trace = true,
show_every = 100,
iterations = 10_000
)
)

# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ
max_int_error_upper = maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j])
for j in 1:(I.cache.N))
@show max_int_error_upper
# max_int_error_upper = 2.465190328815662e-32

# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ
max_int_error_lower = maximum(
abs(I(x -> f(x, j)) - I.cache.rhs_lower[j - I.cache.N])
for
j in (I.cache.N + 1):(I.cache.N + I.cache.n)
)
@show max_int_error_lower
# max_int_error_lower = 9.079315240327527e-30

# Example with eˣ
exp_error = abs(I(Base.exp) - (Base.exp(1) - 1))
@show exp_error;
# exp_error = 6.43662141695295e-18
296 changes: 296 additions & 0 deletions GaussTuranExampleTemp/linear_algebra_fix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,296 @@
using LinearAlgebra
using LinearAlgebra: AdjOrTrans, require_one_based_indexing, _ustrip

# manually hoisting b[j] significantly improves performance as of Dec 2015
# manually eliding bounds checking significantly improves performance as of Dec 2015
# replacing repeated references to A.data[j,j] with [Ajj = A.data[j,j] and references to Ajj]
# does not significantly impact performance as of Dec 2015
# in the transpose and conjugate transpose naive substitution variants,
# accumulating in z rather than b[j,k] significantly improves performance as of Dec 2015
function LinearAlgebra.generic_trimatdiv!(
C::AbstractVecOrMat, uploc, isunitc, tfun::Function,
A::AbstractMatrix, B::AbstractVecOrMat)
require_one_based_indexing(C, A, B)
mA, nA = size(A)
m, n = size(B, 1), size(B, 2)
if nA != m
throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal"))
end
if size(C) != size(B)
throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
end
iszero(mA) && return C
oA = oneunit(eltype(A))
@inbounds if uploc == 'U'
if isunitc == 'N'
if tfun === identity
for k in 1:n
amm = A[m, m]
iszero(amm) && throw(SingularException(m))
Cm = C[m, k] = amm \ B[m, k]
# fill C-column
for i in (m - 1):-1:1
C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm
end
for j in (m - 1):-1:1
ajj = A[j, j]
iszero(ajj) && throw(SingularException(j))
Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
for i in (j - 1):-1:1
C[i, k] -= _ustrip(A[i, j]) * Cj
end
end
end
else # tfun in (adjoint, transpose)
for k in 1:n
for j in 1:m
ajj = A[j, j]
iszero(ajj) && throw(SingularException(j))
Bj = B[j, k]
for i in 1:(j - 1)
Bj -= tfun(A[i, j]) * C[i, k]
end
C[j, k] = tfun(ajj) \ Bj
end
end
end
else # isunitc == 'U'
if tfun === identity
for k in 1:n
Cm = C[m, k] = oA \ B[m, k]
# fill C-column
for i in (m - 1):-1:1
C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm
end
for j in (m - 1):-1:1
Cj = C[j, k]
for i in 1:(j - 1)
C[i, k] -= _ustrip(A[i, j]) * Cj
end
end
end
else # tfun in (adjoint, transpose)
for k in 1:n
for j in 1:m
Bj = B[j, k]
for i in 1:(j - 1)
Bj -= tfun(A[i, j]) * C[i, k]
end
C[j, k] = oA \ Bj
end
end
end
end
else # uploc == 'L'
if isunitc == 'N'
if tfun === identity
for k in 1:n
a11 = A[1, 1]
iszero(a11) && throw(SingularException(1))
C1 = C[1, k] = a11 \ B[1, k]
# fill C-column
for i in 2:m
C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1
end
for j in 2:m
ajj = A[j, j]
iszero(ajj) && throw(SingularException(j))
Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
for i in (j + 1):m
C[i, k] -= _ustrip(A[i, j]) * Cj
end
end
end
else # tfun in (adjoint, transpose)
for k in 1:n
for j in m:-1:1
ajj = A[j, j]
iszero(ajj) && throw(SingularException(j))
Bj = B[j, k]
for i in (j + 1):m
Bj -= tfun(A[i, j]) * C[i, k]
end
C[j, k] = tfun(ajj) \ Bj
end
end
end
else # isunitc == 'U'
if tfun === identity
for k in 1:n
C1 = C[1, k] = oA \ B[1, k]
# fill C-column
for i in 2:m
C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1
end
for j in 2:m
Cj = C[j, k]
for i in (j + 1):m
C[i, k] -= _ustrip(A[i, j]) * Cj
end
end
end
else # tfun in (adjoint, transpose)
for k in 1:n
for j in m:-1:1
Bj = B[j, k]
for i in (j + 1):m
Bj -= tfun(A[i, j]) * C[i, k]
end
C[j, k] = oA \ Bj
end
end
end
end
end
return C
end
# conjugate cases
function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function,
xA::AdjOrTrans, B::AbstractVecOrMat)
A = parent(xA)
require_one_based_indexing(C, A, B)
mA, nA = size(A)
m, n = size(B, 1), size(B, 2)
if nA != m
throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal"))
end
if size(C) != size(B)
throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
end
iszero(mA) && return C
oA = oneunit(eltype(A))
@inbounds if uploc == 'U'
if isunitc == 'N'
for k in 1:n
amm = conj(A[m, m])
iszero(amm) && throw(SingularException(m))
Cm = C[m, k] = amm \ B[m, k]
# fill C-column
for i in (m - 1):-1:1
C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm
end
for j in (m - 1):-1:1
ajj = conj(A[j, j])
iszero(ajj) && throw(SingularException(j))
Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
for i in (j - 1):-1:1
C[i, k] -= _ustrip(conj(A[i, j])) * Cj
end
end
end
else # isunitc == 'U'
for k in 1:n
Cm = C[m, k] = oA \ B[m, k]
# fill C-column
for i in (m - 1):-1:1
C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm
end
for j in (m - 1):-1:1
Cj = C[j, k]
for i in 1:(j - 1)
C[i, k] -= _ustrip(conj(A[i, j])) * Cj
end
end
end
end
else # uploc == 'L'
if isunitc == 'N'
for k in 1:n
a11 = conj(A[1, 1])
iszero(a11) && throw(SingularException(1))
C1 = C[1, k] = a11 \ B[1, k]
# fill C-column
for i in 2:m
C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1
end
for j in 2:m
ajj = conj(A[j, j])
iszero(ajj) && throw(SingularException(j))
Cj = C[j, k] = _ustrip(ajj) \ C[j, k]
for i in (j + 1):m
C[i, k] -= _ustrip(conj(A[i, j])) * Cj
end
end
end
else # isunitc == 'U'
for k in 1:n
C1 = C[1, k] = oA \ B[1, k]
# fill C-column
for i in 2:m
C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1
end
for j in 1:m
Cj = C[j, k]
for i in (j + 1):m
C[i, k] -= _ustrip(conj(A[i, j])) * Cj
end
end
end
end
end
return C
end

function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function,
A::AbstractMatrix, B::AbstractMatrix)
require_one_based_indexing(C, A, B)
m, n = size(A)
if size(B, 1) != n
throw(DimensionMismatch(lazy"right hand side B needs first dimension of size $n, has size $(size(B,1))"))
end
if size(C) != size(A)
throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of left hand side, $(size(A))"))
end
oB = oneunit(eltype(B))
unit = isunitc == 'U'
@inbounds if uploc == 'U'
if tfun === identity
for i in 1:m
for j in 1:n
Aij = A[i, j]
for k in 1:(j - 1)
Aij -= C[i, k] * B[k, j]
end
unit || (iszero(B[j, j]) && throw(SingularException(j)))
C[i, j] = Aij / (unit ? oB : B[j, j])
end
end
else # tfun in (adjoint, transpose)
for i in 1:m
for j in n:-1:1
Aij = A[i, j]
for k in (j + 1):n
Aij -= C[i, k] * tfun(B[j, k])
end
unit || (iszero(B[j, j]) && throw(SingularException(j)))
C[i, j] = Aij / (unit ? oB : tfun(B[j, j]))
end
end
end
else # uploc == 'L'
if tfun === identity
for i in 1:m
for j in n:-1:1
Aij = A[i, j]
for k in (j + 1):n
Aij -= C[i, k] * B[k, j]
end
unit || (iszero(B[j, j]) && throw(SingularException(j)))
C[i, j] = Aij / (unit ? oB : B[j, j])
end
end
else # tfun in (adjoint, transpose)
for i in 1:m
for j in 1:n
Aij = A[i, j]
for k in 1:(j - 1)
Aij -= C[i, k] * tfun(B[j, k])
end
unit || (iszero(B[j, j]) && throw(SingularException(j)))
C[i, j] = Aij / (unit ? oB : tfun(B[j, j]))
end
end
end
end
return C
end
4 changes: 4 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,9 @@ Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[extensions]
Expand All @@ -29,6 +32,7 @@ IntegralsCubaExt = "Cuba"
IntegralsCubatureExt = "Cubature"
IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
IntegralsForwardDiffExt = "ForwardDiff"
IntegralsGaussTuranExt = ["Optim", "TaylorDiff", "PreallocationTools"]
IntegralsMCIntegrationExt = "MCIntegration"
IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]

Expand Down
Loading
Loading