@@ -6,8 +6,8 @@ import Base: (*), convert, copy, eltype, get, getindex, show, showarray, size,
6
6
linearindexing, LinearFast, LinearSlow, ctranspose
7
7
8
8
import 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
11
11
12
12
importall .. SparseArrays
13
13
@@ -1194,8 +1194,9 @@ Ac_mul_B(A::Sparse, B::VecOrMat) = Ac_mul_B(A, Dense(B))
1194
1194
1195
1195
# # Factorization methods
1196
1196
1197
+ # # Compute that symbolic factorization only
1197
1198
function 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[],
1199
1200
postorder:: Bool = true , userperm_only:: Bool = true )
1200
1201
1201
1202
sA = unsafe_load (get (A. p))
@@ -1214,85 +1215,160 @@ function fact_{Tv<:VTypes}(A::Sparse{Tv}, cm::Array{UInt8};
1214
1215
F = analyze_p (A, SuiteSparse_long[p- 1 for p in perm], cm)
1215
1216
end
1216
1217
1217
- factorize_p! (A, shift, F, cm)
1218
1218
return F
1219
1219
end
1220
1220
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 ()
1224
1223
1225
1224
# Makes it an LLt
1226
1225
unsafe_store! (common_final_ll, 1 )
1227
1226
1228
- F = fact_ (A, cm; kws... )
1227
+ # Compute the numerical factorization
1228
+ factorize_p! (A, shift, F, cm)
1229
+
1229
1230
s = unsafe_load (get (F. p))
1230
1231
s. minor < size (A, 1 ) && throw (Base. LinAlg. PosDefException (s. minor))
1231
1232
return F
1232
1233
end
1233
1234
1234
1235
"""
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)
1256
1257
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.
1257
1288
"""
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... )
1259
1293
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 ()
1263
1297
1264
1298
# Makes it an LDLt
1265
1299
unsafe_store! (common_final_ll, 0 )
1266
1300
1301
+ # Compute the numerical factorization
1302
+ factorize_p! (A, shift, F, cm)
1303
+
1267
1304
# Really make sure it's an LDLt by avoiding supernodal factorisation
1268
1305
unsafe_store! (common_supernodal, 0 )
1269
1306
1270
- F = fact_ (A, cm; kws... )
1271
1307
s = unsafe_load (get (F. p))
1272
1308
s. minor < size (A, 1 ) && throw (Base. LinAlg. ArgumentError (" matrix has one or more zero pivots" ))
1273
1309
return F
1274
1310
end
1275
1311
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
1276
1317
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 )
1284
1331
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)
1288
1337
1289
1338
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
1294
1341
end
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... )
1296
1372
1297
1373
# # Solvers
1298
1374
0 commit comments