diff --git a/CHANGELOG.md b/CHANGELOG.md index 186604ce2f493b0087d6d7897978e6f2d256e56d..d2a89468892e21dd4e09a49dfe998b0b3a343ea1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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`, diff --git a/dune/istl/spqr.hh b/dune/istl/spqr.hh index f6403c66c1b93c67b190104f1b78205c841204f4..c2f91d2b8f6483d03a6e0509a6122b4aeb9ee989 100644 --- a/dune/istl/spqr.hh +++ b/dune/istl/spqr.hh @@ -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_); }