remove pointer-based mulpars implementations (#427)

Author: jishnub

src/Caching/blockbanded.jl

@@ -254,7 +254,7 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,BlockBandedMatrix{T}}}, :
          κ, ξ = blockindex.(bi)
 
          st = bs.block_strides[J1]  # the stride of the matrix
-         shft = bs.block_starts[K1,J1]-1 + st*(ξ-1) + κ-1 # the index of the pointer to the j, j entry
+         shft = bs.block_starts[K1,J1]-1 + st*(ξ-1) + κ-1 # the linear index of the j, j entry
 
 
          K_CS = Int(blockcolstop(R, Block(J1))) # last row in J-th blockcolumn
@@ -286,7 +286,7 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,BlockBandedMatrix{T}}}, :
 
          for J = J1+1:min(K1+u,J_end)
              st = bs.block_strides[J]
-             shft = bs.block_starts[K1,J] + κ-2 # the index of the pointer to the j, j entry
+             shft = bs.block_starts[K1,J] + κ-2 # the linear index of the j, j entry
              for ξ_2 = axes(bs.axes[2][Block(J)],1)
                  # we now apply I-2v*v' in place
                  RM = view(R.data, shft + st*(ξ_2-1) .+ (1:M))

src/Caching/matrix.jl

@@ -42,6 +42,9 @@ end
 
 function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
                       B::AbstractVector{T},tolerance,maxlength, inplace::Val = Val(false)) where {RR,T}
+
+    Base.require_one_based_indexing(B)
+
     A = parent(Ac)
     if length(B) > A.QR.ncols
         # upper triangularize extra columns to prepare for \
@@ -50,97 +53,40 @@ function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
 
     H=A.QR.H
     M=size(H,1)
-    m=length(B)
-    Y=_pad!!(inplace)(B,m+M+10)
-
-    k=1
-    yp=view(Y,1:M)
-    while (k ≤ m+M || norm(yp) > tolerance )
-        if k > maxlength
-            @warn "Maximum length $maxlength reached."
-            break
-        end
-
-        if k+M-1>length(Y)
-            pad!(Y,2*(k+M))
-        end
-        if k > A.QR.ncols
-            # upper triangularize extra columns to prepare for \
-            resizedata!(A.QR,:,k+M+50)
-            H=A.QR.H
-        end
-
-        wp=view(H,:,k)
-        yp=view(Y,k:k+M-1)
-
-        dt=dot(wp,yp)
-        axpy!(-2*dt,wp,yp)
-        k+=1
-    end
-    nz = findlast(!iszero, Y)
-    resize!(Y,nz === nothing ? k : min(k, nz))  # chop off zeros
-end
-
-
-# BLAS apply Q
-
-function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
-                       B::AbstractVector{T},
-                       tolerance,maxlength, inplace::Val = Val(false)) where {RR,T<:BlasFloat}
-
-    A = parent(Ac)
-
-    A_dim = size(A,1)
-
-    if length(B) > A.QR.ncols
-        # upper triangularize extra columns to prepare for \
-        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
-    end
 
-    if size(A.QR.H,1) == 1  # diagonal scaling, avoid growing
-        diagv=view(A.QR.H,1,1:length(B))
+    if M == 1  # diagonal scaling, avoid growing
         ret = inplace isa Val{true} ? B : similar(B)
-        @simd for k=1:length(ret)
-            @inbounds ret[k]=(1-2A.QR.H[1,k]^2)*B[k]
+        @simd for k in eachindex(ret)
+            @inbounds ret[k]=(1-2H[1,k]^2)*B[k]
         end
         return ret
     end
 
-    H=A.QR.H
-    h=pointer(H)
-    M=size(H,1)
-    st=stride(H,2)
-
-    sz=sizeof(T)
-
     m=length(B)
-    Y=_pad!!(inplace)(B,m+M+10)
-    y=pointer(Y)
+    Y = _pad!!(inplace)(B,m+M+10)
 
     k=1
-    yp=y
-    while (k ≤ min(m+M,A_dim) || BLAS.nrm2(M,yp,1) > tolerance )
+    yp=view(Y,1:M)
+    while (k ≤ m+M || norm(yp) > tolerance )
         if k > maxlength
             @warn "Maximum length $maxlength reached."
             break
         end
 
         if k+M-1>length(Y)
             pad!(Y,2*(k+M))
-            y=pointer(Y)
         end
         if k > A.QR.ncols
             # upper triangularize extra columns to prepare for \
             resizedata!(A.QR,:,2*(k+M))
             H=A.QR.H
-            h=pointer(H)
         end
 
-        wp=h+sz*st*(k-1)
-        yp=y+sz*(k-1)
+        wp=view(H,:,k)
+        yp=view(Y, range(k, length=M))
 
-        dt = dot(M,wp,1,yp,1)
-        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
+        dt=dot(wp,yp)
+        axpy!(-2dt,wp,yp)
         k+=1
     end
     nz = findlast(!iszero, Y)

src/Caching/ragged.jl

@@ -193,60 +193,3 @@ function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,RaggedMatrix{T},T},T}
     nz = findlast(!iszero, Y)
     resize!(Y,nz === nothing ? k : min(k, nz))  # chop off zeros
 end
-
-
-
-# BLAS apply Q
-
-function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,RaggedMatrix{T},T},T}},
-           B::AbstractVector{T},tolerance,maxlength) where {RR,T<:BlasFloat}
-    A = parent(Ac)
-    if length(B) > A.QR.ncols
-        # upper triangularize extra columns to prepare for \
-        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
-    end
-
-    H=A.QR.H
-    h=pointer(H.data)
-
-    M=size(H,1)
-    m=length(B)
-    Y=pad(B,m+M+10)
-
-    sz=sizeof(T)
-
-    k=1
-    y=pointer(Y)
-
-    yp=y
-    while (k ≤ m || BLAS.nrm2(M,yp,1) > tolerance )
-        if k > maxlength
-            @warn "Maximum length $maxlength reached."
-            break
-        end
-        if k > A.QR.ncols
-            # upper triangularize extra columns to prepare for \
-            resizedata!(A.QR,:,k+M+50)
-            H=A.QR.H
-            h=pointer(H.data)
-        end
-
-
-        M=colstop(H,k)
-
-        if k+M-1>length(Y)
-            pad!(Y,2*(k+M))
-            y=pointer(Y)
-        end
-
-
-        wp=h + sz*(H.cols[k]-1)
-        yp=y+sz*(k-1)
-
-        dt = dot(M,wp,1,yp,1)
-        BLAS.axpy!(M,-2dt,wp,1,yp,1)
-        k+=1
-    end
-    nz = findlast(!iszero, Y)
-    resize!(Y, nz === nothing ? k : min(k, nz))  # chop off zeros
-end

src/LinearAlgebra/helper.jl

@@ -237,6 +237,7 @@ function pad(v::AbstractVector,n::Integer,m::Integer)
 end
 
 function pad(A::AbstractMatrix,n::Integer,m::Integer)
+    Base.require_one_based_indexing(A)
     T=eltype(A)
 	if n <= size(A,1) && m <= size(A,2)
         A[1:n,1:m]