# How could I define my own division operation with QR decomposition?

**URL:** <https://discourse.julialang.org/t/how-could-i-define-my-own-division-operation-with-qr-decomposition/26015>\
**Category:** Numerics\
**Created:** [July 4, 2019, 4:28pm UTC](https://discourse.julialang.org/t/how-could-i-define-my-own-division-operation-with-qr-decomposition/26015 "2019-07-04T16:28:05Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Noel\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noel_araujo/32/4862_2.png) [@Noel\_Araujo](https://discourse.julialang.org/u/Noel_Araujo)\
**Post date:** [July 4, 2019, 4:28pm UTC](https://discourse.julialang.org/t/how-could-i-define-my-own-division-operation-with-qr-decomposition/26015/1 "2019-07-04T16:28:05Z")

</div>

(_Feel free to skip this part and look at the minimal example below_)

After a long time, I had to continue again with some old code, that [stillyslalom](https://discourse.julialang.org/u/stillyslalom) already posted [here](https://discourse.julialang.org/t/how-do-i-use-matlab-jl-and-vpa/18595), and because is a different question, I created a new post.  
So, in the end of the post we have the code to compute eigenvectors, given the eigenvalues:

```julia
using LinearAlgebra, GenericLinearAlgebra
import LinearAlgebra.eigen

function eigen(A::Array{T,2}) where T <: Union{BigFloat,Complex{BigFloat}}
    ni, nj = size(A)
    Λ = eigvals(A)
    if norm(imag(Λ)) < eps(BigFloat)
        Λ = real(Λ)
    end
    
    V = rand(eltype(Λ), nj, nj)
    bkp = similar(V, nj)
    for j = 1:nj
        Δ = BigFloat(Inf)
        bk = view(V,:,j)
        f = qr(A - Λ[j]*I)
        i = 0
        while Δ > eps(BigFloat) && i < 10
            bkp .= f\bk
            ck = norm(bkp)
            bkp ./= ck
            Δ = norm(bkp - bk)
            bk .= bkp
            i += 1
        end
    end
    return Λ, V
end

```

The eigenvector part seems to be a matrix multiplication operation. Then, for fun, I was debugging the code in hope to implement a GPU version of it.

My struggle is with the line `bkp .= f\bk`.

A minimal example with Float64 type:

```julia
using LinearAlgebra
A = rand(4,4)
f = qr(A)
b = rand(4)
C = f\b

```

`f` has type of `LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}`, and stores matrices `Q` and `R`.

**My questions is** : can I convert the division operation into standard matrix multiplication between `Q`,`R` and `b`?  
(I don’t mind about high tunned algorithm, the GPU will do muliplication operations faster anyway)

* * *

* * *

* * *

Extra information for reference:

- In the end. Looking at the internals of `qr.jl` file, the division operation becomes a Lapack call `LAPACK.gemqrt!('L','T',A.factors,A.T,B)`.

- A simplest code that I would like to have in GPU is:

```julia
using ArrayFire
Agpu = AFArray(A)
fgpu = qr(Agpu)
bgpu = AFArray(b)
cgpu = fgpu\bgpu

```

If I run this code, I got a message telling me that `adjoint` is not defined to `AFArray` type. So, I defined it as:  
`ArrayFire.adjoint(A::AFArray) = A'`

but I got a StackOverFlow message 😢