Skip to content

Commit

Permalink
Based on discussions in #2212.
Browse files Browse the repository at this point in the history 
Minimal invasive changes to the Factorization objects for 0.1.
All existing functionality is retained,	with only a few	renamings for clarity.

Replaces the *d routines with *fact routines:
         lud ->lufact
         chold -> cholfact
         cholpd-> cholpfact
         qrd ->qrfact
         qrpd -> qrpfact

Replaces * and Ac_mul_B on QRDense objects:
         qmulQR performs Q*A
         qTmulQR performs Q'*A

Updates to exports, documentation, tests, and deprecations.

This is a self-contained commit that can be reverted if necessary
close #2062
  • Loading branch information
@ViralBShah
ViralBShah committed Feb 10, 2013
1 parent ce4c54d commit 69e407b
Show file tree
Hide file tree
Showing 5 changed files with 100 additions and 85 deletions.
5 changes: 5 additions & 0 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,11 @@ end
@deprecate cov_pearson cov
@deprecate areduce reducedim
@deprecate tmpnam tempname
@deprecate lud lufact
@deprecate chold cholfact
@deprecate cholpd cholpfact
@deprecate qrd qrfact
@deprecate qrpd qrpfact

export grow!
function grow!(a, d)
Expand Down
22 changes: 12 additions & 10 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -517,10 +517,10 @@ export

# linear algebra
chol,
chold,
chold!,
cholpd,
cholpd!,
cholfact,
cholfact!,
cholpfact,
cholpfact!,
cond,
cross,
ctranspose,
Expand Down Expand Up @@ -549,18 +549,20 @@ export
ldltd,
linreg,
lu,
lud,
lud!,
lufact,
lufact!,
norm,
normfro,
null,
pinv,
qr,
qrd!,
qrd,
qrfact!,
qrfact,
qrp,
qrpd!,
qrpd,
qrpfact!,
qrpfact,
qmulQR,
qTmulQR,
randsym,
rank,
rref,
Expand Down
90 changes: 45 additions & 45 deletions base/linalg_dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -425,14 +425,14 @@ function inv(C::CholeskyDense)
end

## Should these functions check that the matrix is Hermitian?
chold!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
chold!{T<:BlasFloat}(A::Matrix{T}) = chold!(A, 'U')
chold{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = chold!(copy(A), UL)
chold{T<:Number}(A::Matrix{T}, UL::BlasChar) = chold(float64(A), UL)
chold{T<:Number}(A::Matrix{T}) = chold(A, 'U')
cholfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
cholfact!{T<:BlasFloat}(A::Matrix{T}) = cholfact!(A, 'U')
cholfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholfact!(copy(A), UL)
cholfact{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholfact(float64(A), UL)
cholfact{T<:Number}(A::Matrix{T}) = cholfact(A, 'U')

## Matlab (and R) compatible
chol(A::Matrix, UL::BlasChar) = factors(chold(A, UL))
chol(A::Matrix, UL::BlasChar) = factors(cholfact(A, UL))
chol(A::Matrix) = chol(A, 'U')
chol(x::Number) = imag(x) == 0 && real(x) > 0 ? sqrt(x) : error("Argument not positive-definite")

Expand Down Expand Up @@ -489,18 +489,18 @@ function inv(C::CholeskyPivotedDense)
end

## Should these functions check that the matrix is Hermitian?
cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyPivotedDense{T}(A, UL, tol)
cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(A, UL, -1.)
cholpd!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(A, 'U', tol)
cholpd!{T<:BlasFloat}(A::Matrix{T}) = cholpd!(A, 'U', -1.)
cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd(float64(A), UL, tol)
cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpd(float64(A), UL, -1.)
cholpd{T<:Number}(A::Matrix{T}, tol::Real) = cholpd(float64(A), true, tol)
cholpd{T<:Number}(A::Matrix{T}) = cholpd(float64(A), 'U', -1.)
cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd!(copy(A), UL, tol)
cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(copy(A), UL, -1.)
cholpd{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(copy(A), 'U', tol)
cholpd{T<:BlasFloat}(A::Matrix{T}) = cholpd!(copy(A), 'U', -1.)
cholpfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyPivotedDense{T}(A, UL, tol)
cholpfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpfact!(A, UL, -1.)
cholpfact!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpfact!(A, 'U', tol)
cholpfact!{T<:BlasFloat}(A::Matrix{T}) = cholpfact!(A, 'U', -1.)
cholpfact{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpfact(float64(A), UL, tol)
cholpfact{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpfact(float64(A), UL, -1.)
cholpfact{T<:Number}(A::Matrix{T}, tol::Real) = cholpfact(float64(A), true, tol)
cholpfact{T<:Number}(A::Matrix{T}) = cholpfact(float64(A), 'U', -1.)
cholpfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpfact!(copy(A), UL, tol)
cholpfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpfact!(copy(A), UL, -1.)
cholpfact{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpfact!(copy(A), 'U', tol)
cholpfact{T<:BlasFloat}(A::Matrix{T}) = cholpfact!(copy(A), 'U', -1.)

type LUDense{T} <: Factorization{T}
lu::Matrix{T}
Expand Down Expand Up @@ -530,16 +530,16 @@ function factors{T}(lu::LUDense{T})
L, U, P
end

function lud!{T<:BlasFloat}(A::Matrix{T})
function lufact!{T<:BlasFloat}(A::Matrix{T})
lu, ipiv, info = LAPACK.getrf!(A)
LUDense{T}(lu, ipiv, info)
end

lud{T<:BlasFloat}(A::Matrix{T}) = lud!(copy(A))
lud{T<:Number}(A::Matrix{T}) = lud(float64(A))
lufact{T<:BlasFloat}(A::Matrix{T}) = lufact!(copy(A))
lufact{T<:Number}(A::Matrix{T}) = lufact(float64(A))

## Matlab-compatible
lu{T<:Number}(A::Matrix{T}) = factors(lud(A))
lu{T<:Number}(A::Matrix{T}) = factors(lufact(A))
lu(x::Number) = (one(x), x, [1])

function det{T}(lu::LUDense{T})
Expand Down Expand Up @@ -571,31 +571,31 @@ end
size(A::QRDense) = size(A.hh)
size(A::QRDense,n) = size(A.hh,n)

qrd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
qrd{T<:BlasFloat}(A::StridedMatrix{T}) = qrd!(copy(A))
qrd{T<:Real}(A::StridedMatrix{T}) = qrd(float64(A))
qrfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
qrfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrfact!(copy(A))
qrfact{T<:Real}(A::StridedMatrix{T}) = qrfact(float64(A))

function factors{T<:BlasFloat}(qrd::QRDense{T})
aa = copy(qrd.hh)
function factors{T<:BlasFloat}(qrfact::QRDense{T})
aa = copy(qrfact.hh)
R = triu(aa[1:min(size(aa)),:]) # must be *before* call to orgqr!
LAPACK.orgqr!(aa, qrd.tau, size(aa,2)), R
LAPACK.orgqr!(aa, qrfact.tau, size(aa,2)), R
end

qr{T<:Number}(x::StridedMatrix{T}) = factors(qrd(x))
qr{T<:Number}(x::StridedMatrix{T}) = factors(qrfact(x))
qr(x::Number) = (one(x), x)

## Multiplication by Q from the QR decomposition
(*){T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
qmulQR{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
LAPACK.ormqr!('L', 'N', A.hh, size(A.hh,2), A.tau, copy(B))

## Multiplication by Q' from the QR decomposition
Ac_mul_B{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
qTmulQR{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
LAPACK.ormqr!('L', iscomplex(A.tau)?'C':'T', A.hh, size(A.hh,2), A.tau, copy(B))

## Least squares solution. Should be more careful about cases with m < n
function (\){T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T})
n = length(A.tau)
ans, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(A'*B)[1:n,:])
ans, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(qTmulQR(A,B))[1:n,:])
if info > 0; throw(LAPACK.SingularException(info)); end
return ans
end
Expand All @@ -615,28 +615,28 @@ end
size(x::QRPivotedDense) = size(x.hh)
size(x::QRPivotedDense,d) = size(x.hh,d)
## Multiplication by Q from the QR decomposition
(*){T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
qmulQR{T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
LAPACK.ormqr!('L', 'N', A.hh, size(A,2), A.tau, copy(B))
## Multiplication by Q' from the QR decomposition
Ac_mul_B{T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
qTmulQR{T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
LAPACK.ormqr!('L', iscomplex(A.tau)?'C':'T', A.hh, size(A,2), A.tau, copy(B))

qrpd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPivotedDense{T}(LAPACK.geqp3!(A)...)
qrpd{T<:BlasFloat}(A::StridedMatrix{T}) = qrpd!(copy(A))
qrpd{T<:Real}(x::StridedMatrix{T}) = qrpd(float64(x))
qrpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPivotedDense{T}(LAPACK.geqp3!(A)...)
qrpfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrpfact!(copy(A))
qrpfact{T<:Real}(x::StridedMatrix{T}) = qrpfact(float64(x))

function factors{T<:BlasFloat}(x::QRPivotedDense{T})
aa = copy(x.hh)
R = triu(aa[1:min(size(aa)),:])
LAPACK.orgqr!(aa, x.tau, size(aa,2)), R, x.jpvt
end

qrp{T<:BlasFloat}(x::StridedMatrix{T}) = factors(qrpd(x))
qrp{T<:BlasFloat}(x::StridedMatrix{T}) = factors(qrpfact(x))
qrp{T<:Real}(x::StridedMatrix{T}) = qrp(float64(x))

function (\){T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T})
n = length(A.tau)
x, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(A'*B)[1:n,:])
x, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(qTmulQR(A,B))[1:n,:])
if info > 0; throw(LAPACK.SingularException(info)); end
isa(B, Vector) ? x[invperm(A.jpvt)] : x[:,invperm(A.jpvt)]
end
Expand Down Expand Up @@ -674,7 +674,7 @@ function det(A::Matrix)
m, n = size(A)
if m != n; error("det only defined for square matrices"); end
if istriu(A) | istril(A); return prod(diag(A)); end
return det(lud(A))
return det(lufact(A))
end
det(x::Number) = x

Expand Down Expand Up @@ -921,7 +921,7 @@ function cond(A::StridedMatrix, p)
elseif p == 1 || p == Inf
m, n = size(A)
if m != n; error("Use 2-norm for non-square matrices"); end
cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lud(A).lu, norm(A, p))
cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lufact(A).lu, norm(A, p))
else
error("Norm type must be 1, 2 or Inf")
end
Expand Down Expand Up @@ -1217,16 +1217,16 @@ end

#show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))

lud!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
lud{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
lu(A::Tridiagonal) = factors(lud(A))
lufact!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
lufact{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
lu(A::Tridiagonal) = factors(lufact(A))

function det{T}(lu::LUTridiagonal{T})
n = length(lu.d)
prod(lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
end

det(A::Tridiagonal) = det(lud(A))
det(A::Tridiagonal) = det(lufact(A))

(\){T<:BlasFloat}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
LAPACK.gttrs!('N', lu.dl, lu.d, lu.du, lu.du2, lu.ipiv, copy(B))
Expand Down
50 changes: 29 additions & 21 deletions doc/stdlib/base.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1860,57 +1860,65 @@ Linear algebra functions in Julia are largely implemented by calling functions f

Compute the LU factorization of ``A``, such that ``A[P,:] = L*U``.

.. function:: lud(A) -> LUDense
.. function:: lufact(A) -> LUDense

Compute the LU factorization of ``A`` and return a ``LUDense`` object. ``factors(lud(A))`` returns the triangular matrices containing the factorization. The following functions are available for ``LUDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
Compute the LU factorization of ``A`` and return a ``LUDense`` object. ``factors(lufact(A))`` returns the triangular matrices containing the factorization. The following functions are available for ``LUDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.

.. function:: lud!(A) -> LUDense
.. function:: lufact!(A) -> LUDense

``lud!`` is the same as ``lud`` but saves space by overwriting the input A, instead of creating a copy.
``lufact!`` is the same as ``lufact`` but saves space by overwriting the input A, instead of creating a copy.

.. function:: chol(A, [LU]) -> F

Compute Cholesky factorization of a symmetric positive-definite matrix ``A`` and return the matrix ``F``. If ``LU`` is ``L`` (Lower), ``A = L*L'``. If ``LU`` is ``U`` (Upper), ``A = R'*R``.

.. function:: chold(A, [LU]) -> CholeskyDense
.. function:: cholfact(A, [LU]) -> CholeskyDense

Compute the Cholesky factorization of a symmetric positive-definite matrix ``A`` and return a ``CholeskyDense`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(chold(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.PosDefException`` error is thrown in case the matrix is not positive definite.
Compute the Cholesky factorization of a symmetric positive-definite matrix ``A`` and return a ``CholeskyDense`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(cholfact(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.PosDefException`` error is thrown in case the matrix is not positive definite.

.. function: chold!(A, [LU]) -> CholeskyDense
.. function: cholfact!(A, [LU]) -> CholeskyDense
``chold!`` is the same as ``chold`` but saves space by overwriting the input A, instead of creating a copy.
``cholfact!`` is the same as ``cholfact`` but saves space by overwriting the input A, instead of creating a copy.
.. function:: cholpd(A, [LU]) -> CholeskyPivotedDense
.. function:: cholpfact(A, [LU]) -> CholeskyPivotedDense

Compute the pivoted Cholesky factorization of a symmetric positive semi-definite matrix ``A`` and return a ``CholeskyDensePivoted`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(cholpd(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDensePivoted`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.RankDeficientException`` error is thrown in case the matrix is rank deficient.
Compute the pivoted Cholesky factorization of a symmetric positive semi-definite matrix ``A`` and return a ``CholeskyDensePivoted`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(cholpfact(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDensePivoted`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.RankDeficientException`` error is thrown in case the matrix is rank deficient.

.. function:: cholpd!(A, [LU]) -> CholeskyPivotedDense
.. function:: cholpfact!(A, [LU]) -> CholeskyPivotedDense

``cholpd!`` is the same as ``cholpd`` but saves space by overwriting the input A, instead of creating a copy.
``cholpfact!`` is the same as ``cholpfact`` but saves space by overwriting the input A, instead of creating a copy.

.. function:: qr(A) -> Q, R

Compute the QR factorization of ``A`` such that ``A = Q*R``. Also see ``qrd``.

.. function:: qrd(A)
.. function:: qrfact(A)

Compute the QR factorization of ``A`` and return a ``QRDense`` object. ``factors(qrd(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDense`` objects: ``size``, ``factors``, ``*``, ``Ac_mul_B``, ``\``.
Compute the QR factorization of ``A`` and return a ``QRDense`` object. ``factors(qrfact(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDense`` objects: ``size``, ``factors``, ``qmulQR``, ``qTmulQR``, ``\``.

.. function:: qrd!(A)
.. function:: qrfact!(A)

``qrd!`` is the same as ``qrd`` but saves space by overwriting the input A, instead of creating a copy.
``qrfact!`` is the same as ``qrfact`` but saves space by overwriting the input A, instead of creating a copy.

.. function:: qrp(A) -> Q, R, P

Compute the QR factorization of ``A`` with pivoting, such that ``A*I[:,P] = Q*R``, where ``I`` is the identity matrix. Also see ``qrpd``.
Compute the QR factorization of ``A`` with pivoting, such that ``A*I[:,P] = Q*R``, where ``I`` is the identity matrix. Also see ``qrpfact``.

.. function:: qrpd(A) -> QRPivotedDense
.. function:: qrpfact(A) -> QRPivotedDense

Compute the QR factorization of ``A`` with pivoting and return a ``QRDensePivoted`` object. ``factors(qrpd(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDensePivoted`` objects: ``size``, ``factors``, ``*``, ``Ac_mul_B``, ``\``.
Compute the QR factorization of ``A`` with pivoting and return a ``QRDensePivoted`` object. ``factors(qrpfact(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDensePivoted`` objects: ``size``, ``factors``, ``qmulQR``, ``qTmulQR``, ``\``.

.. function:: qrpd!(A) -> QRPivotedDense
.. function:: qrpfact!(A) -> QRPivotedDense

``qrpd!`` is the same as ``qrpd`` but saves space by overwriting the input A, instead of creating a copy.
``qrpfact!`` is the same as ``qrpfact`` but saves space by overwriting the input A, instead of creating a copy.

.. function:: qmulQR(QR, A)

Perform Q*A efficiently, where Q is a an orthogonal matrix defined as the product of k elementary reflectors from the QR decomposition.

.. function:: qTmulQR(QR, A)

Perform Q'*A efficiently, where Q is a an orthogonal matrix defined as the product of k elementary reflectors from the QR decomposition.
.. function:: eig(A) -> D, V

Expand Down
18 changes: 9 additions & 9 deletions test/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,18 +8,18 @@ for elty in (Float32, Float64, Complex64, Complex128)
apd = a'*a # symmetric positive-definite
b = convert(Vector{elty}, b)

capd = chold(apd) # upper Cholesky factor
capd = cholfact(apd) # upper Cholesky factor
r = factors(capd)
@assert_approx_eq r'*r apd
@assert_approx_eq b apd * (capd\b)
@assert_approx_eq apd * inv(capd) eye(elty, n)
@assert_approx_eq a*(capd\(a'*b)) b # least squares soln for square a
@assert_approx_eq det(capd) det(apd)

l = factors(chold(apd, 'L')) # lower Cholesky factor
l = factors(cholfact(apd, 'L')) # lower Cholesky factor
@assert_approx_eq l*l' apd

cpapd = cholpd(apd) # pivoted Choleksy decomposition
cpapd = cholpfact(apd) # pivoted Choleksy decomposition
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.LR))).<=0.) # diagonal should be non-increasing
@assert_approx_eq b apd * (cpapd\b)
Expand All @@ -32,7 +32,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
@assert_approx_eq inv(bc2) * apd eye(elty, n)
@assert_approx_eq apd * (bc2\b) b

lua = lud(a) # LU decomposition
lua = lufact(a) # LU decomposition
l,u,p = lu(a)
L,U,P = factors(lua)
@test l == L && u == U && p == P
Expand All @@ -41,18 +41,18 @@ for elty in (Float32, Float64, Complex64, Complex128)
@assert_approx_eq a * inv(lua) eye(elty, n)
@assert_approx_eq a*(lua\b) b

qra = qrd(a) # QR decomposition
qra = qrfact(a) # QR decomposition
q,r = factors(qra)
@assert_approx_eq q'*q eye(elty, n)
@assert_approx_eq q*q' eye(elty, n)
Q,R = qr(a)
@test q == Q && r == R
@assert_approx_eq q*r a
@assert_approx_eq qra*b Q*b
@assert_approx_eq qra'*b Q'*b
@assert_approx_eq qmulQR(qra,b) Q*b
@assert_approx_eq qTmulQR(qra,b) Q'*b
@assert_approx_eq a*(qra\b) b

qrpa = qrpd(a) # pivoted QR decomposition
qrpa = qrpfact(a) # pivoted QR decomposition
q,r,p = factors(qrpa)
@assert_approx_eq q'*q eye(elty, n)
@assert_approx_eq q*q' eye(elty, n)
Expand Down Expand Up @@ -294,7 +294,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
@assert_approx_eq solve(T,v) invFv
B = convert(Matrix{elty}, B)
@assert_approx_eq solve(T, B) F\B
Tlu = lud(T)
Tlu = lufact(T)
x = Tlu\v
@assert_approx_eq x invFv
@assert_approx_eq det(T) det(F)
Expand Down

0 comments on commit 69e407b

Please sign in to comment.