Skip to content

Special case for 2x3 matrix solve #1292

Open
@weymouth

Description

@weymouth

I'm doing some coordinate transformation stuff that requires the solution of a 2 by 3 system. Right now, all rectangular SMatrix systems are done by LinearAlgebra and the allocations mean I can't run this on the GPU. I coded this version which seems to be working and non-allocating.

@inline function (\)(a::SMatrix{2,3}, b::SVector{2})
    # columns of a'
    a1,a2 = a[1,:],a[2,:]
    
    # Q,R decomposition
    r11 = norm(a1)
    q1 = a1/r11
    r12 = q1'*a2
    p = a2-r12*q1
    r22 = norm(p) < eps(r11) ? one(r11) : norm(p)
    q2 = p/r22

    # forward substitution to solve R'v = b
    v1 = b[1]/r11
    v2 = (b[2]-r12*v1)/r22

    # return solution x = Qv 
    return q1*v1+q2*v2
end

I'm happy to add this to solve.jl and do a PR if it would be helpful.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions