diff --git a/dune/solvers/iterationsteps/cgstep.cc b/dune/solvers/iterationsteps/cgstep.cc index 13fa1432b21f1703bba20aacabd873ab78d3ecb6..7b62b7ad50d2efa3d8fdf77cb38d04c9b3b299b0 100644 --- a/dune/solvers/iterationsteps/cgstep.cc +++ b/dune/solvers/iterationsteps/cgstep.cc @@ -9,10 +9,10 @@ template <class MatrixType, class VectorType> void CGStep<MatrixType, VectorType>::preprocess() { // Compute the residual (r starts out as the rhs) - matrix_.mmv(x_,r_); + matrix_->mmv(x_,r_); if (preconditioner_) { - preconditioner_->setMatrix(matrix_); + preconditioner_->setMatrix(*matrix_); preconditioner_->apply(p_, r_); } else p_ = r_; @@ -25,7 +25,7 @@ void CGStep<MatrixType, VectorType>::iterate() { VectorType q(x_); - matrix_.mv(p_, q); // q_0 = Ap_0 + matrix_->mv(p_, q); // q_0 = Ap_0 const double alpha = r_squared_old_ / (p_ * q); // alpha_0 = r_0*r_0/p_0*Ap_0 x_.axpy(alpha, p_); // x_1 = x_0 + alpha_0 p_0 r_.axpy(-alpha, q); // r_1 = r_0 - alpha_0 Ap_0 diff --git a/dune/solvers/iterationsteps/cgstep.hh b/dune/solvers/iterationsteps/cgstep.hh index 12e9dbb74c0e624969ed430d7e0821091f74b329..7f65fb0975db360c322ab6f2511aef561239d55e 100644 --- a/dune/solvers/iterationsteps/cgstep.hh +++ b/dune/solvers/iterationsteps/cgstep.hh @@ -18,7 +18,7 @@ namespace Dune { CGStep(const MatrixType& matrix, VectorType& x, const VectorType& rhs) - : p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix), + : p_(rhs.size()), r_(rhs), x_(x), matrix_(Dune::stackobject_to_shared_ptr(matrix)), preconditioner_(nullptr) {} @@ -33,7 +33,7 @@ namespace Dune { //! Set linear operator, solution and right hand side virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) { - matrix_ = mat; + matrix_ = Dune::stackobject_to_shared_ptr(mat); x_ = x; r_ = rhs; } @@ -48,7 +48,7 @@ namespace Dune { VectorType p_; // search direction VectorType r_; // residual VectorType& x_; - const MatrixType& matrix_; + std::shared_ptr<const MatrixType> matrix_; double r_squared_old_; Preconditioner<MatrixType, VectorType>* preconditioner_; };