-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Move extraneous utilities to separate file
- Loading branch information
1 parent
5b2bbea
commit 07a6ed6
Showing
4 changed files
with
165 additions
and
162 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,164 @@ | ||
| ## --- Various useful internal utility functions | ||
|
|
||
| # Utility function for agepoint and Tpoint buffers | ||
| function collectto!(buffer, a, b, c) | ||
| i₀ = firstindex(buffer) | ||
| copyto!(buffer, i₀, a, 1, length(a)) | ||
| i₀ += length(a) | ||
| copyto!(buffer, i₀, b, 1, length(b)) | ||
| i₀ += length(b) | ||
| copyto!(buffer, i₀, c, 1, length(c)) | ||
| i₀ += length(c) | ||
| return view(buffer, firstindex(buffer):i₀-1) | ||
| end | ||
|
|
||
| # Utility functions for checking maximum reheating or cooling rate | ||
| function maxdiff(x::AbstractVector) | ||
| i₀ = firstindex(x) | ||
| δₘ = x[i₀+1] - x[i₀] | ||
| if length(x) > 2 | ||
| last = x[i₀+1] | ||
| @inbounds for i ∈ (i₀+2):(i₀+length(x)-1) | ||
| δᵢ = x[i] - last | ||
| if δᵢ > δₘ | ||
| δₘ = δᵢ | ||
| end | ||
| last = x[i] | ||
| end | ||
| end | ||
| return δₘ | ||
| end | ||
| function mindiff(x::AbstractVector) | ||
| i₀ = firstindex(x) | ||
| δₘ = x[i₀+1] - x[i₀] | ||
| if length(x) > 2 | ||
| last = x[i₀+1] | ||
| @inbounds for i ∈ (i₀+2):(i₀+length(x)-1) | ||
| δᵢ = x[i] - last | ||
| if δᵢ < δₘ | ||
| δₘ = δᵢ | ||
| end | ||
| last = x[i] | ||
| end | ||
| end | ||
| return δₘ | ||
| end | ||
| function maxabsdiff(x::AbstractVector) | ||
| i₀ = firstindex(x) | ||
| δₘ = abs(x[i₀+1] - x[i₀]) | ||
| if length(x) > 2 | ||
| last = x[i₀+1] | ||
| @inbounds for i ∈ (i₀+2):(i₀+length(x)-1) | ||
| δᵢ = abs(x[i] - last) | ||
| if δᵢ > δₘ | ||
| δₘ = δᵢ | ||
| end | ||
| last = x[i] | ||
| end | ||
| end | ||
| return δₘ | ||
| end | ||
|
|
||
| # Check if point k is distinct from other points in list within ± δ | ||
| function isdistinct(points::DenseArray, npoints::Int, k::Int, δ::Number) | ||
| @inbounds for i = 1:npoints | ||
| if i!=k && abs(points[i] - points[k]) < δ | ||
| return false | ||
| end | ||
| end | ||
| return true | ||
| end | ||
|
|
||
| function pointsininterval(points::DenseArray, npoints::Int, min::Number, max::Number) | ||
| n = 0 | ||
| @inbounds for i = 1:npoints | ||
| if min < points[i] < max | ||
| n += 1 | ||
| end | ||
| end | ||
| return n | ||
| end | ||
|
|
||
|
|
||
|
|
||
| ## --- Modify LinearAlgebra.lu! to reuse du2 for pivoting | ||
|
|
||
| using LinearAlgebra: BlasInt, checknonsingular | ||
| # Modified from LinearAlgebra stdlib to reuse `du2` and `ipiv` | ||
| function lu!(A::Tridiagonal{T,V}, ipiv::Vector{BlasInt}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T,V} | ||
|
|
||
| # Extract values | ||
| n = size(A, 1) | ||
| @assert length(ipiv) == n | ||
|
|
||
| # initialize variables | ||
| info = 0 | ||
| dl = A.dl | ||
| d = A.d | ||
| du = A.du | ||
| if dl === du | ||
| throw(ArgumentError("off-diagonals of `A` must not alias")) | ||
| end | ||
| # Check if Tridiagonal matrix already has du2 for pivoting | ||
| d2len = max(0, n-2) # Proper length of a second off-diagonal | ||
| has_du2_defined = isdefined(A, :du2) && length(A.du2) == d2len | ||
| du2 = (has_du2_defined ? A.du2 : similar(d, d2len))::V | ||
| fill!(du2, 0) | ||
|
|
||
| @inbounds begin | ||
| for i = 1:n | ||
| ipiv[i] = i | ||
| end | ||
| for i = 1:n-2 | ||
| # pivot or not? | ||
| if pivot === NoPivot() || abs(d[i]) >= abs(dl[i]) | ||
| # No interchange | ||
| if d[i] != 0 | ||
| fact = dl[i]/d[i] | ||
| dl[i] = fact | ||
| d[i+1] -= fact*du[i] | ||
| du2[i] = 0 | ||
| end | ||
| else | ||
| # Interchange | ||
| fact = d[i]/dl[i] | ||
| d[i] = dl[i] | ||
| dl[i] = fact | ||
| tmp = du[i] | ||
| du[i] = d[i+1] | ||
| d[i+1] = tmp - fact*d[i+1] | ||
| du2[i] = du[i+1] | ||
| du[i+1] = -fact*du[i+1] | ||
| ipiv[i] = i+1 | ||
| end | ||
| end | ||
| if n > 1 | ||
| i = n-1 | ||
| if pivot === NoPivot() || abs(d[i]) >= abs(dl[i]) | ||
| if d[i] != 0 | ||
| fact = dl[i]/d[i] | ||
| dl[i] = fact | ||
| d[i+1] -= fact*du[i] | ||
| end | ||
| else | ||
| fact = d[i]/dl[i] | ||
| d[i] = dl[i] | ||
| dl[i] = fact | ||
| tmp = du[i] | ||
| du[i] = d[i+1] | ||
| d[i+1] = tmp - fact*d[i+1] | ||
| ipiv[i] = i+1 | ||
| end | ||
| end | ||
| # check for a zero on the diagonal of U | ||
| for i = 1:n | ||
| if d[i] == 0 | ||
| info = i | ||
| break | ||
| end | ||
| end | ||
| end | ||
| B = has_du2_defined ? A : Tridiagonal{T,V}(dl, d, du, du2) | ||
| check && checknonsingular(info, pivot) | ||
| return LU{T,Tridiagonal{T,V},typeof(ipiv)}(B, ipiv, convert(BlasInt, info)) | ||
| end |