Merge pull request #52 from ctkelley/main

Address JOSS reviewer comments: I

.github/workflows/ci.yml

@@ -26,9 +26,9 @@ jobs:
           - '~1.11.0-0'
           - 'nightly'
         os:
-#          - ubuntu-latest
+          - ubuntu-latest
           - macOS-latest
-#          - windows-latest
+          - windows-latest
         arch:
           - x64
     steps:

Project.toml

@@ -6,12 +6,14 @@ version = "0.1.3"
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 LinearAlgebra = "1"
 OhMyThreads = "0.5, 0.6"
+Reexport = "1.2.2"
 SIAMFANLEquations = "1"
 SparseArrays = "1"
 Test = "1"

README.md

@@ -182,7 +182,6 @@ julia> using MultiPrecisionArrays
 julia> using MultiPrecisionArrays.Examples
 
 julia> using BenchmarkTools
-
 julia> N=4096;
 
 julia> G=Gmat(N);
@@ -347,14 +346,17 @@ There are more examples for this in the [paper](https://github.com/ctkelley/Mult
 
 ## Dependencies
 
-As of now you need to install these packages
+As of now these packages are the only direct dependencies
 
 - OhMyThreads
+- Reexport
 - SIAMFANLEquations
 
-and use LinearAlgebra and SparseArrays from Base. I use the Krylov solvers and examples from SIAMFANLEquations. OhMyThreads is here because 
+as we use LinearAlgebra and SparseArrays from Base. I use the Krylov solvers and examples from SIAMFANLEquations. OhMyThreads is here because 
 threading with ```OhMyThreads.@tasks``` makes the LU for Float16 run much faster. Once LU for Float16 is in LAPACK/BLAS, I will eliminate that dependency.
 
+I use Reexport to make the identity ```I``` part of the Examples submodule. This makes it easier to run the examples in the docs.
+
 ## Endorsement
 
 I have used this in my own work and will be adding links to that stuff as I finish it. 

src/Examples/Gmat.jl

@@ -1,12 +1,12 @@
 """
 Gmat: discrete Green's operator for -d^2/dx^2 in 1D
 """
-function Gmat(n,T=Float64)
-    onet=T(1.0)
-    h = onet / (n-onet)
+function Gmat(n, T = Float64)
+    onet = T(1.0)
+    h = onet / (n - onet)
     Ns = collect(0:1:n-1)
-    X = h*Ns
-    X[n]=onet
+    X = h * Ns
+    X[n] = onet
     G = [greens(x, y, onet) for x in X, y in X]
     G .*= h
     return G

src/Factorizations/hlu!.jl

@@ -25,7 +25,7 @@ function hlu!(A::AbstractMatrix{T}) where {T}
     # Extract values and make sure the problem is square
     m, n = size(A)
     # Small n? Revert to normal lu
-    (n < 128) && (AF=lu!(A); return AF )
+    (n < 128) && (AF = lu!(A); return AF)
     minmn = min(m, n)
     # Initialize variables
     info = 0
@@ -71,13 +71,13 @@ function hlu!(A::AbstractMatrix{T}) where {T}
                 info = k
             end
             # Update the rest
-            ntasks=min(nthreads(), 1 + floor(Int,(n-k)/8))
-                  tforeach(k+1:n; ntasks=ntasks) do j
-                    Akj = -A[k, j]
-                    @inbounds @simd ivdep for i = k+1:m
-                        A[i, j] += A[i, k] * Akj
-                    end # i loop
-                end #j loop
+            ntasks = min(nthreads(), 1 + floor(Int, (n - k) / 8))
+            tforeach(k+1:n; ntasks = ntasks) do j
+                Akj = -A[k, j]
+                @inbounds @simd ivdep for i = k+1:m
+                    A[i, j] += A[i, k] * Akj
+                end # i loop
+            end #j loop
         end
     end
     checknonsingular(info, pivot)

src/Factorizations/mpblu!.jl

@@ -10,17 +10,17 @@ storage for the things BiCGSTAB needs.
 You get a factorization
 object as output and can use ```\\``` to solve linear systems.
 """
-function mpblu!(MPBA::MPBArray; residterm=residtermdefault)
-AL=MPBA.AL
-AH=MPBA.AH
-VStore=MPBA.VStore
-KStore=MPBA.KStore
-res=MPBA.residual
-sol=MPBA.sol
-TF=eltype(AL)
-(TF == Float16) ? ALF = hlu!(AL) : ALF = lu!(AL)
-MPF=MPBFact(AH, AL, ALF, VStore, KStore, res, sol, true, residterm)
-return MPF
+function mpblu!(MPBA::MPBArray; residterm = residtermdefault)
+    AL = MPBA.AL
+    AH = MPBA.AH
+    VStore = MPBA.VStore
+    KStore = MPBA.KStore
+    res = MPBA.residual
+    sol = MPBA.sol
+    TF = eltype(AL)
+    (TF == Float16) ? ALF = hlu!(AL) : ALF = lu!(AL)
+    MPF = MPBFact(AH, AL, ALF, VStore, KStore, res, sol, true, residterm)
+    return MPF
 end
 
 """
@@ -46,21 +46,21 @@ structure is immutable and MPF.AF.info cannot be changed.
 Reassigning MPG works and resuses almost all of the storage in the
 original array.
 """
-function mpblu!(MPG::MPBFact, A::AbstractArray{TW,2}) where TW <: Real
-TF=eltype(MPG.AH)
-(TF == TW) || error("Precision error in mplu!")
-AH=MPG.AH
-AH = A
-TF = eltype(MPG.AL)
-AL=MPG.AL
-AL .= TF.(A)
-(TF == Float16) ? AF = hlu!(AL) : AF = lu!(AL)
-MPG.AF.ipiv .= AF.ipiv
-VStore=MPG.VStore 
-KStore=MPG.KStore
-residterm=MPG.residterm 
-MPG=MPBFact(AH, AL, AF, VStore, KStore, MPG.residual, MPG.sol, true, residterm)
-return MPG
+function mpblu!(MPG::MPBFact, A::AbstractArray{TW,2}) where {TW<:Real}
+    TF = eltype(MPG.AH)
+    (TF == TW) || error("Precision error in mplu!")
+    AH = MPG.AH
+    AH = A
+    TF = eltype(MPG.AL)
+    AL = MPG.AL
+    AL .= TF.(A)
+    (TF == Float16) ? AF = hlu!(AL) : AF = lu!(AL)
+    MPG.AF.ipiv .= AF.ipiv
+    VStore = MPG.VStore
+    KStore = MPG.KStore
+    residterm = MPG.residterm
+    MPG = MPBFact(AH, AL, AF, VStore, KStore, MPG.residual, MPG.sol, true, residterm)
+    return MPG
 end
 
 
@@ -150,10 +150,14 @@ julia> xp=Float64.(A)\\b; norm(xp-mout3.sol,Inf)
 ```
 
 """
-function mpblu(A::AbstractArray{TW,2}; residterm=residtermdefault,
-          TF=Float32, TR=nothing) where TW <: Real
-(TW==Float32) ? TF=Float16 : TF=TF
-MPBA=MPBArray(A; TF=TF, TR=TR)
-MPBF=mpblu!(MPBA; residterm=residterm)
-return MPBF
+function mpblu(
+    A::AbstractArray{TW,2};
+    residterm = residtermdefault,
+    TF = Float32,
+    TR = nothing,
+) where {TW<:Real}
+    (TW == Float32) ? TF = Float16 : TF = TF
+    MPBA = MPBArray(A; TF = TF, TR = TR)
+    MPBF = mpblu!(MPBA; residterm = residterm)
+    return MPBF
 end

src/Factorizations/mpglu!.jl

@@ -15,17 +15,17 @@ The kwarg ```residterm``` sets the termination criterion.
 small residuals.  ```residterm == false``` terminates the iteration on
 small normwise backward errors. Look at the docs for details.
 """
-function mpglu!(MPGA::MPGArray; residterm=residtermdefault)
-AL=MPGA.AL
-AH=MPGA.AH
-VStore=MPGA.VStore
-KStore=MPGA.KStore
-res=MPGA.residual
-sol=MPGA.sol
-TF=eltype(AL)
-(TF == Float16) ? ALF = hlu!(AL) : ALF = lu!(AL)
-MPF=MPGEFact(AH, AL, ALF, VStore, KStore, res, sol, true,residterm)
-return MPF
+function mpglu!(MPGA::MPGArray; residterm = residtermdefault)
+    AL = MPGA.AL
+    AH = MPGA.AH
+    VStore = MPGA.VStore
+    KStore = MPGA.KStore
+    res = MPGA.residual
+    sol = MPGA.sol
+    TF = eltype(AL)
+    (TF == Float16) ? ALF = hlu!(AL) : ALF = lu!(AL)
+    MPF = MPGEFact(AH, AL, ALF, VStore, KStore, res, sol, true, residterm)
+    return MPF
 end
 
 """
@@ -51,21 +51,21 @@ structure is immutable and MPF.AF.info cannot be changed.
 Reassigning MPG works and resuses almost all of the storage in the
 original array
 """
-function mpglu!(MPG::MPGEFact, A::AbstractArray{TW,2}) where TW <: Real
-TF=eltype(MPG.AH)
-(TF == TW) || error("Precision error in mplu!")
-AH=MPG.AH
-AH = A
-TF = eltype(MPG.AL)
-AL=MPG.AL
-AL .= TF.(A)
-(TF == Float16) ? AF = hlu!(AL) : AF = lu!(AL)
-MPG.AF.ipiv .= AF.ipiv
-VStore=MPG.VStore 
-KStore=MPG.KStore
-residterm=MPG.residterm
-MPG=MPGEFact(AH, AL, AF, VStore, KStore, MPG.residual, MPG.sol, true, residterm)
-return MPG
+function mpglu!(MPG::MPGEFact, A::AbstractArray{TW,2}) where {TW<:Real}
+    TF = eltype(MPG.AH)
+    (TF == TW) || error("Precision error in mplu!")
+    AH = MPG.AH
+    AH = A
+    TF = eltype(MPG.AL)
+    AL = MPG.AL
+    AL .= TF.(A)
+    (TF == Float16) ? AF = hlu!(AL) : AF = lu!(AL)
+    MPG.AF.ipiv .= AF.ipiv
+    VStore = MPG.VStore
+    KStore = MPG.KStore
+    residterm = MPG.residterm
+    MPG = MPGEFact(AH, AL, AF, VStore, KStore, MPG.residual, MPG.sol, true, residterm)
+    return MPG
 end
 
 
@@ -170,26 +170,30 @@ julia> xp=Float64.(A)\\b; norm(xp-mout3.sol,Inf)
 
 ```
 """
-function mpglu(A::AbstractArray{TW,2}; TF=Float32, TR=nothing, 
-          residterm=residtermdefault, basissize=10) where TW <: Real
-(TW==Float32) ? TF=Float16 : TF=TF
-MPGA=MPGArray(A; basissize=basissize, TF=TF, TR=TR)
-MPGF=mpglu!(MPGA; residterm=residterm)
-return MPGF
+function mpglu(
+    A::AbstractArray{TW,2};
+    TF = Float32,
+    TR = nothing,
+    residterm = residtermdefault,
+    basissize = 10,
+) where {TW<:Real}
+    (TW == Float32) ? TF = Float16 : TF = TF
+    MPGA = MPGArray(A; basissize = basissize, TF = TF, TR = TR)
+    MPGF = mpglu!(MPGA; residterm = residterm)
+    return MPGF
 end
 
 
 #
 # Factor a heavy MPArray and set it up for GMRES with \
 # If you want to use it with IR (why?) then set gmresok=false
 #
-function mpglu!(MPH::MPHArray; gmresok = true, basissize=10, 
-              residterm=residtermdefault)
+function mpglu!(MPH::MPHArray; gmresok = true, basissize = 10, residterm = residtermdefault)
     AH = MPH.AH
     TD = eltype(AH)
     res = MPH.residual
     sol = MPH.sol
-    n=length(res)
+    n = length(res)
     AStore = MPH.AStore
     AL = MPH.AL
     TF = eltype(AL)
@@ -203,8 +207,8 @@ function mpglu!(MPH::MPHArray; gmresok = true, basissize=10,
     AStore .= TD.(AL)
     AF = LU(AStore, ALF.ipiv, ALF.info)
     if gmresok
-        VStore=zeros(TD, n, basissize)
-        KStore=kstore(n,"gmres")
+        VStore = zeros(TD, n, basissize)
+        KStore = kstore(n, "gmres")
         MPF = MPGHFact(AH, AL, AF, VStore, KStore, res, sol, true, residterm)
     else
         MPF = MPHFact(AH, AL, AF, res, sol, true, residterm)
@@ -216,5 +220,5 @@ end
 # for CI.
 #
 function mphlu!(MPH::MPHArray)
-    MPF = mpglu!(MPH; gmresok = false, residterm=residtermdefault)
+    MPF = mpglu!(MPH; gmresok = false, residterm = residtermdefault)
 end