Skip to content

Commit

Permalink
Rename CHOLMOD.update! -> CHOLMOD.cholfact!/CHOLMOD.ldltfact!
Browse files Browse the repository at this point in the history
Move cholfact documentation to source
  • Loading branch information
andreasnoack committed Dec 16, 2015
1 parent 79b2418 commit 96c9998
Show file tree
Hide file tree
Showing 6 changed files with 214 additions and 113 deletions.
1 change: 0 additions & 1 deletion base/docs/helpdb.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,3 @@ include("helpdb/Cartesian.jl")
include("helpdb/Base.jl")
include("helpdb/Dates.jl")
include("helpdb/Pkg.jl")

47 changes: 0 additions & 47 deletions base/docs/helpdb/Base.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,43 +31,6 @@ tab-delimited text to `f` by either `writedlm(f, [x y])` or by `writedlm(f, zip(
"""
writedlm

"""
cholfact(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) -> Cholesky
Compute the Cholesky factorization of a dense symmetric positive (semi)definite matrix `A`
and return either a `Cholesky` if `pivot==Val{false}` or `CholeskyPivoted` if
`pivot==Val{true}`. `LU` may be `:L` for using the lower part or `:U` for the upper part.
The default is to use `:U`. The triangular matrix can be obtained from the factorization `F`
with: `F[:L]` and `F[:U]`. The following functions are available for `Cholesky` objects:
`size`, `\\`, `inv`, `det`. For `CholeskyPivoted` there is also defined a `rank`. If
`pivot==Val{false}` a `PosDefException` exception is thrown in case the matrix is not
positive definite. The argument `tol` determines the tolerance for determining the rank. For
negative values, the tolerance is the machine precision.
"""
cholfact(A, LU=:U, pivot=Val{false})

"""
cholfact(A; shift=0, perm=Int[]) -> CHOLMOD.Factor
Compute the Cholesky factorization of a sparse positive definite matrix `A`. A fill-reducing
permutation is used. `F = cholfact(A)` is most frequently used to solve systems of equations
with `F\\b`, but also the methods `diag`, `det`, `logdet` are defined for `F`. You can also
extract individual factors from `F`, using `F[:L]`. However, since pivoting is on by
default, the factorization is internally represented as `A == P'*L*L'*P` with a permutation
matrix `P`; using just `L` without accounting for `P` will give incorrect answers. To
include the effects of permutation, it's typically preferable to extact "combined" factors
like `PtL = F[:PtL]` (the equivalent of `P'*L`) and `LtP = F[:UP]` (the equivalent of
`L'*P`).
Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead
of `A`. If the `perm` argument is nonempty, it should be a permutation of `1:size(A,1)`
giving the ordering to use (instead of CHOLMOD's default AMD ordering).
The function calls the C library CHOLMOD and many other functions from the library are
wrapped but not exported.
"""
cholfact(A)

"""
digamma(x)
Expand Down Expand Up @@ -5889,16 +5852,6 @@ indices to the parent array on the fly without checking bounds.
"""
sub

"""
cholfact!(A [,LU=:U [,pivot=Val{false}]][;tol=-1.0]) -> Cholesky
`cholfact!` is the same as [`cholfact`](:func:`cholfact`), but saves space by overwriting
the input `A`, instead of creating a copy. `cholfact!` can also reuse the symbolic
factorization from a different matrix `F` with the same structure when used as:
`cholfact!(F::CholmodFactor, A)`.
"""
cholfact!

"""
expanduser(path::AbstractString) -> AbstractString
Expand Down
44 changes: 41 additions & 3 deletions base/linalg/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,8 @@ end
"""
chol(A::AbstractMatrix) -> U
Compute the Cholesky factorization of a positive definite matrix `A` and return the
UpperTriangular matrix `U` such that `A = U'U`.
Compute the Cholesky factorization of a positive definite matrix `A`
and return the UpperTriangular matrix `U` such that `A = U'U`.
"""
function chol{T}(A::AbstractMatrix{T})
S = promote_type(typeof(chol(one(T))), Float32)
Expand Down Expand Up @@ -122,20 +122,58 @@ function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{true}}, uplo::
A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
return CholeskyPivoted{T,StridedMatrix{T}}(A, uplochar, piv, rank, tol, info)
end

"""
cholfact!(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
The same as `cholfact`, but saves space by overwriting the input `A`, instead of creating a copy.
"""
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}})
if uplo == :U
Cholesky(chol!(A, UpperTriangular).data, 'U')
else
Cholesky(chol!(A, LowerTriangular).data, 'L')
end
end
cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) = throw(ArgumentError("generic pivoted Cholesky fectorization is not implemented yet"))

"""
cholfact!(A::StridedMatrix, uplo::Symbol, Val{true}) -> PivotedCholesky
The same as `cholfact`, but saves space by overwriting the input `A`, instead of creating a copy.
"""
cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
throw(ArgumentError("generic pivoted Cholesky fectorization is not implemented yet"))
cholfact!(A::StridedMatrix, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false})

"""
cholfact(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
Compute the Cholesky factorization of a dense symmetric positive definite matrix `A`
and return a `Cholesky` factorization.
The `uplo` argument may be `:L` for using the lower part or `:U` for the upper part of `A`.
The default is to use `:U`.
The triangular Cholesky factor can be obtained from the factorization `F` with: `F[:L]` and `F[:U]`.
The following functions are available for `Cholesky` objects: `size`, `\`, `inv`, `det`.
A `PosDefException` exception is thrown in case the matrix is not positive definite.
"""
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{false})

"""
cholfact(A::StridedMatrix, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted
Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix `A`
and return a `CholeskyPivoted` factorization.
The `uplo` argument may be `:L` for using the lower part or `:U` for the upper part of `A`.
The default is to use `:U`.
The triangular Cholesky factor can be obtained from the factorization `F` with: `F[:L]` and `F[:U]`.
The following functions are available for `PivotedCholesky` objects: `size`, `\`, `inv`, `det`, and `rank`.
The argument `tol` determines the tolerance for determining the rank.
For negative values, the tolerance is the machine precision.
"""
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{true}; tol = tol)

cholfact{T}(A::StridedMatrix{T}, uplo::Symbol = :U) = cholfact(A, uplo, Val{false})

function cholfact(x::Number, uplo::Symbol=:U)
Expand Down
174 changes: 125 additions & 49 deletions base/sparse/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@ import Base: (*), convert, copy, eltype, get, getindex, show, showarray, size,
linearindexing, LinearFast, LinearSlow, ctranspose

import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
cholfact, det, diag, ishermitian, isposdef,
issym, ldltfact, logdet
cholfact, cholfact!, det, diag, ishermitian, isposdef,
issym, ldltfact, ldltfact!, logdet

importall ..SparseArrays

Expand Down Expand Up @@ -1194,8 +1194,9 @@ Ac_mul_B(A::Sparse, B::VecOrMat) = Ac_mul_B(A, Dense(B))

## Factorization methods

## Compute that symbolic factorization only
function fact_{Tv<:VTypes}(A::Sparse{Tv}, cm::Array{UInt8};
shift::Real=0.0, perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[],
perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[],
postorder::Bool=true, userperm_only::Bool=true)

sA = unsafe_load(get(A.p))
Expand All @@ -1214,85 +1215,160 @@ function fact_{Tv<:VTypes}(A::Sparse{Tv}, cm::Array{UInt8};
F = analyze_p(A, SuiteSparse_long[p-1 for p in perm], cm)
end

factorize_p!(A, shift, F, cm)
return F
end

function cholfact(A::Sparse; kws...)
cm = defaults(common()) # setting the common struct to default values. Should only be done when creating new factorization.
set_print_level(cm, 0) # no printing from CHOLMOD by default
function cholfact!{Tv}(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0)
cm = common()

# Makes it an LLt
unsafe_store!(common_final_ll, 1)

F = fact_(A, cm; kws...)
# Compute the numerical factorization
factorize_p!(A, shift, F, cm)

s = unsafe_load(get(F.p))
s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor))
return F
end

"""
ldltfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; shift=0, perm=Int[]) -> CHOLMOD.Factor
Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix. A
fill-reducing permutation is used. `F = ldltfact(A)` is most frequently used to
solve systems of equations `A*x = b` with `F\\b`. The returned factorization
object `F` also supports the methods `diag`, `det`, and `logdet`. You can
extract individual factors from `F` using `F[:L]`. However, since pivoting is
on by default, the factorization is internally represented as `A == P'*L*D*L'*P`
with a permutation matrix `P`; using just `L` without accounting for `P` will
give incorrect answers. To include the effects of permutation, it's typically
preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of
`P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`). The complete list of
supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
Setting optional `shift` keyword argument computes the factorization of
`A+shift*I` instead of `A`. If the `perm` argument is nonempty, it should be a
permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's
default AMD ordering).
The function calls the C library CHOLMOD and many other functions from the
library are wrapped but not exported.
cholfact!(F::Factor, A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
shift = 0.0) -> CHOLMOD.Factor
Compute the LDLt factorization of `A`, reusing the symbolic factorization `F`.
"""
cholfact!(F::Factor, A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
shift = 0.0) =
cholfact!(F, Sparse(A); shift = shift)

function cholfact(A::Sparse; shift::Real=0.0,
perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[])

cm = defaults(common())
set_print_level(cm, 0)

# Compute the symbolic factorization
F = fact_(A, cm; perm = perm)

# Compute the numerical factorization
cholfact!(F, A; shift = shift)

s = unsafe_load(get(F.p))
s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor))
return F
end

"""
cholfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,
SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},
SuiteSparse_long}}}; shift = 0.0, perm=Int[]) -> CHOLMOD.Factor
Compute the Cholesky factorization of a sparse positive definite matrix `A`.
A fill-reducing permutation is used.
`F = cholfact(A)` is most frequently used to solve systems of equations with `F\b`,
but also the methods `diag`, `det`, `logdet` are defined for `F`.
You can also extract individual factors from `F`, using `F[:L]`.
However, since pivoting is on by default,
the factorization is internally represented as `A == P'*L*L'*P` with a permutation matrix `P`;
using just `L` without accounting for `P` will give incorrect answers.
To include the effects of permutation,
it's typically preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of `P'*L`)
and `LtP = F[:UP]` (the equivalent of `L'*P`).
Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead of `A`.
If the `perm` argument is nonempty,
it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering).
The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
"""
ldltfact(A::SparseMatrixCSC; shift=0, perm=Int[])
cholfact(A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
kws...) = cholfact(Sparse(A); kws...)

function ldltfact(A::Sparse; kws...)
cm = defaults(common()) # setting the common struct to default values. Should only be done when creating new factorization.
set_print_level(cm, 0) # no printing from CHOLMOD by default

function ldltfact!{Tv}(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0)
cm = common()

# Makes it an LDLt
unsafe_store!(common_final_ll, 0)

# Compute the numerical factorization
factorize_p!(A, shift, F, cm)

# Really make sure it's an LDLt by avoiding supernodal factorisation
unsafe_store!(common_supernodal, 0)

F = fact_(A, cm; kws...)
s = unsafe_load(get(F.p))
s.minor < size(A, 1) && throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots"))
return F
end

"""
ldltfact!(F::Factor, A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
shift = 0.0) -> CHOLMOD.Factor
for f in (:cholfact, :ldltfact)
@eval begin
$f(A::SparseMatrixCSC; kws...) = $f(Sparse(A); kws...)
$f(A::Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}; kws...) = $f(Sparse(A); kws...)
$f(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}; kws...) = $f(Sparse(A); kws...)
end
end
Compute the LDLt factorization of `A`, reusing the symbolic factorization `F`.
"""
ldltfact!(F::Factor, A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
shift = 0.0) =
ldltfact!(F, Sparse(A), shift = shift)

function ldltfact(A::Sparse; shift::Real=0.0,
perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[])

cm = defaults(common())
set_print_level(cm, 0)

function update!{Tv<:VTypes}(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0)
cm = defaults(common()) # setting the common struct to default values. Should only be done when creating new factorization.
set_print_level(cm, 0) # no printing from CHOLMOD by default
# Compute the symbolic factorization
F = fact_(A, cm; perm = perm)

# Compute the numerical factorization
ldltfact!(F, A; shift = shift)

s = unsafe_load(get(F.p))
if s.is_ll!=0
unsafe_store!(common_final_ll, 1) # Makes it an LLt
end
factorize_p!(A, shift, F, cm)
s.minor < size(A, 1) && throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots"))
return F
end
update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}; kws...) = update!(F, Sparse(A); kws...)

"""
ldltfact(::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
shift = 0.0, perm=Int[]) -> CHOLMOD.Factor
Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix.
A fill-reducing permutation is used.
`F = ldltfact(A)` is most frequently used to solve systems of equations `A*x = b` with `F\b`.
The returned factorization object `F` also supports the methods `diag`,
`det`, and `logdet`. You can extract individual factors from `F` using `F[:L]`.
However, since pivoting is on by default,
the factorization is internally represented as `A == P'*L*D*L'*P` with a permutation matrix `P`;
using just `L` without accounting for `P` will give incorrect answers.
To include the effects of permutation,
it's typically preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of
`P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`).
The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead of `A`.
If the `perm` argument is nonempty,
it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering).
The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
"""
ldltfact(A::Union{SparseMatrixCSC,
Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
kws...) = ldltfact(Sparse(A); kws...)

## Solvers

Expand Down
Loading

0 comments on commit 96c9998

Please sign in to comment.