Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reorganize and cleanup to match Julia v1 conventions #1

Merged
merged 3 commits into from
Mar 24, 2023
Merged

Conversation

Omar-Elrefaei
Copy link
Member

@Omar-Elrefaei Omar-Elrefaei commented Feb 24, 2023

Implement StridedArray interface for compatibility with BLAS/LAPACK operations.

This allows all QuasiTriangular matrisies to work with LinearAlgebra.BLAS operations or our ones in ExtendedMul (in fact both do a ccall to libblas)

I think this practically makes obsolete all A_mul_B functions except the ones that specify extra arguments (ma, nc, etc...) for submatrix operations.

Reference: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-strided-arrays

for compatibility with BLAS/LAPACK operations
Copy link
Member

@MichelJuillard MichelJuillard left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is more complicated than that. We could only use BLAS/LAPACK operations if the storage of QuasiUpperTriangular matrix was quasi upper triangular. It is not the way that Triangular matrices work in Julia

using LinearAlgebra
julia> UT = UpperTriangular(rand(3,3))
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 0.249811  0.439284  0.218579
  ⋅                0.378951  0.888622
  ⋅                ⋅                 0.467821

julia> UT.data
3×3 Matrix{Float64}:
 0.249811  0.439284  0.218579
 0.689945  0.378951  0.888622
 0.540507  0.720874  0.467821

Only specialized algorithms that ignore the lower triangle can produce correct results. BLAS/LAPACK operations on UT.data would return wrong results.

In addition for QuasiUpperTriangular, we have to make sure that there is no contiguous nonzero elements on the first sub-diagonal.

@Omar-Elrefaei Omar-Elrefaei changed the title Implement StridedArray interface Reorganize and cleanup to match Julia v1 conventions Mar 23, 2023
@Omar-Elrefaei
Copy link
Member Author

Omar-Elrefaei commented Mar 23, 2023

Decisions log

  • Tests: mul!(C, Q, B) ≈ mul!(C, T, B) put our Quasi version on the left hand side to assert that we are not mutating B mistakingly (as some of the original functions where).
  • Deleted optimized mul implementations that used BLAS for a more generic version for now.
    • Bring back from git when needed. (Specialize methods on things like StridedArray, BlasFloat)
  • Remove the in-place functions A_mul_B!(A, B) that mutated A or B. There should be lmul! or rmul!. To be implemented later if needed. I actually leveraged their implementation to write the 3-arg mul!(C, A, B) in this pr.

Conventions

  • Formatting with JuliaFormatter SciMl style
  • Capitalize matrices A as opposed to vectors a
  • Use adjA to represent A'. As used in stdlib and to be able to capitalize the A. (as opposed to aAdj)

Learning notes:

  • If QuasiTriangular extends LinearAlgebra then QuasiTriangular.mul! and LinearAlgebra.mul! dispatch to the same method table no matter the type of inputs.
  • Since QuasiUpperTriangular is a subtype of AbstractMatrix, any functions not defined "safely" to the ones defined in LinearAlgebra.
    • This is true since getindex is overridden to ignore any values below the subdiagonal, and setindex is overridden to disallow assigning values below the subdiagonal.

@MichelJuillard MichelJuillard merged commit 3b37ff7 into main Mar 24, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants