axpy! error:parent must be contiguous

[andreasvarga]

After a routine update of MatrixEquations, the nightly run fails with strange errors see here, while the tests on Julia 1.6 and 1.7 run error free. All errors contain

ArgumentError: Parent must be contiguous.

occured in

[5] axpy!(alpha::Float64, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, y::SubArray{Float64, 1, Base.ReshapedArray{Float64, 1, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{UnitRange{Int64}}, false})
      @ LinearAlgebra.BLAS /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:496

[andreasvarga]

As for my part, I fixed the problem, making some vectors contiguous (avoiding reshaping of views). This probably will make the code even more efficient. With my last run, I see no more errors of type axpy! error:parent must be contiguous.

However, I assume that somebody from the Julia team is trying to also fix the problem, since I obtained new error messages (which were not present before):

ArgumentError: only support vector like inputs
  Stacktrace:
    [1] vec_pointer_stride
      @ /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:152 [inlined]
    [2] vec_pointer_stride
      @ /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:151 [inlined]
    [3] axpy!(alpha::Float64, x::SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, y::SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
      @ LinearAlgebra.BLAS /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:496

[dkarrasch]

From your fix in andreasvarga/MatrixEquations.jl@d3289c6, it looks like you were having zero-stride cases, which may have undefined behavior depending on which BLAS engine is used... or something like that. @N5N3 will probably know.

[N5N3]

The first error comes from strides for ReshapedArray. It currently only support "dense" array and StridedVector as parent. From the error message, your parent array isa SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}. IIUC you want to vec() it, so it's very likely the strides is undefined (and undefinable).

As for the second error, it hit the layout check before BLAS call. BLAS.axpy! is a level I function, thus it only support vector-like inputs. e.g.

julia> BLAS.vec_pointer_stride(view(randn(3,3),:,1:2))
(Ptr{Float64} @0x000000004778dcc0, 1)

julia> BLAS.vec_pointer_stride(view(randn(3,3),1:3,1:2))
(Ptr{Float64} @0x000000004dfe08c0, 1)

julia> BLAS.vec_pointer_stride(view(randn(3,3),1:2,1:2))
ERROR: ArgumentError: only support vector like inputs
Stacktrace:
 [1] vec_pointer_stride
   @ C:\Users\MYJ\AppData\Local\Programs\Julia-1.9.0-DEV\share\julia\stdlib\v1.9\LinearAlgebra\src\blas.jl:152 [inlined]
 [2] vec_pointer_stride(x::SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
   @ LinearAlgebra.BLAS C:\Users\MYJ\AppData\Local\Programs\Julia-1.9.0-DEV\share\julia\stdlib\v1.9\LinearAlgebra\src\blas.jl:151
 [3] top-level scope
   @ REPL[6]:1

I didn't look into your code. But I guess you'd better use LinearAlgebra.axpy! instead.

Edit: some special cases.

julia> BLAS.vec_pointer_stride(view(randn(3,3),1:2,1:1))
ERROR: ArgumentError: only support vector like inputs
Stacktrace:
 [1] vec_pointer_stride
   @ C:\Users\MYJ\AppData\Local\Programs\Julia-1.9.0-DEV\share\julia\stdlib\v1.9\LinearAlgebra\src\blas.jl:152 [inlined]
 [2] vec_pointer_stride(x::SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
   @ LinearAlgebra.BLAS C:\Users\MYJ\AppData\Local\Programs\Julia-1.9.0-DEV\share\julia\stdlib\v1.9\LinearAlgebra\src\blas.jl:151
 [3] top-level scope
   @ REPL[18]:1

julia> BLAS.vec_pointer_stride(view(randn(3,3),1:2,1))
(Ptr{Float64} @0x0000000048b16fc0, 1)

Maybe we should ingore all 1-sized dim in the check.

[andreasvarga]

The first error has been fixed.

I didn't look into your code. But I guess you'd better use LinearAlgebra.axpy! instead.

I tried it, but I get the same errors: ERROR: ArgumentError: only support vector like inputs

A line where this happens is
LinearAlgebra.axpy!(-α[1,1], v, rbar)
where rbar and v are views in Wr and Wv, respectively, defined as below:

           Wr = Matrix{T1}(undef,n,2)
            Wv = similar(Wr)
            dll = 1:dl
            ir1 = 1:n-j+1
            rbar = view(Wr,ir1,dll)
            v = view(Wv,ir1,dll)

The error occurs only for dl = 1.

[N5N3]

Well, that's the case above, which should be fixable.
As for your case, IIUC, axpy! is used only when dll == 1.
Call axpy!(-α[1,1], view(Wv,ir1,1), view(Wv,ir1,1)) should be enough. (Although I wonder why LinearAlgebra.axpy! not work.) It turns out LinearAlgebra.axpy! === BLAS.axpy!, my fault.

[andreasvarga]

Thanks. Finally it is running correctly on the nightly version. From my point of view, this issue can be closed. N5N3 ***@***.***> schrieb am Mo., 7. März 2022, 18:21:

…

Well, that's the case above, which should be fixable. As for your case, IIUC, axpy! is used only when dll == 1. Call axpy!(-α[1,1], view(Wv,ir1,1), view(Wv,ir1,1)) should be enough. (Although I wonder why LinearAlgebra.axpy! not work.) — Reply to this email directly, view it on GitHub <#916>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALJDHEE7S73JFHONLGXGJXDU6Y3LBANCNFSM5QCYOCWA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[ViralBShah]

Will get closed when the PR is merged.