Use LinearAlgbera.axpy! uniformly

src/ApproxFunBase.jl

@@ -41,7 +41,8 @@ import Combinatorics: multiexponents
 import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, cross,
               qr, qr!, rank, isdiag, istril, istriu, issymmetric,
               Tridiagonal, diagm, diagm_container, factorize,
-              nullspace, Hermitian, Symmetric, adjoint, transpose, char_uplo
+              nullspace, Hermitian, Symmetric, adjoint, transpose, char_uplo,
+              axpy!
 
 import SparseArrays: blockdiag
 

src/Caching/almostbanded.jl

@@ -202,7 +202,7 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
 
     kr=co.datasize[1]+1:n
     jr=max(1,kr[1]-l):n+u
-    BLAS.axpy!(1.0,view(co.op.ops[ind],kr .- r,jr),
+    axpy!(1.0,view(co.op.ops[ind],kr .- r,jr),
                     view(co.data.bands,kr,jr))
 
     co.datasize=(n,n+u)
@@ -252,7 +252,7 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
     jr=max(ncols+1,kr[1]-l):n+u
     io∞=InterlaceOperator(io.ops[r∞,d∞])
 
-    BLAS.axpy!(1.0,view(io∞,kr.-r,jr.-ncols),view(co.data.bands,kr,jr))
+    axpy!(1.0,view(io∞,kr.-r,jr.-ncols),view(co.data.bands,kr,jr))
 
     co.datasize=(n,n+u)
     co
@@ -318,7 +318,7 @@ function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
             dind=R.u+1+k-j
             v=view(R.data,dind:dind+M-1,j)
             dt=dot(wp,v)
-            LinearAlgebra.axpy!(-2*dt,wp,v)
+            axpy!(-2*dt,wp,v)
         end
 
         # scale banded/filled entries
@@ -331,15 +331,15 @@ function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
                 @inbounds dt=muladd(conj(W[ℓ-k+1,k]),
                                     unsafe_getindex(MO.data.fill,ℓ,j),dt)
             end
-            LinearAlgebra.axpy!(-2*dt,wp2,v)
+            axpy!(-2*dt,wp2,v)
         end
 
         # scale filled entries
 
         for j=1:size(F,2)
             v=view(F,k:k+M-1,j) # the k,jth entry of F
             dt=dot(wp,v)
-            LinearAlgebra.axpy!(-2*dt,wp,v)
+            axpy!(-2*dt,wp,v)
         end
     end
     QR.ncols=col

src/Caching/banded.jl

@@ -19,7 +19,7 @@ function resizedata!(B::CachedOperator{T,<:BandedMatrix{T}},n::Integer,m_in::Int
 
         kr=B.datasize[1]+1:n
         jr=max(B.datasize[1]+1-B.data.l,1):min(n+B.data.u,M)
-        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
+        axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
 
         B.datasize = (n,m)
     end
@@ -72,7 +72,7 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T},
             dind=R.u+1+k-j
             v=view(R.data,dind:dind+M-1,j)
             dt=dot(wp,v)
-            LinearAlgebra.axpy!(-2*dt,wp,v)
+            axpy!(-2*dt,wp,v)
         end
 
         # scale banded/filled entries
@@ -81,7 +81,7 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T},
             v=view(R.data,1:M-p,j)  # shift down each time
             wp2=view(wp,p+1:M)
             dt=dot(wp2,v)
-            LinearAlgebra.axpy!(-2*dt,wp2,v)
+            axpy!(-2*dt,wp2,v)
         end
     end
     QR.ncols=col

src/Caching/bandedblockbanded.jl

@@ -31,7 +31,7 @@ function resizedata!(B::CachedOperator{T,<:BandedBlockBandedMatrix{T}}, ::Colon,
         jr=B.datasize[2]+1:col
         kr=colstart(B.data,jr[1]):colstop(B.data,jr[end])
 
-        isempty(kr) || BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
+        isempty(kr) || axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
 
         B.datasize = (last(kr),col)
     end

src/Caching/blockbanded.jl

@@ -180,7 +180,7 @@ QROperator(R::CachedOperator{T,BlockBandedMatrix{T}}) where {T} =
 #         for j=k:rowstop(R,k)
 #             v=view(R,k:cs,j)
 #             dt=dot(wp,v)
-#             LinearAlgebra.axpy!(-2*dt,wp,v)
+#             axpy!(-2*dt,wp,v)
 #         end
 #
 #         # scale banded/filled entries
@@ -189,7 +189,7 @@ QROperator(R::CachedOperator{T,BlockBandedMatrix{T}}) where {T} =
 #             v=view(R,csrt:cs,j)  # shift down each time
 #             wp2=view(wp,csrt-k+1:cs-k+1)
 #             dt=dot(wp2,v)
-#             LinearAlgebra.axpy!(-2*dt,wp2,v)
+#             axpy!(-2*dt,wp2,v)
 #         end
 #     end
 #     QR.ncols=col
@@ -309,7 +309,7 @@ function resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T},
 
                  dt = dot(view(W.data, w_j:w_j+M-1) ,
                                     view(R.data, shft + st*(ξ_2-1) +1:shft + st*(ξ_2-1) +M))
-                 BLAS.axpy!(-2*dt, view(W.data, w_j:w_j+M-1) ,
+                 axpy!(-2*dt, view(W.data, w_j:w_j+M-1) ,
                                     view(R.data, shft + st*(ξ_2-1) +1:shft + st*(ξ_2-1) +M))
              end
          end

src/Caching/matrix.jl

@@ -74,7 +74,7 @@ function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
         yp=view(Y,k:k+M-1)
 
         dt=dot(wp,yp)
-        LinearAlgebra.axpy!(-2*dt,wp,yp)
+        axpy!(-2*dt,wp,yp)
         k+=1
     end
     nz = findlast(!iszero, Y)

src/Caching/ragged.jl

@@ -41,7 +41,7 @@ function resizedata!(B::CachedOperator{T,RaggedMatrix{T}},::Colon,n::Integer) wh
 
         jr=B.datasize[2]+1:n
         kr=1:K
-        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
+        axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))
 
         B.datasize = (K,n)
 
@@ -84,7 +84,7 @@ function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T},
                 kr=J:J+length(wp)-1
                 v=view(MO.data,kr,j)
                 dt=dot(wp,v)
-                LinearAlgebra.axpy!(-2*dt,wp,v)
+                axpy!(-2*dt,wp,v)
             end
         end
     end
@@ -115,7 +115,7 @@ function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T},
         for j=k:MO.datasize[2]
             v=view(MO.data,kr,j)
             dt=dot(wp,v)
-            LinearAlgebra.axpy!(-2*dt,wp,v)
+            axpy!(-2*dt,wp,v)
         end
     end
     QR.ncols=col
@@ -221,7 +221,7 @@ for ArrTyp in (:AbstractVector, :AbstractMatrix)
             for k=n:-1:1
                 @inbounds ck = A.cols[k]
                 @inbounds u[k,c] /= A.data[ck+k-1]
-                BLAS.axpy!(-u[k,c], view(A.data,ck:ck+k-2), view(u,1:k-1,c))
+                axpy!(-u[k,c], view(A.data,ck:ck+k-2), view(u,1:k-1,c))
             end
         end
         u
@@ -270,7 +270,7 @@ function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,RaggedMatrix{T},T},T}
         yp=view(Y,k-1+(cr))
 
         dt=dot(wp,yp)
-        LinearAlgebra.axpy!(-2*dt,wp,yp)
+        axpy!(-2*dt,wp,yp)
         k+=1
     end
     nz = findlast(!iszero, Y)

src/LinearAlgebra/RaggedMatrix.jl

@@ -153,19 +153,19 @@ function mul!(y::Vector, A::RaggedMatrix, b::Vector)
     fill!(y,zero(T))
     for j=1:m
         kr=A.cols[j]:A.cols[j+1]-1
-        BLAS.axpy!(b[j],view(A.data,kr),view(y,1:length(kr)))
+        axpy!(b[j],view(A.data,kr),view(y,1:length(kr)))
     end
     y
 end
 
 
-function BLAS.axpy!(a, X::RaggedMatrix, Y::RaggedMatrix)
+function axpy!(a, X::RaggedMatrix, Y::RaggedMatrix)
     if size(X) ≠ size(Y)
         throw(BoundsError())
     end
 
     if X.cols == Y.cols
-        BLAS.axpy!(a,X.data,Y.data)
+        axpy!(a,X.data,Y.data)
     else
         for j = 1:size(X,2)
             Xn = colstop(X,j)
@@ -176,7 +176,7 @@ function BLAS.axpy!(a, X::RaggedMatrix, Y::RaggedMatrix)
                 end
             end
             cs = min(Xn,Yn)
-            BLAS.axpy!(a,view(X.data,X.cols[j]:X.cols[j]+cs-1),
+            axpy!(a,view(X.data,X.cols[j]:X.cols[j]+cs-1),
                          view(Y.data,Y.cols[j]:Y.cols[j]+cs-1))
         end
     end
@@ -188,7 +188,7 @@ colstop(X::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}},
                                             first(parentindices(X)[1]) + 1,
                             size(X,1))
 
-function BLAS.axpy!(a,X::RaggedMatrix,
+function axpy!(a,X::RaggedMatrix,
                     Y::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}}) where T
     if size(X) ≠ size(Y)
         throw(BoundsError())
@@ -213,7 +213,7 @@ function BLAS.axpy!(a,X::RaggedMatrix,
         end
 
 
-        BLAS.axpy!(a,view(X.data,kr),
+        axpy!(a,view(X.data,kr),
                     view(P.data,(P.cols[j + jsh] + ksh-1) .+ (1:length(kr))))
     end
 
@@ -235,7 +235,7 @@ function unsafe_mul!(Y::RaggedMatrix,A::RaggedMatrix,B::RaggedMatrix)
     fill!(Y.data,0)
 
     for j=1:size(B,2),k=1:colstop(B,j)
-        LinearAlgebra.axpy!(B[k,j], view(A,1:colstop(A,k),k),
+        axpy!(B[k,j], view(A,1:colstop(A,k),k),
             view(Y.data,Y.cols[j] .- 1 .+ (1:colstop(A,k))))
     end
 

src/Operators/Operator.jl

@@ -484,7 +484,7 @@ end
 
 
 # Convenience for wrapper ops
-unwrap_axpy!(α,P,A) = BLAS.axpy!(α,view(parent(P).op,P.indexes[1],P.indexes[2]),A)
+unwrap_axpy!(α,P,A) = axpy!(α,view(parent(P).op,P.indexes[1],P.indexes[2]),A)
 iswrapper(_) = false
 haswrapperstructure(_) = false
 
@@ -539,7 +539,7 @@ macro wrappergetindex(Wrap)
         Base.getindex(OP::$Wrap,k::Colon, j::ApproxFunBase.InfRanges) = view(OP, k, j)
         Base.getindex(OP::$Wrap,k::Colon, j::Colon) = view(OP, k, j)
 
-        BLAS.axpy!(α,P::ApproxFunBase.SubOperator{T,OP},A::AbstractMatrix) where {T,OP<:$Wrap} =
+        axpy!(α,P::ApproxFunBase.SubOperator{T,OP},A::AbstractMatrix) where {T,OP<:$Wrap} =
             ApproxFunBase.unwrap_axpy!(α,P,A)
 
         ApproxFunBase.mul_coefficients(A::$Wrap,b) = ApproxFunBase.mul_coefficients(A.op,b)
@@ -729,7 +729,7 @@ const WrapperOperator = Union{SpaceOperator,MultiplicationWrapper,DerivativeWrap
 ## BLAS and matrix routines
 # We assume that copy may be overriden
 
-BLAS.axpy!(a, X::Operator, Y::AbstractMatrix) = (Y .= a .* AbstractMatrix(X) .+ Y)
+axpy!(a, X::Operator, Y::AbstractMatrix) = (Y .= a .* AbstractMatrix(X) .+ Y)
 copyto!(dest::AbstractMatrix, src::Operator) = copyto!(dest, AbstractMatrix(src))
 
 # this is for operators that implement copy via axpy!
@@ -751,7 +751,7 @@ RaggedMatrix(::Type{Zeros}, V::Operator) =
 
 
 convert_axpy!(::Type{MT}, S::Operator) where {MT <: AbstractMatrix} =
-        BLAS.axpy!(one(eltype(S)), S, MT(Zeros, S))
+        axpy!(one(eltype(S)), S, MT(Zeros, S))
 
 
 

src/Operators/general/InterlaceOperator.jl

@@ -298,7 +298,7 @@ function getindex(L::InterlaceOperator{T},kr::UnitRange) where T
         if !isempty(ret_kr)
             sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]
 
-            LinearAlgebra.axpy!(1.0,L.ops[ν][sub_kr],view(ret,ret_kr))
+            axpy!(1.0,L.ops[ν][sub_kr],view(ret,ret_kr))
         end
     end
     ret
@@ -349,7 +349,7 @@ for TYP in (:BandedMatrix, :BlockBandedMatrix, :BandedBlockBandedMatrix, :Ragged
                 if !isempty(ret_kr)
                     sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]
 
-                    BLAS.axpy!(1.0,view(L.ops[ν],sub_kr,jr),view(ret,ret_kr,:))
+                    axpy!(1.0,view(L.ops[ν],sub_kr,jr),view(ret,ret_kr,:))
                 end
             end
             ret
@@ -380,7 +380,7 @@ for TYP in (:BandedMatrix, :BlockBandedMatrix, :BandedBlockBandedMatrix, :Ragged
                     sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]
                     sub_jr=cd[ret_jr[1]][2]:cd[ret_jr[end]][2]
 
-                    BLAS.axpy!(1.0,view(L.ops[ν,μ],sub_kr,sub_jr),
+                    axpy!(1.0,view(L.ops[ν,μ],sub_kr,sub_jr),
                                    view(ret,ret_kr,ret_jr))
                 end
             end

src/Operators/general/algebra.jl

@@ -116,9 +116,9 @@ for TYP in (:RaggedMatrix, :Matrix, :BandedMatrix,
     end
 end
 
-function BLAS.axpy!(α, P::SubOperator{<:Any,<:PlusOperator}, A::AbstractMatrix)
+function axpy!(α, P::SubOperator{<:Any,<:PlusOperator}, A::AbstractMatrix)
     for op in parent(P).ops
-        BLAS.axpy!(α, view(op, P.indexes[1], P.indexes[2]), A)
+        axpy!(α, view(op, P.indexes[1], P.indexes[2]), A)
     end
 
     A
@@ -216,7 +216,7 @@ end
 
 
 
-BLAS.axpy!(α, S::SubOperator{T,OP}, A::AbstractMatrix) where {T,OP<:ConstantTimesOperator} =
+axpy!(α, S::SubOperator{T,OP}, A::AbstractMatrix) where {T,OP<:ConstantTimesOperator} =
     unwrap_axpy!(α * parent(S).λ, S, A)