diff --git a/base/linalg.jl b/base/linalg.jl index cb8f8d36421f3..66fbd72440bad 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -224,7 +224,6 @@ include("linalg/symmetric.jl") include("linalg/diagonal.jl") include("linalg/bidiag.jl") include("linalg/uniformscaling.jl") -include("linalg/rectfullpacked.jl") include("linalg/givens.jl") include("linalg/special.jl") include("linalg/bitarray.jl") diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 97c2df36642e8..95d79e65dbe4a 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -4085,204 +4085,6 @@ for (trsen, tgsen, elty) in end end -### Rectangular full packed format - -# Symmetric rank-k operation for matrix in RFP format. -for (fn, elty, relty) in ((:dsfrk_, :Float64, :Float64), - (:ssfrk_, :Float32, :Float32), - (:zhfrk_, :Complex128, :Float64), - (:chfrk_, :Complex64, :Float32)) - @eval begin - function sfrk!(transr::Char, uplo::Char, trans::Char, alpha::Real, A::StridedMatrix{$elty}, beta::Real, C::StridedVector{$elty}) - chkuplo(uplo) - chkstride1(A) - if trans == 'N' || trans == 'n' - n, k = size(A) - elseif trans == 'T' || trans == 't' - k, n = size(A) - end - lda = max(1, stride(A, 2)) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, - Ptr{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, - Ptr{$relty}, Ptr{$elty}), - &transr, &uplo, &trans, &n, - &k, &alpha, A, &lda, - &beta, C) - C - end - end -end - -# Cholesky factorization of a real symmetric positive definite matrix A -for (fn, elty) in ((:dpftrf_, :Float64), - (:spftrf_, :Float32), - (:zpftrf_, :Complex128), - (:cpftrf_, :Complex64)) - @eval begin - function pftrf!(transr::Char, uplo::Char, A::StridedVector{$elty}) - chkuplo(uplo) - n = round(Int,div(sqrt(8length(A)), 2)) - info = Array(BlasInt, 1) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}), - &transr, &uplo, &n, A, - info) - @assertargsok - @assertnonsingular - A - end - end -end - -# Computes the inverse of a (real) symmetric positive definite matrix A using the Cholesky factorization -for (fn, elty) in ((:dpftri_, :Float64), - (:spftri_, :Float32), - (:zpftri_, :Complex128), - (:cpftri_, :Complex64)) - @eval begin - function pftri!(transr::Char, uplo::Char, A::StridedVector{$elty}) - chkuplo(uplo) - n = round(Int,div(sqrt(8length(A)), 2)) - info = Array(BlasInt, 1) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}), - &transr, &uplo, &n, A, - info) - @assertargsok - @assertnonsingular - A - end - end -end - -# DPFTRS solves a system of linear equations A*X = B with a symmetric positive definite matrix A using the Cholesky factorization -for (fn, elty) in ((:dpftrs_, :Float64), - (:spftrs_, :Float32), - (:zpftrs_, :Complex128), - (:cpftrs_, :Complex64)) - @eval begin - function pftrs!(transr::Char, uplo::Char, A::StridedVector{$elty}, B::StridedVecOrMat{$elty}) - chkuplo(uplo) - chkstride1(B) - n = round(Int,div(sqrt(8length(A)), 2)) - if n != size(B, 1) - throw(DimensionMismatch("B has first dimension $(size(B,1)) but needs $n")) - end - nhrs = size(B, 2) - ldb = max(1, stride(B, 2)) - info = Array(BlasInt, 1) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, - Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &transr, &uplo, &n, &nhrs, - A, B, &ldb, info) - @assertargsok - @assertposdef - B - end - end -end - -# Solves a matrix equation (one operand is a triangular matrix in RFP format) -for (fn, elty) in ((:dtfsm_, :Float64), - (:stfsm_, :Float32), - (:ztfsm_, :Complex128), - (:ctfsm_, :Complex64)) - @eval begin - function pftrs!(transr::Char, side::Char, uplo::Char, trans::Char, diag::Char, alpha::Real, A::StridedVector{$elty}, B::StridedMatrix{$elty}) - chkuplo(uplo) - chkside(side) - chkdiag(diag) - chkstride1(B) - m, n = size(B) - if round(Int, div(sqrt(8length(A)), 2)) != m - throw(DimensionMismatch("First dimension of B must equal $(round(Int, div(sqrt(8length(A)), 2))), got $m")) - end - ldb = max(1, stride(B, 2)) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, - Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}), - &transr, &side, &uplo, &trans, - &diag, &m, &n, &alpha, - A, B, &ldb) - B - end - end -end - -# Computes the inverse of a triangular matrix A stored in RFP format. -for (fn, elty) in ((:dtftri_, :Float64), - (:stftri_, :Float32), - (:ztftri_, :Complex128), - (:ctftri_, :Complex64)) - @eval begin - function tftri!(transr::Char, uplo::Char, diag::Char, A::StridedVector{$elty}) - chkuplo(uplo) - chkdiag(diag) - n = round(Int,div(sqrt(8length(A)), 2)) - info = Array(BlasInt, 1) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, - Ptr{$elty}, Ptr{BlasInt}), - &transr, &uplo, &diag, &n, - A, info) - @assertargsok - @assertnonsingular - A - end - end -end - -# Copies a triangular matrix from the rectangular full packed format (TF) to the standard full format (TR) -for (fn, elty) in ((:dtfttr_, :Float64), - (:stfttr_, :Float32), - (:ztfttr_, :Complex128), - (:ctfttr_, :Complex64)) - @eval begin - function tfttr!(transr::Char, uplo::Char, Arf::StridedVector{$elty}) - chkuplo(uplo) - n = round(Int,div(sqrt(8length(Arf)), 2)) - info = Array(BlasInt, 1) - A = similar(Arf, $elty, n, n) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &transr, &uplo, &n, Arf, - A, &n, info) - @assertargsok - A - end - end -end - -# Copies a triangular matrix from the standard full format (TR) to the rectangular full packed format (TF). -for (fn, elty) in ((:dtrttf_, :Float64), - (:strttf_, :Float32), - (:ztrttf_, :Complex128), - (:ctrttf_, :Complex64)) - @eval begin - function trttf!(transr::Char, uplo::Char, A::StridedMatrix{$elty}) - chkuplo(uplo) - chkstride1(A) - n = size(A, 1) - lda = max(1, stride(A, 2)) - info = Array(BlasInt, 1) - Arf = similar(A, $elty, div(n*(n+1), 2)) - ccall(($(blasfunc(fn)), liblapack), Void, - (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - &transr, &uplo, &n, A, - &lda, Arf, info) - @assertargsok - Arf - end - end -end - # Solves the real Sylvester matrix equation: op(A)*X +- X*op(B) = scale*C and A and B are both upper quasi triangular. for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), (:strsyl_, :Float32, :Float32), diff --git a/base/linalg/rectfullpacked.jl b/base/linalg/rectfullpacked.jl deleted file mode 100644 index cfb5654144d89..0000000000000 --- a/base/linalg/rectfullpacked.jl +++ /dev/null @@ -1,40 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -# Rectangular Full Packed Matrices - -type SymmetricRFP{T<:BlasFloat} <: AbstractMatrix{T} - data::Vector{T} - transr::Char - uplo::Char -end - -function Ac_mul_A_RFP{T<:BlasFloat}(A::Matrix{T}) - n = size(A, 2) - C = LAPACK.sfrk!('N', 'U', 'T', 1.0, A, 0.0, Array(T, div(n*(n+1),2))) - SymmetricRFP(C, 'N', 'U') -end - -type TriangularRFP{T<:BlasFloat} <: AbstractMatrix{T} - data::Vector{T} - transr::Char - uplo::Char -end -TriangularRFP(A::Matrix) = TriangularRFP(trttf!('N', 'U', A), 'N', 'U') - -full(A::TriangularRFP) = (A.uplo=='U' ? triu! : tril!)(LAPACK.tfttr!(A.transr, A.uplo, A.data)) - -type CholeskyDenseRFP{T<:BlasFloat} <: Factorization{T} - data::Vector{T} - transr::Char - uplo::Char -end - -cholfact!{T<:BlasFloat}(A::SymmetricRFP{T}) = CholeskyDenseRFP(LAPACK.pftrf!(A.transr, A.uplo, copy(A.data)), A.transr, A.uplo) -cholfact{T<:BlasFloat}(A::SymmetricRFP{T}) = cholfact!(copy(A)) - -copy(A::SymmetricRFP) = SymmetricRFP(copy(A.data), A.transr, A.uplo) - -# Least squares -\(A::CholeskyDenseRFP, B::VecOrMat) = LAPACK.pftrs!(A.transr, A.uplo, A.data, copy(B)) - -inv(A::CholeskyDenseRFP)=LAPACK.pftri!(A.transr, A.uplo, copy(A.data))