From d8705234dad37c8fd73149ac142bd3af0f615ff1 Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@mathe-macht-spass.de>
Date: Thu, 29 Sep 2022 20:03:23 +0000
Subject: [PATCH] Merge branch 'fix-lsq-bug' into 'master'
Fixes segfault for least squares problems with SPQR
See merge request core/dune-istl!489
(cherry picked from commit 9a72de2d686b15fbc13a320d9e47dcfd5b322973)
d5eba453 Allow SPQR solver to work with non-square matrices
---
CHANGELOG.md | 3 +++
dune/istl/spqr.hh | 15 ++++++++++-----
2 files changed, 13 insertions(+), 5 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d2528d488..426ff793c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -35,6 +35,9 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
- 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`,
diff --git a/dune/istl/spqr.hh b/dune/istl/spqr.hh
index 4cce4468b..92a1d1626 100644
--- a/dune/istl/spqr.hh
+++ b/dune/istl/spqr.hh
@@ -286,20 +286,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_);
}
--
GitLab