Comments (6)
Hi,
No, unfortunately we don't (yet) support LDLT.
from strumpack.
Hi @pghysels I'd like to better understand the steps to develop symmetric LDLT from the existing general LU. Thanks!
From the book Multifrontal Methods Parallelism, Memory Usage and Numerical Aspects (written by MUMPS's author; MUMPS has BLR LU and LDLT implementations), there should be those differences at the abstract algorithm level:
- The eliminate tree structure should be simpler: in the symmetric case, the transitive reduction of the graph of explicit dependencies is the symmetric elimination tree (Chapter 1.2.3 Elimination graph structures)
- Need 2x2 pivot to maintain the symmetry of the frontal matrices (Chapter 1.3.2 Numerical accuracy and pivoting)
- The call to BLAS kernel is different (Chapter 2.1.4 Factorization of type 1 nodes)
The above are for full-rank exact factorization. Regarding the low-rank compression of frontal matrices, I assume that the compression algorithm (RRQR) should stay exact the same (symmetric vs general)?
What I am not sure is how those textbook modifications map to Strumpack's code structure --
- I would naively assume that the low-rank compression code (
src/sparse/fronts/FrontalMatrix*.cpp
) should stay the same for symmetric and nonsymmetric matrices? - The major changes should be maded in
src/sparse/EliminationTree(MPI).hpp
and the matrix container (src/sparse/CSRMatrix(MPI).cpp
), I think? - I notice a
symm_sparse
bool attribute defined inCSRMatrix
andCompressedSparseMatrix
, but don't see how this tag is used during factorization
from strumpack.
STRUMPACK uses a symmetric pattern, so no changes would be needed in the symbolic phase.
Since we use a symmetric pattern, we need to symmetrize the pattern before passing it to say METIS, or before doing the symbolic phase. The user can specify if the pattern is already symmetric with the symm_sparse
input. But this doesn't really save much.
I think the main changes would just be to replace _getrf
with _sytrf
, and then replace the _trsm
calls (I'm not sure by what). See here and the lines below that.
Then the forward/backward solve needs to be updated accordingly. The code for distributed memory (calling ScaLAPACK p_getrf
) needs to modified, as well as the GPU code.
Keeping the low rank algorithms symmetric will be more complicated I'm afraid.
There is also the extend-add phase which passes the Schur complement (22 block of the front) to the parent in the tree. This phase now assumes we have the whole front, and not just the lower/upper triangular part.
from strumpack.
Thanks for the explanation!
Keeping the low rank algorithms symmetric will be more complicated I'm afraid.
Could you elaborate more on this point?
From the MUMPS BLR paper, the only few comments on the symmetric BLR factorization are:
- "Note that Algorithm 1 (Dense BLR LU factorization) can be easily adapted to the symmetric (definite and indefinite) case, where operations on U are not performed in order to take advantage of the symmetry."
- "The flexibility of the BLR format makes it compatible with standard threshold partial pivoting, ... In the case of symmetric matrices, only diagonal pivoting is allowed in order to maintain symmetry, so that pivots must be part of AI,I ; two-by-two pivots are used when necessary."
At least for the BLR case, the modification looks minor? For the hierachical compressions involving tree structures (HSS, HODLR), the changes will be more complicated indeed.
from strumpack.
then replace the
_trsm
calls (I'm not sure by what).
I think trsm
is the only LAPACK triangular solve routine -- the triangular factors should not be aware of the symmetry of the original matrix.
The one for LDLT should be sytrs
, which works for both diagonal or 2x2 block diagonal D
from strumpack.
The problem is that dsytrf
does not return the triangular factor. So we cannot use trsm
directly. sytrs
solves a system, so it solves with both the triangular factor and it's transpose.
from strumpack.
Related Issues (20)
- Fortran Interface for Finite Element nonlinear solver HOT 4
- Windows 10 support HOT 6
- Iterate Cluster Tree HOT 5
- mpi_type for std::pair leads to memory leaks HOT 3
- MSVC: Cannot convert from strumpack::EliminationTreeMPIDist<scalar_t,integer_t> *' to 'const strumpack::EliminationTree<scalar_t,integer_t> *' HOT 1
- Including ButterflyPACK HOT 2
- Segmentation fault during forward solve with CUDA HOT 4
- STRUMPACK with > 1 GPU per node HOT 11
- Crash with updated PETSc interface using version 7.1.3 HOT 2
- Does STRUMPACK support the 'long double' data type? HOT 1
- Symmetric support HOT 1
- Access elements in HODLRMatrix HOT 2
- Inconsistant diagonal entries for dense matrix and corresponding HODLR matrix HOT 7
- namespace "thrust" has no member "stable_sort_by_key" HOT 1
- Input matrix NAN in ComputeRange HOT 8
- HODLR OpenMP support HOT 1
- std::length_error', when total nonzeros is higher than maximum of integer4, Fortran use of STRUMPACK_MPIdist HOT 23
- HSSMatrix HSSMatrix log determinant HOT 2
- construct_partially_matrix_free for MPI HOT 2
- [error] Cmake based build with intel mkl library HOT 35
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from strumpack.