Skip to content
Snippets Groups Projects
Commit 9a72de2d authored by Christian Engwer's avatar Christian Engwer
Browse files

Merge branch 'fix-lsq-bug' into 'master'

Fixes segfault for least squares problems with SPQR

See merge request core/dune-istl!489
parents 59188482 d5eba453
No related branches found
No related tags found
1 merge request!489Fixes segfault for least squares problems with SPQR
Pipeline #54019 passed
......@@ -30,6 +30,9 @@
- You can now use `std::tuple_element` to get the types of `MultiTypeBlockVector` entries
and `MultiTypeBlockMatrix` rows.
- The SPQR solver can now work with non-square matrices (a bug which caused a segfault when previously
attempting to do it was found and resolved).
## Deprecations and removals
- The deprecated ILU functions `bilu_backsolve`, `bilu0_decomposition`, `bilu_backsolve`,
......
......@@ -284,20 +284,25 @@ namespace Dune {
/** @brief Computes the QR decomposition. */
void decompose()
{
const std::size_t dimMat(spqrMatrix_.N());
const std::size_t nnz(spqrMatrix_.getColStart()[dimMat]);
const std::size_t nrows(spqrMatrix_.N());
const std::size_t ncols(spqrMatrix_.M());
const std::size_t nnz(spqrMatrix_.getColStart()[ncols]);
// initialise the matrix A (sorted, packed, unsymmetric, real entries)
A_ = cholmod_l_allocate_sparse(dimMat, dimMat, nnz, 1, 1, 0, 1, cc_);
A_ = cholmod_l_allocate_sparse(nrows, ncols, nnz, 1, 1, 0, 1, cc_);
// copy all the entries of Ap, Ai, Ax
for(std::size_t k = 0; k != (dimMat+1); ++k)
for(std::size_t k = 0; k != (ncols+1); ++k)
(static_cast<long int *>(A_->p))[k] = spqrMatrix_.getColStart()[k];
for(std::size_t k = 0; k != nnz; ++k)
{
(static_cast<long int*>(A_->i))[k] = spqrMatrix_.getRowIndex()[k];
(static_cast<T*>(A_->x))[k] = spqrMatrix_.getValues()[k];
}
// initialise the vector B
B_ = cholmod_l_allocate_dense(dimMat, 1, dimMat, A_->xtype, cc_);
B_ = cholmod_l_allocate_dense(nrows, 1, nrows, A_->xtype, cc_);
// compute factorization of A
spqrfactorization_=SuiteSparseQR_factorize<T>(SPQR_ORDERING_DEFAULT,SPQR_DEFAULT_TOL,A_,cc_);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment