We use LAPACK's QR with column pivoting.
See http://www.netlib.org/lapack/lug/node42.html. LAPACK uses blocking for
BLAS3 but that is not necessary for a generic implementation. So it's the
task is just to sort the columns by norm at each step. If you want to try
an implementation you can look for the reflector! and reflectorApply!
functions in LinAlg.
On Wednesday, September 7, 2016 at 7:53:41 AM UTC-4, Stuart Brorson wrote:
>
> Thank you, Michele & Chris for your suggestion to use qrfact. Good
> idea. Unfortunately, I neglected to mention that my matrix is a
> BigFloat array, and qrfact with pivoting fails like this:
>
> julia> Z = qrfact(A, true)
> ERROR: MethodError: `qrfact` has no method matching
> qrfact(::Array{BigFloat,2}, ::Bool)
>
> Hmmmm.... I guess I'll need to figure out how to implement the
> QR-with-pivot algorithm in Golub and Van Loan (ugh!). Or is there a
> clear reference to the algorithm Julia uses?
>
> Stuart
>
>
> On Tue, 6 Sep 2016, Chris Rackauckas wrote:
>
> > Use qrfact() (or the in-place version). It doesn't return [Q,R,P]
> because
> > Julia is smarter than that. It returns a specially designed type for the
> QR
> > return. It holds all of the those parts so you can get Q, R, and P out,
> but
> > most of the time you won't need to. For example, if A is the return of
> > qrfact, then A\b will automatically dispatch to do the fast algorithm
> for a
> > QR-factored matrix (which is the most common use for accessing Q and R).
> >
> > On Tuesday, September 6, 2016 at 6:37:59 PM UTC-7, Stuart Brorson wrote:
> >>
> >> Hello Julia users,
> >>
> >> Matlab has a variant of the QR decomposition invoked like this:
> >>
> >> [Q,R,P] = qr(A)
> >>
> >> This variant of qr() returns matrix R where the diagonal elements are
> >> sorted from largest to smallest magnitude as you go from upper left to
> >> lower right. The matrix P is the permutation matrix which permutes
> >> the rows/cols of A to give this ordering of the diagonal elements of
> >> R. That is, Q*R = A*P.
> >>
> >> I tried doing the naive, analogous thing in Julia, but get an error:
> >>
> >> julia> A = rand(3,3)
> >> 3x3 Array{Float64,2}:
> >> 0.243071 0.454947 0.89657
> >> 0.112843 0.802457 0.375417
> >> 0.154241 0.0182734 0.992542
> >>
> >> julia> Q,R,P = qr(A)
> >> ERROR: BoundsError: attempt to access (
> >> 3x3 Array{Float64,2}:
> >> -0.786117 0.0985642 -0.610168
> >> -0.364946 -0.870763 0.329523
> >> -0.498833 0.481723 0.720492,
> >>
> >> 3x3 Array{Float64,2}:
> >> -0.309204 -0.659611 -1.33693
> >> 0.0 -0.645106 0.2396
> >> 0.0 0.0 0.29177)
> >> at index [3]
> >> in indexed_next at tuple.jl:21
> >>
> >> My question: What's the best way to get the equivalent of Matlab's
> >> [Q,R,P] = qr(A) in Julia? Should I write my own qr() (not too
> >> difficult)? Or just do some row/col permutation of the output of
> >> regular qr()?
> >>
> >> Thanks for any advice,
> >>
> >> Stuart
> >>
> >> p.s. I am using Version 0.4.3-pre+6 (2015-12-11 00:38 UTC)
> >>
> >
>
