From 96c999892fbc372b1a76042d3e0922d5d2e3ed40 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 13 Dec 2015 07:45:56 -0500 Subject: [PATCH] Rename CHOLMOD.update! -> CHOLMOD.cholfact!/CHOLMOD.ldltfact! Move cholfact documentation to source --- base/docs/helpdb.jl | 1 - base/docs/helpdb/Base.jl | 47 ---------- base/linalg/cholesky.jl | 44 +++++++++- base/sparse/cholmod.jl | 174 +++++++++++++++++++++++++++----------- doc/stdlib/linalg.rst | 51 +++++++++-- test/sparsedir/cholmod.jl | 10 +-- 6 files changed, 214 insertions(+), 113 deletions(-) diff --git a/base/docs/helpdb.jl b/base/docs/helpdb.jl index 47d440b8cfbd3..d3b15e0474d3f 100644 --- a/base/docs/helpdb.jl +++ b/base/docs/helpdb.jl @@ -11,4 +11,3 @@ include("helpdb/Cartesian.jl") include("helpdb/Base.jl") include("helpdb/Dates.jl") include("helpdb/Pkg.jl") - diff --git a/base/docs/helpdb/Base.jl b/base/docs/helpdb/Base.jl index 50700fded25fa..b2540652cb4e8 100644 --- a/base/docs/helpdb/Base.jl +++ b/base/docs/helpdb/Base.jl @@ -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) @@ -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 diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 91623ffc4bbff..6ffe52bc049a6 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -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) @@ -122,6 +122,12 @@ 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') @@ -129,13 +135,45 @@ function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) 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) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 63ab13d7edb66..03ef1f62a86b1 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -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 @@ -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)) @@ -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 diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 8223d2a3ed0bd..1da5ecaa2fab4 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -137,27 +137,50 @@ Linear algebra functions in Julia are largely implemented by calling functions f Compute the square root of a non-negative number ``x``\ . -.. function:: cholfact(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) -> Cholesky +.. function:: cholfact(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky .. Docstring generated from Julia source - 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. + 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. -.. function:: cholfact(A; shift=0, perm=Int[]) -> CHOLMOD.Factor +.. function:: cholfact(A::StridedMatrix, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted .. Docstring generated from Julia source - 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``\ ). + 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. + +.. function:: cholfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64, + SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64}, + SuiteSparse_long}}}; shift = 0.0, perm=Int[]) -> CHOLMOD.Factor + + .. Docstring generated from Julia source + + 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``\ , 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. -.. function:: cholfact!(A [,LU=:U [,pivot=Val{false}]][;tol=-1.0]) -> Cholesky +.. function:: 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 + + .. Docstring generated from Julia source + + Compute the LDLt factorization of ``A``\ , reusing the symbolic factorization ``F``\ . + +.. function:: cholfact!(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky .. Docstring generated from Julia source - ``cholfact!`` is the same as :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)``\ . + The same as ``cholfact``\ , but saves space by overwriting the input ``A``\ , instead of creating a copy. + +.. function:: cholfact!(A::StridedMatrix, uplo::Symbol, Val{true}) -> PivotedCholesky + + .. Docstring generated from Julia source + + The same as ``cholfact``\ , but saves space by overwriting the input ``A``\ , instead of creating a copy. .. currentmodule:: Base.LinAlg @@ -193,16 +216,28 @@ Linear algebra functions in Julia are largely implemented by calling functions f Compute an ``LDLt`` factorization of a real symmetric tridiagonal matrix such that ``A = L*Diagonal(d)*L'`` where ``L`` is a unit lower triangular matrix and ``d`` is a vector. The main use of an ``LDLt`` factorization ``F = ldltfact(A)`` is to solve the linear system of equations ``Ax = b`` with ``F\b``\ . -.. function:: ldltfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; shift=0, perm=Int[]) -> CHOLMOD.Factor +.. function:: ldltfact(::Union{SparseMatrixCSC, + Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}}, + Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; + shift = 0.0, perm=Int[]) -> CHOLMOD.Factor .. Docstring generated from Julia source - 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``\ . + 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``\ . 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. +.. function:: 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 + + .. Docstring generated from Julia source + + Compute the LDLt factorization of ``A``\ , reusing the symbolic factorization ``F``\ . + .. function:: ldltfact!(::SymTridiagonal) -> LDLt .. Docstring generated from Julia source diff --git a/test/sparsedir/cholmod.jl b/test/sparsedir/cholmod.jl index 5a41ebcdaaa07..60edcb642faf6 100644 --- a/test/sparsedir/cholmod.jl +++ b/test/sparsedir/cholmod.jl @@ -409,27 +409,27 @@ for elty in (Float64, Complex{Float64}) @test F1 == F2 end - ### update! + ### cholfact!/ldltfact! F = cholfact(A1pd) CHOLMOD.change_factor!(elty, false, false, true, true, F) @test unsafe_load(F.p).is_ll == 0 CHOLMOD.change_factor!(elty, true, false, true, true, F) - @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + @test_approx_eq CHOLMOD.Sparse(cholfact!(copy(F), A1pd)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality @test size(F, 2) == 5 @test size(F, 3) == 1 @test_throws ArgumentError size(F, 0) F = cholfact(A1pdSparse, shift=2) @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty}) - @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, shift=2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + @test_approx_eq CHOLMOD.Sparse(cholfact!(copy(F), A1pd, shift=2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality F = ldltfact(A1pd) @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty}) - @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + @test_approx_eq CHOLMOD.Sparse(ldltfact!(copy(F), A1pd)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality F = ldltfact(A1pdSparse, shift=2) @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty}) - @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, shift=2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + @test_approx_eq CHOLMOD.Sparse(ldltfact!(copy(F), A1pd, shift=2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality @test isa(CHOLMOD.factor_to_sparse!(F), CHOLMOD.Sparse) @test_throws CHOLMOD.CHOLMODException CHOLMOD.factor_to_sparse!(F)