Merge pull request #261 from jbmuir/fiveargmul

Author: haampie

.github/workflows/ci.yml

@@ -19,7 +19,7 @@ jobs:
       matrix:
         version:
           - '1'
-          - '1.5'
+          - '1.3'
           # - 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -11,4 +11,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 RecipesBase = "0.6, 0.7, 0.8, 1.0"
-julia = "1"
+julia = "1.3"

src/bicgstabl.jl

@@ -130,10 +130,9 @@ function iterate(it::BiCGStabIterable, iteration::Int=start(it))
     F = lu!(view(it.M, L, L))
     ldiv!(it.γ, F, view(it.M, L, 1))
 
-    # This could even be BLAS 3 when combined.
-    BLAS.gemv!('N', -one(T), view(it.us, :, L), it.γ, one(T), view(it.us, :, 1))
-    BLAS.gemv!('N', one(T), view(it.rs, :, 1 : it.l), it.γ, one(T), it.x)
-    BLAS.gemv!('N', -one(T), view(it.rs, :, L), it.γ, one(T), view(it.rs, :, 1))
+    mul!(view(it.us, :, 1), view(it.us, :, L), it.γ, -one(T), one(T))
+    mul!(it.x, view(it.rs, :, 1 : it.l), it.γ, one(T), one(T))
+    mul!(view(it.rs, :, 1), view(it.rs, :, L), it.γ, -one(T), one(T))
 
     it.ω = it.γ[it.l]
     it.residual = norm(view(it.rs, :, 1))

src/gmres.jl

@@ -278,9 +278,7 @@ end
 
 function update_solution!(x, y, arnoldi::ArnoldiDecomp{T}, Pr::Identity, k::Int, Ax) where {T}
     # Update x ← x + V * y
-
-    # TODO: find the SugarBLAS alternative
-    BLAS.gemv!('N', one(T), view(arnoldi.V, :, 1 : k - 1), y, one(T), x)
+    mul!(x, view(arnoldi.V, :, 1 : k - 1), y, one(T), one(T))
 end
 
 function update_solution!(x, y, arnoldi::ArnoldiDecomp{T}, Pr, k::Int, Ax) where {T}

src/minres.jl

@@ -1,6 +1,6 @@
 export minres_iterable, minres, minres!
 using Printf
-import LinearAlgebra: BLAS.axpy!, givensAlgorithm
+import LinearAlgebra: axpy!, givensAlgorithm
 import Base: iterate
 
 mutable struct MINRESIterable{matT, solT, vecT <: DenseVector, smallVecT <: DenseVector, rotT <: Number, realT <: Real}

src/orthogonalize.jl

@@ -12,7 +12,7 @@ struct ModifiedGramSchmidt <: OrthogonalizationMethod end
 function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}, ::Type{DGKS}) where {T}
     # Orthogonalize using BLAS-2 ops
     mul!(h, adjoint(V), w)
-    BLAS.gemv!('N', -one(T), V, h, one(T), w)
+    mul!(w, V, h, -one(T), one(T))
     nrm = norm(w)
 
     # Constant used by ARPACK.
@@ -26,7 +26,7 @@ function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T},
         correction = adjoint(V)*w
         projection_size = norm(correction)
         # w = w - V * correction
-        BLAS.gemv!('N', -one(T), V, correction, one(T), w)
+        mul!(w, V, correction, -one(T), one(T))
         h .+= correction
         nrm = norm(w)
     end
@@ -40,7 +40,7 @@ end
 function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}, ::Type{ClassicalGramSchmidt}) where {T}
     # Orthogonalize using BLAS-2 ops
     mul!(h, adjoint(V), w)
-    BLAS.gemv!('N', -one(T), V, h, one(T), w)
+    mul!(w, V, h, -one(T), one(T))
     nrm = norm(w)
 
     # Normalize

src/simple.jl

@@ -35,7 +35,7 @@ function iterate(p::PowerMethodIterable, iteration::Int=start(p))
 
     # (Previous) residual vector r = Ax - λx
     copyto!(p.r, p.Ax)
-    BLAS.axpy!(-p.θ, p.x, p.r)
+    axpy!(-p.θ, p.x, p.r)
 
     # Normed residual
     p.residual = norm(p.r)