@@ -6,8 +6,8 @@ import Base: (*), convert, copy, eltype, get, getindex, show, showarray, size,
66 linearindexing, LinearFast, LinearSlow, ctranspose
77
88import Base. LinAlg: (\ ), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
9- cholfact, det, diag, ishermitian, isposdef,
10- issym, ldltfact, logdet
9+ cholfact, cholfact!, det, diag, ishermitian, isposdef,
10+ issym, ldltfact, ldltfact!, logdet
1111
1212importall .. SparseArrays
1313
@@ -1194,8 +1194,9 @@ Ac_mul_B(A::Sparse, B::VecOrMat) = Ac_mul_B(A, Dense(B))
11941194
11951195# # Factorization methods
11961196
1197+ # # Compute that symbolic factorization only
11971198function fact_ {Tv<:VTypes} (A:: Sparse{Tv} , cm:: Array{UInt8} ;
1198- shift :: Real = 0.0 , perm:: AbstractVector{SuiteSparse_long} = SuiteSparse_long[],
1199+ perm:: AbstractVector{SuiteSparse_long} = SuiteSparse_long[],
11991200 postorder:: Bool = true , userperm_only:: Bool = true )
12001201
12011202 sA = unsafe_load (get (A. p))
@@ -1214,85 +1215,160 @@ function fact_{Tv<:VTypes}(A::Sparse{Tv}, cm::Array{UInt8};
12141215 F = analyze_p (A, SuiteSparse_long[p- 1 for p in perm], cm)
12151216 end
12161217
1217- factorize_p! (A, shift, F, cm)
12181218 return F
12191219end
12201220
1221- function cholfact (A:: Sparse ; kws... )
1222- cm = defaults (common ()) # setting the common struct to default values. Should only be done when creating new factorization.
1223- set_print_level (cm, 0 ) # no printing from CHOLMOD by default
1221+ function cholfact! {Tv} (F:: Factor{Tv} , A:: Sparse{Tv} ; shift:: Real = 0.0 )
1222+ cm = common ()
12241223
12251224 # Makes it an LLt
12261225 unsafe_store! (common_final_ll, 1 )
12271226
1228- F = fact_ (A, cm; kws... )
1227+ # Compute the numerical factorization
1228+ factorize_p! (A, shift, F, cm)
1229+
12291230 s = unsafe_load (get (F. p))
12301231 s. minor < size (A, 1 ) && throw (Base. LinAlg. PosDefException (s. minor))
12311232 return F
12321233end
12331234
12341235"""
1235- ldltfact(:: Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; shift=0, perm=Int[]) -> CHOLMOD.Factor
1236-
1237- Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix. A
1238- fill-reducing permutation is used. `F = ldltfact(A)` is most frequently used to
1239- solve systems of equations `A*x = b` with `F \\ b`. The returned factorization
1240- object `F` also supports the methods `diag `, `det`, and `logdet`. You can
1241- extract individual factors from `F` using `F[:L]`. However, since pivoting is
1242- on by default, the factorization is internally represented as `A == P'*L*D*L'*P`
1243- with a permutation matrix `P`; using just `L` without accounting for `P` will
1244- give incorrect answers. To include the effects of permutation, it's typically
1245- preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of
1246- `P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`). The complete list of
1247- supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
1248-
1249- Setting optional `shift` keyword argument computes the factorization of
1250- `A+shift*I` instead of `A`. If the `perm` argument is nonempty, it should be a
1251- permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's
1252- default AMD ordering).
1253-
1254- The function calls the C library CHOLMOD and many other functions from the
1255- library are wrapped but not exported.
1236+ cholfact!(F::Factor, A:: Union{SparseMatrixCSC,
1237+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1238+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1239+ shift = 0.0) -> CHOLMOD.Factor
1240+
1241+ Compute the LDLt factorization of `A `, reusing the symbolic factorization `F`.
1242+ """
1243+ cholfact! (F :: Factor , A :: Union {SparseMatrixCSC,
1244+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1245+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1246+ shift = 0.0 ) =
1247+ cholfact! (F, Sparse (A); shift = shift)
1248+
1249+ function cholfact (A :: Sparse ; shift :: Real = 0.0 ,
1250+ perm :: AbstractVector{SuiteSparse_long} = SuiteSparse_long[])
1251+
1252+ cm = defaults ( common ())
1253+ set_print_level (cm, 0 )
1254+
1255+ # Compute the symbolic factorization
1256+ F = fact_ (A, cm; perm = perm)
12561257
1258+ # Compute the numerical factorization
1259+ cholfact! (F, A; shift = shift)
1260+
1261+ s = unsafe_load (get (F. p))
1262+ s. minor < size (A, 1 ) && throw (Base. LinAlg. PosDefException (s. minor))
1263+ return F
1264+ end
1265+
1266+ """
1267+ cholfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,
1268+ SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},
1269+ SuiteSparse_long}}}; shift = 0.0, perm=Int[]) -> CHOLMOD.Factor
1270+
1271+ Compute the Cholesky factorization of a sparse positive definite matrix `A`.
1272+ A fill-reducing permutation is used.
1273+ `F = cholfact(A)` is most frequently used to solve systems of equations with `F\b `,
1274+ but also the methods `diag`, `det`, `logdet` are defined for `F`.
1275+ You can also extract individual factors from `F`, using `F[:L]`.
1276+ However, since pivoting is on by default,
1277+ the factorization is internally represented as `A == P'*L*L'*P` with a permutation matrix `P`;
1278+ using just `L` without accounting for `P` will give incorrect answers.
1279+ To include the effects of permutation,
1280+ it's typically preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of `P'*L`)
1281+ and `LtP = F[:UP]` (the equivalent of `L'*P`).
1282+
1283+ Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead of `A`.
1284+ If the `perm` argument is nonempty,
1285+ it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering).
1286+
1287+ The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
12571288"""
1258- ldltfact (A:: SparseMatrixCSC ; shift= 0 , perm= Int[])
1289+ cholfact (A:: Union {SparseMatrixCSC,
1290+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1291+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1292+ kws... ) = cholfact (Sparse (A); kws... )
12591293
1260- function ldltfact (A :: Sparse ; kws ... )
1261- cm = defaults ( common ()) # setting the common struct to default values. Should only be done when creating new factorization.
1262- set_print_level (cm, 0 ) # no printing from CHOLMOD by default
1294+
1295+ function ldltfact! {Tv} (F :: Factor{Tv} , A :: Sparse{Tv} ; shift :: Real = 0.0 )
1296+ cm = common ()
12631297
12641298 # Makes it an LDLt
12651299 unsafe_store! (common_final_ll, 0 )
12661300
1301+ # Compute the numerical factorization
1302+ factorize_p! (A, shift, F, cm)
1303+
12671304 # Really make sure it's an LDLt by avoiding supernodal factorisation
12681305 unsafe_store! (common_supernodal, 0 )
12691306
1270- F = fact_ (A, cm; kws... )
12711307 s = unsafe_load (get (F. p))
12721308 s. minor < size (A, 1 ) && throw (Base. LinAlg. ArgumentError (" matrix has one or more zero pivots" ))
12731309 return F
12741310end
12751311
1312+ """
1313+ ldltfact!(F::Factor, A::Union{SparseMatrixCSC,
1314+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1315+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1316+ shift = 0.0) -> CHOLMOD.Factor
12761317
1277- for f in (:cholfact , :ldltfact )
1278- @eval begin
1279- $ f (A:: SparseMatrixCSC ; kws... ) = $ f (Sparse (A); kws... )
1280- $ f (A:: Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}} ; kws... ) = $ f (Sparse (A); kws... )
1281- $ f (A:: Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}} ; kws... ) = $ f (Sparse (A); kws... )
1282- end
1283- end
1318+ Compute the LDLt factorization of `A`, reusing the symbolic factorization `F`.
1319+ """
1320+ ldltfact! (F:: Factor , A:: Union {SparseMatrixCSC,
1321+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1322+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1323+ shift = 0.0 ) =
1324+ ldltfact! (F, Sparse (A), shift = shift)
1325+
1326+ function ldltfact (A:: Sparse ; shift:: Real = 0.0 ,
1327+ perm:: AbstractVector{SuiteSparse_long} = SuiteSparse_long[])
1328+
1329+ cm = defaults (common ())
1330+ set_print_level (cm, 0 )
12841331
1285- function update! {Tv<:VTypes} (F:: Factor{Tv} , A:: Sparse{Tv} ; shift:: Real = 0.0 )
1286- cm = defaults (common ()) # setting the common struct to default values. Should only be done when creating new factorization.
1287- set_print_level (cm, 0 ) # no printing from CHOLMOD by default
1332+ # Compute the symbolic factorization
1333+ F = fact_ (A, cm; perm = perm)
1334+
1335+ # Compute the numerical factorization
1336+ ldltfact! (F, A; shift = shift)
12881337
12891338 s = unsafe_load (get (F. p))
1290- if s. is_ll!= 0
1291- unsafe_store! (common_final_ll, 1 ) # Makes it an LLt
1292- end
1293- factorize_p! (A, shift, F, cm)
1339+ s. minor < size (A, 1 ) && throw (Base. LinAlg. ArgumentError (" matrix has one or more zero pivots" ))
1340+ return F
12941341end
1295- update! {T<:VTypes} (F:: Factor{T} , A:: SparseMatrixCSC{T} ; kws... ) = update! (F, Sparse (A); kws... )
1342+
1343+ """
1344+ ldltfact(::Union{SparseMatrixCSC,
1345+ Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},
1346+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1347+ shift = 0.0, perm=Int[]) -> CHOLMOD.Factor
1348+
1349+ Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix.
1350+ A fill-reducing permutation is used.
1351+ `F = ldltfact(A)` is most frequently used to solve systems of equations `A*x = b` with `F\b `.
1352+ The returned factorization object `F` also supports the methods `diag`,
1353+ `det`, and `logdet`. You can extract individual factors from `F` using `F[:L]`.
1354+ However, since pivoting is on by default,
1355+ the factorization is internally represented as `A == P'*L*D*L'*P` with a permutation matrix `P`;
1356+ using just `L` without accounting for `P` will give incorrect answers.
1357+ To include the effects of permutation,
1358+ it's typically preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of
1359+ `P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`).
1360+ The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
1361+
1362+ Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead of `A`.
1363+ If the `perm` argument is nonempty,
1364+ it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering).
1365+
1366+ The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
1367+ """
1368+ ldltfact (A:: Union {SparseMatrixCSC,
1369+ Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
1370+ Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}};
1371+ kws... ) = ldltfact (Sparse (A); kws... )
12961372
12971373# # Solvers
12981374
0 commit comments