diff --git a/dune-solvers/solvers/tcgsolver.cc b/dune-solvers/solvers/tcgsolver.cc index 511d1016ddbaacaadb06a944dc941f209ead4de0..e67c3af70638a9eb040e22aa771f31ba5eaba28d 100644 --- a/dune-solvers/solvers/tcgsolver.cc +++ b/dune-solvers/solvers/tcgsolver.cc @@ -69,7 +69,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() *x_ = eta; // return eta - std::cout << "Abort1, radius = " << std::sqrt((*x_)*(*x_)) << std::endl; + std::cout << "Abort1, radius = " << std::sqrt(trustRegionScalarProduct(*x_,*x_)) << std::endl; if (this->historyBuffer_!="") this->writeIterate(*x_, i); @@ -92,7 +92,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() *x_ = eta_j; // return eta + tau*delta_j x_->axpy(tau, delta_j); - std::cout << "Abort2, radius = " << std::sqrt((*x_)*(*x_)) << std::endl; + std::cout << "Abort2, radius = " << std::sqrt(trustRegionScalarProduct(*x_,*x_)) << std::endl; if (this->historyBuffer_!="") this->writeIterate(*x_, i); @@ -113,7 +113,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() *x_ = eta_jPlus1; eta_j = eta_jPlus1; - std::cout << "Radius = " << sqrt((*x_)*(*x_)) << std::endl; + std::cout << "Radius = " << sqrt(trustRegionScalarProduct(*x_,*x_)) << std::endl; // ///////////////////////////////////////////// // write iteration to file, if requested diff --git a/dune-solvers/solvers/tcgsolver.hh b/dune-solvers/solvers/tcgsolver.hh index 7144e66719e656769bad55bad1ca4e36210ae752..e0fea8e26586991b3883c03500c3935a11fc4026 100644 --- a/dune-solvers/solvers/tcgsolver.hh +++ b/dune-solvers/solvers/tcgsolver.hh @@ -28,6 +28,8 @@ class TruncatedCGSolver : public IterativeSolver<VectorType> // There must be one positive and one negative root assert(root1*root2 <= 0); + std::cout << "root1: " << root1 << ", root2: " << root2 << std::endl; + return std::max(root1, root2); } @@ -71,13 +73,14 @@ public: TruncatedCGSolver(int maxIterations, double tolerance, Norm<VectorType>* errorNorm, + const MatrixType* trustRegionNormMatrix, NumProc::VerbosityMode verbosity, bool useRelativeError=true) : IterativeSolver<VectorType>(tolerance,maxIterations,verbosity,useRelativeError), matrix_(NULL), x_(NULL), rhs_(NULL), errorNorm_(errorNorm), trustRegionRadius_(0), - trustRegionNormMatrix_(NULL) + trustRegionNormMatrix_(trustRegionNormMatrix) {} void setProblem(const MatrixType& matrix,