[Dune] [Dune-Commit] dune-istl r832 - trunk/istl
Oliver Sander
sander at mi.fu-berlin.de
Tue Nov 27 12:05:02 CET 2007
>
> Don't get disturbed by the other changes, that are just spacings.
>
Good to know that the other ones are just spacings. But how
am I to find the relevant 2% of this commit message?
If you want to change the indentation please use a separate
commit.
Thanks,
Oliver
>
>
> Modified: trunk/istl/solvers.hh
> ===================================================================
> --- trunk/istl/solvers.hh 2007-11-05 10:01:49 UTC (rev 831)
> +++ trunk/istl/solvers.hh 2007-11-26 21:01:35 UTC (rev 832)
> @@ -19,7 +19,7 @@
> @ingroup ISTL
> */
> /** @addtogroup ISTL_Solvers
> - @{
> + @{
> */
>
>
> @@ -42,35 +42,35 @@
> struct InverseOperatorResult
> {
> /** \brief Default constructor */
> - InverseOperatorResult ()
> - {
> - clear();
> - }
> + InverseOperatorResult ()
> + {
> + clear();
> + }
>
> /** \brief Resets all data */
> - void clear ()
> - {
> - iterations = 0;
> - reduction = 0;
> - converged = false;
> - conv_rate = 1;
> - elapsed = 0;
> - }
> + void clear ()
> + {
> + iterations = 0;
> + reduction = 0;
> + converged = false;
> + conv_rate = 1;
> + elapsed = 0;
> + }
>
> /** \brief Number of iterations */
> - int iterations;
> + int iterations;
>
> /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */
> - double reduction;
> + double reduction;
>
> /** \brief True if convergence criterion has been met */
> - bool converged;
> + bool converged;
>
> /** \brief Convergence rate (average reduction per step) */
> - double conv_rate;
> + double conv_rate;
>
> /** \brief Elapsed time in seconds */
> - double elapsed;
> + double elapsed;
> };
>
>
> @@ -99,13 +99,13 @@
> typedef typename X::field_type field_type;
>
> /**
> - \brief Apply inverse operator,
> -
> - \warning Note: right hand side b may be overwritten!
> + \brief Apply inverse operator,
> +
> + \warning Note: right hand side b may be overwritten!
>
> - \param x The left hand side to store the result in.
> - \param b The right hand side
> - \param res Object to store the statistics about applying the operator.
> + \param x The left hand side to store the result in.
> + \param b The right hand side
> + \param res Object to store the statistics about applying the operator.
> */
> virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
>
> @@ -123,6 +123,41 @@
>
> //! \brief Destructor
> virtual ~InverseOperator () {}
> +
> + protected:
> + // spacing values
> + enum { iterationSpacing = 5 , normSpacing = 16 };
> +
> + //! helper function for printing header of solver output
> + void printHeader(std::ostream& s) const
> + {
> + s << std::setw(iterationSpacing) << " Iter";
> + s << std::setw(normSpacing) << "Defect";
> + s << std::setw(normSpacing) << "Rate" << std::endl;
> + }
> +
> + //! helper function for printing solver output
> + template <class DataType>
> + void printOutput(std::ostream& s,
> + const int iter,
> + const DataType& norm,
> + const DataType& norm_old) const
> + {
> + const DataType rate = norm/norm_old;
> + s << std::setw(iterationSpacing) << iter << " ";
> + s << std::setw(normSpacing) << norm << " ";
> + s << std::setw(normSpacing) << rate << std::endl;
> + }
> +
> + //! helper function for printing solver output
> + template <class DataType>
> + void printOutput(std::ostream& s,
> + const int iter,
> + const DataType& norm) const
> + {
> + s << std::setw(iterationSpacing) << iter << " ";
> + s << std::setw(normSpacing) << norm << std::endl;
> + }
> };
>
>
> @@ -169,7 +204,7 @@
> */
> template<class L, class P>
> LoopSolver (L& op, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -177,28 +212,28 @@
> }
>
> /**
> - \brief Set up loop solver
> -
> - \param op The operator we solve.
> - \param sp The scalar product to use, e. g. SeqScalarproduct.
> - \param prec The preconditioner to apply in each iteration of the loop.
> - Has to inherit from Preconditioner.
> - \param reduction The relative defect reduction to achieve when applying
> - the operator.
> - \param maxit The maximum number of iteration steps allowed when applying
> - the operator.
> - \param verbose The verbosity level.
> -
> - Verbose levels are:
> - <ul>
> - <li> 0 : print nothing </li>
> - <li> 1 : print initial and final defect and statistics </li>
> - <li> 2 : print line for each iteration </li>
> - </ul>
> + \brief Set up loop solver
> +
> + \param op The operator we solve.
> + \param sp The scalar product to use, e. g. SeqScalarproduct.
> + \param prec The preconditioner to apply in each iteration of the loop.
> + Has to inherit from Preconditioner.
> + \param reduction The relative defect reduction to achieve when applying
> + the operator.
> + \param maxit The maximum number of iteration steps allowed when applying
> + the operator.
> + \param verbose The verbosity level.
> +
> + Verbose levels are:
> + <ul>
> + <li> 0 : print nothing </li>
> + <li> 1 : print initial and final defect and statistics </li>
> + <li> 2 : print line for each iteration </li>
> + </ul>
> */
> template<class L, class S, class P>
> LoopSolver (L& op, S& sp, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -206,92 +241,96 @@
> }
>
>
> - //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
> - virtual void apply (X& x, X& b, InverseOperatorResult& res)
> - {
> - // clear solver statistics
> - res.clear();
> + //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
> + virtual void apply (X& x, X& b, InverseOperatorResult& res)
> + {
> + // clear solver statistics
> + res.clear();
>
> - // start a timer
> - Timer watch;
> + // start a timer
> + Timer watch;
>
> - // prepare preconditioner
> - _prec.pre(x,b);
> + // prepare preconditioner
> + _prec.pre(x,b);
>
> - // overwrite b with defect
> - _op.applyscaleadd(-1,x,b);
> + // overwrite b with defect
> + _op.applyscaleadd(-1,x,b);
>
> - // compute norm, \todo parallelization
> - double def0 = _sp.norm(b);
> + // compute norm, \todo parallelization
> + double def0 = _sp.norm(b);
>
> - // printing
> - if (_verbose>0)
> - {
> - std::cout << "=== LoopSolver" << std::endl;
> - if (_verbose>1)
> - std::cout << " Iter Defect Rate" << std::endl;
> - if (_verbose>1)
> - std::cout << "0 " << def0 << std::endl;
> - }
> + // printing
> + if (_verbose>0)
> + {
> + std::cout << "=== LoopSolver" << std::endl;
> + if (_verbose>1)
> + {
> + this->printHeader(std::cout);
> + this->printOutput(std::cout,0,def0);
> + }
> + }
>
> - // allocate correction vector
> - X v(x);
> + // allocate correction vector
> + X v(x);
>
> - // iteration loop
> - int i=1; double def=def0;
> - for ( ; i<=_maxit; i++ )
> - {
> - v = 0; // clear correction
> - _prec.apply(v,b); // apply preconditioner
> - x += v; // update solution
> - _op.applyscaleadd(-1,v,b); // update defect
> - double defnew=_sp.norm(b);// comp defect norm
> - if (_verbose>1) // print
> - std::cout << i << " " << defnew << " " << defnew/def << std::endl;
> - def = defnew; // update norm
> - if (def<def0*_reduction || def<1E-30) // convergence check
> - {
> - res.converged = true;
> - break;
> - }
> - }
> + // iteration loop
> + int i=1; double def=def0;
> + for ( ; i<=_maxit; i++ )
> + {
> + v = 0; // clear correction
> + _prec.apply(v,b); // apply preconditioner
> + x += v; // update solution
> + _op.applyscaleadd(-1,v,b); // update defect
> + double defnew=_sp.norm(b);// comp defect norm
> + if (_verbose>1) // print
> + this->printOutput(std::cout,i,defnew,def);
> + //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
> + def = defnew; // update norm
> + if (def<def0*_reduction || def<1E-30) // convergence check
> + {
> + res.converged = true;
> + break;
> + }
> + }
>
> - // print
> - if (_verbose==1)
> - std::cout << i << " " << def << std::endl;
> -
> - // postprocess preconditioner
> - _prec.post(x);
> + // print
> + if (_verbose==1)
> + this->printOutput(std::cout,i,def);
> +
> + // postprocess preconditioner
> + _prec.post(x);
>
> - // fill statistics
> - res.iterations = i;
> - res.reduction = def/def0;
> - res.conv_rate = pow(res.reduction,1.0/i);
> - res.elapsed = watch.elapsed();
> + // fill statistics
> + res.iterations = i;
> + res.reduction = def/def0;
> + res.conv_rate = pow(res.reduction,1.0/i);
> + res.elapsed = watch.elapsed();
>
> - // final print
> - if (_verbose>0)
> - std::cout << "=== rate=" << res.conv_rate
> - << ", T=" << res.elapsed
> - << ", TIT=" << res.elapsed/i
> - << ", IT=" << i << std::endl;
> - }
> + // final print
> + if (_verbose>0)
> + {
> + std::cout << "=== rate=" << res.conv_rate
> + << ", T=" << res.elapsed
> + << ", TIT=" << res.elapsed/i
> + << ", IT=" << i << std::endl;
> + }
> + }
>
> - //! \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
> - virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
> - {
> - _reduction = reduction;
> - (*this).apply(x,b,res);
> - }
> + //! \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
> + virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
> + {
> + _reduction = reduction;
> + (*this).apply(x,b,res);
> + }
>
> private:
> - SeqScalarProduct<X> ssp;
> - LinearOperator<X,X>& _op;
> - Preconditioner<X,X>& _prec;
> - ScalarProduct<X>& _sp;
> - double _reduction;
> - int _maxit;
> - int _verbose;
> + SeqScalarProduct<X> ssp;
> + LinearOperator<X,X>& _op;
> + Preconditioner<X,X>& _prec;
> + ScalarProduct<X>& _sp;
> + double _reduction;
> + int _maxit;
> + int _verbose;
> };
>
>
> @@ -314,7 +353,7 @@
> */
> template<class L, class P>
> GradientSolver (L& op, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -327,7 +366,7 @@
> */
> template<class L, class S, class P>
> GradientSolver (L& op, S& sp, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -348,42 +387,44 @@
>
> X p(x); // create local vectors
> X q(b);
> -
> +
> double def0 = _sp.norm(b);// compute norm
>
> if (_verbose>0) // printing
> - {
> - std::cout << "=== GradientSolver" << std::endl;
> - if (_verbose>1) {
> - std::cout << " Iter Defect Rate" << std::endl;
> - std::cout << "0 " << def0 << std::endl;
> - }
> - }
> + {
> + std::cout << "=== GradientSolver" << std::endl;
> + if (_verbose>1)
> + {
> + this->printHeader(std::cout);
> + this->printOutput(std::cout,0,def0);
> + }
> + }
>
> int i=1; double def=def0; // loop variables
> field_type lambda;
> - for ( ; i<=_maxit; i++ )
> - {
> - p = 0; // clear correction
> - _prec.apply(p,b); // apply preconditioner
> - _op.apply(p,q); // q=Ap
> - lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization
> - x.axpy(lambda,p); // update solution
> - b.axpy(-lambda,q); // update defect
> -
> - double defnew=_sp.norm(b);// comp defect norm
> - if (_verbose>1) // print
> - std::cout << i << " " << defnew << " " << defnew/def << std::endl;
> - def = defnew; // update norm
> - if (def<def0*_reduction || def<1E-30) // convergence check
> - {
> - res.converged = true;
> - break;
> - }
> - }
> -
> + for ( ; i<=_maxit; i++ )
> + {
> + p = 0; // clear correction
> + _prec.apply(p,b); // apply preconditioner
> + _op.apply(p,q); // q=Ap
> + lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization
> + x.axpy(lambda,p); // update solution
> + b.axpy(-lambda,q); // update defect
> +
> + double defnew=_sp.norm(b);// comp defect norm
> + if (_verbose>1) // print
> + this->printOutput(std::cout,i,defnew,def);
> +
> + def = defnew; // update norm
> + if (def<def0*_reduction || def<1E-30) // convergence check
> + {
> + res.converged = true;
> + break;
> + }
> + }
> +
> if (_verbose==1) // printing for non verbose
> - std::cout << i << " " << def << std::endl;
> + this->printOutput(std::cout,i,def);
>
> _prec.post(x); // postprocess preconditioner
> res.iterations = i; // fill statistics
> @@ -409,13 +450,13 @@
> }
>
> private:
> - SeqScalarProduct<X> ssp;
> - LinearOperator<X,X>& _op;
> - Preconditioner<X,X>& _prec;
> - ScalarProduct<X>& _sp;
> - double _reduction;
> - int _maxit;
> - int _verbose;
> + SeqScalarProduct<X> ssp;
> + LinearOperator<X,X>& _op;
> + Preconditioner<X,X>& _prec;
> + ScalarProduct<X>& _sp;
> + double _reduction;
> + int _maxit;
> + int _verbose;
> };
>
>
> @@ -463,112 +504,118 @@
> */
> virtual void apply (X& x, X& b, InverseOperatorResult& res)
> {
> - res.clear(); // clear solver statistics
> - Timer watch; // start a timer
> - _prec.pre(x,b); // prepare preconditioner
> - _op.applyscaleadd(-1,x,b); // overwrite b with defect
> -
> - X p(x); // the search direction
> - X q(x); // a temporary vector
> -
> - double def0 = _sp.norm(b);// compute norm
> - if (def0<1E-30) // convergence check
> - {
> - res.converged = true;
> - res.iterations = 0; // fill statistics
> - res.reduction = 0;
> - res.conv_rate = 0;
> - res.elapsed=0;
> - if (_verbose>0) // final print
> - std::cout << "=== rate=" << res.conv_rate
> - << ", T=" << res.elapsed << ", TIT=" << res.elapsed
> - << ", IT=0" << std::endl;
> - return;
> - }
> + res.clear(); // clear solver statistics
> + Timer watch; // start a timer
> + _prec.pre(x,b); // prepare preconditioner
> + _op.applyscaleadd(-1,x,b); // overwrite b with defect
> +
> + X p(x); // the search direction
> + X q(x); // a temporary vector
> +
> + double def0 = _sp.norm(b);// compute norm
> + if (def0<1E-30) // convergence check
> + {
> + res.converged = true;
> + res.iterations = 0; // fill statistics
> + res.reduction = 0;
> + res.conv_rate = 0;
> + res.elapsed=0;
> + if (_verbose>0) // final print
> + std::cout << "=== rate=" << res.conv_rate
> + << ", T=" << res.elapsed << ", TIT=" << res.elapsed
> + << ", IT=0" << std::endl;
> + return;
> + }
>
> - if (_verbose>0) // printing
> - {
> - std::cout << "=== CGSolver" << std::endl;
> - if (_verbose>1) {
> - std::cout << " Iter Defect Rate" << std::endl;
> - std::cout << "0 " << def0 << std::endl;
> - }
> - }
> + if (_verbose>0) // printing
> + {
> + std::cout << "=== CGSolver" << std::endl;
> + if (_verbose>1) {
> + this->printHeader(std::cout);
> + this->printOutput(std::cout,0,def0);
> + }
> + }
>
> - // some local variables
> - double def=def0; // loop variables
> - field_type rho,rholast,lambda,alpha,beta;
> + // some local variables
> + double def=def0; // loop variables
> + field_type rho,rholast,lambda,alpha,beta;
>
> - // determine initial search direction
> - p = 0; // clear correction
> - _prec.apply(p,b); // apply preconditioner
> - rholast = _sp.dot(p,b); // orthogonalization
> + // determine initial search direction
> + p = 0; // clear correction
> + _prec.apply(p,b); // apply preconditioner
> + rholast = _sp.dot(p,b); // orthogonalization
>
> - // the loop
> - int i=1;
> - for ( ; i<=_maxit; i++ )
> - {
> - // minimize in given search direction p
> - _op.apply(p,q); // q=Ap
> - alpha = _sp.dot(p,q); // scalar product
> - lambda = rholast/alpha; // minimization
> - x.axpy(lambda,p); // update solution
> - b.axpy(-lambda,q); // update defect
> + // the loop
> + int i=1;
> + for ( ; i<=_maxit; i++ )
> + {
> + // minimize in given search direction p
> + _op.apply(p,q); // q=Ap
> + alpha = _sp.dot(p,q); // scalar product
> + lambda = rholast/alpha; // minimization
> + x.axpy(lambda,p); // update solution
> + b.axpy(-lambda,q); // update defect
>
> - // convergence test
> - double defnew=_sp.norm(b);// comp defect norm
> - if (_verbose>1) // print
> - std::cout << i << " " << defnew << " " << defnew/def << std::endl;
> - def = defnew; // update norm
> - if (def<def0*_reduction || def<1E-30 || i==_maxit) // convergence check
> - {
> - res.converged = true;
> - break;
> - }
> + // convergence test
> + double defnew=_sp.norm(b);// comp defect norm
>
> - // determine new search direction
> - q = 0; // clear correction
> - _prec.apply(q,b); // apply preconditioner
> - rho = _sp.dot(q,b); // orthogonalization
> - beta = rho/rholast; // scaling factor
> - p *= beta; // scale old search direction
> - p += q; // orthogonalization with correction
> - rholast = rho; // remember rho for recurrence
> - }
> + if (_verbose>1) // print
> + this->printOutput(std::cout,i,defnew,def);
>
> - if (_verbose==1) // printing for non verbose
> - std::cout << i << " " << def << std::endl;
> - _prec.post(x); // postprocess preconditioner
> - res.iterations = i; // fill statistics
> - res.reduction = def/def0;
> - res.conv_rate = pow(res.reduction,1.0/i);
> - res.elapsed = watch.elapsed();
> - if (_verbose>0) // final print
> - std::cout << "=== rate=" << res.conv_rate
> - << ", T=" << res.elapsed
> - << ", TIT=" << res.elapsed/i
> - << ", IT=" << i << std::endl;
> - }
> + def = defnew; // update norm
> + if (def<def0*_reduction || def<1E-30 || i==_maxit) // convergence check
> + {
> + res.converged = true;
> + break;
> + }
>
> + // determine new search direction
> + q = 0; // clear correction
> + _prec.apply(q,b); // apply preconditioner
> + rho = _sp.dot(q,b); // orthogonalization
> + beta = rho/rholast; // scaling factor
> + p *= beta; // scale old search direction
> + p += q; // orthogonalization with correction
> + rholast = rho; // remember rho for recurrence
> + }
> +
> + if (_verbose==1) // printing for non verbose
> + this->printOutput(std::cout,i,def);
> +
> + _prec.post(x); // postprocess preconditioner
> + res.iterations = i; // fill statistics
> + res.reduction = def/def0;
> + res.conv_rate = pow(res.reduction,1.0/i);
> + res.elapsed = watch.elapsed();
> +
> + if (_verbose>0) // final print
> + {
> + std::cout << "=== rate=" << res.conv_rate
> + << ", T=" << res.elapsed
> + << ", TIT=" << res.elapsed/i
> + << ", IT=" << i << std::endl;
> + }
> + }
> +
> /*!
> \brief Apply inverse operator with given reduction factor.
>
> \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
> */
> - virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
> - {
> - _reduction = reduction;
> - (*this).apply(x,b,res);
> - }
> + virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
> + {
> + _reduction = reduction;
> + (*this).apply(x,b,res);
> + }
>
> private:
> - SeqScalarProduct<X> ssp;
> - LinearOperator<X,X>& _op;
> - Preconditioner<X,X>& _prec;
> - ScalarProduct<X>& _sp;
> - double _reduction;
> - int _maxit;
> - int _verbose;
> + SeqScalarProduct<X> ssp;
> + LinearOperator<X,X>& _op;
> + Preconditioner<X,X>& _prec;
> + ScalarProduct<X>& _sp;
> + double _reduction;
> + int _maxit;
> + int _verbose;
> };
>
>
> @@ -591,7 +638,7 @@
> */
> template<class L, class P>
> BiCGSTABSolver (L& op, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -604,7 +651,7 @@
> */
> template<class L, class S, class P>
> BiCGSTABSolver (L& op, S& sp, P& prec,
> - double reduction, int maxit, int verbose) :
> + double reduction, int maxit, int verbose) :
> _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> IsTrue< static_cast<int>(L::category) == static_cast<int>(P::category) >::yes();
> @@ -618,80 +665,83 @@
> */
> virtual void apply (X& x, X& b, InverseOperatorResult& res)
> {
> - const double EPSILON=1e-80;
> + const double EPSILON=1e-80;
>
> - int it;
> - field_type rho, rho_new, alpha, beta, h, omega;
> - field_type norm, norm_old, norm_0;
> + int it;
> + field_type rho, rho_new, alpha, beta, h, omega;
> + field_type norm, norm_old, norm_0;
>
> - //
> - // get vectors and matrix
> - //
> - X& r=b;
> - X p(x);
> - X v(x);
> - X t(x);
> - X y(x);
> - X rt(x);
> + //
> + // get vectors and matrix
> + //
> + X& r=b;
> + X p(x);
> + X v(x);
> + X t(x);
> + X y(x);
> + X rt(x);
>
> - //
> - // begin iteration
> - //
> + //
> + // begin iteration
> + //
>
> - // r = r - Ax; rt = r
> - res.clear(); // clear solver statistics
> - Timer watch; // start a timer
> - _op.applyscaleadd(-1,x,r); // overwrite b with defect
> - _prec.pre(x,r); // prepare preconditioner
> + // r = r - Ax; rt = r
> + res.clear(); // clear solver statistics
> + Timer watch; // start a timer
> + _op.applyscaleadd(-1,x,r); // overwrite b with defect
> + _prec.pre(x,r); // prepare preconditioner
>
> - rt=r;
> + rt=r;
>
> - norm = norm_old = norm_0 = _sp.norm(r);
> + norm = norm_old = norm_0 = _sp.norm(r);
>
> - p=0;
> - v=0;
> + p=0;
> + v=0;
>
> - rho = 1;
> - alpha = 1;
> - omega = 1;
> + rho = 1;
> + alpha = 1;
> + omega = 1;
>
> - if (_verbose>0) // printing
> - {
> - std::cout << "=== BiCGSTABSolver" << std::endl;
> - if (_verbose>1) {
> - std::cout << " Iter Defect Rate" << std::endl;
> - std::cout << "0 " << norm_0 << std::endl;
> - }
> - }
> + if (_verbose>0) // printing
> + {
> + std::cout << "=== BiCGSTABSolver" << std::endl;
> + if (_verbose>1)
> + {
> + this->printHeader(std::cout);
> + this->printOutput(std::cout,0,norm_0);
> + //std::cout << " Iter Defect Rate" << std::endl;
> + //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
> + }
> + }
>
> - if ( norm < (_reduction * norm_0) || norm<1E-30)
> - {
> - res.converged = 1;
> - _prec.post(x); // postprocess preconditioner
> - res.iterations = 0; // fill statistics
> - res.reduction = 0;
> - res.conv_rate = 0;
> - res.elapsed = watch.elapsed();
> - return;
> - }
> + if ( norm < (_reduction * norm_0) || norm<1E-30)
> + {
> + res.converged = 1;
> + _prec.post(x); // postprocess preconditioner
> + res.iterations = 0; // fill statistics
> + res.reduction = 0;
> + res.conv_rate = 0;
> + res.elapsed = watch.elapsed();
> + return;
> + }
>
> - //
> - // iteration
> - //
> + //
> + // iteration
> + //
>
> - it = 0;
> + it = 0;
>
> - while ( true )
> - {
> - //
> - // preprocess, set vecsizes etc.
> - //
> + while ( true )
> + {
> + //
> + // preprocess, set vecsizes etc.
> + //
>
> - // rho_new = < rt , r >
> - rho_new = _sp.dot(rt,r);
> + // rho_new = < rt , r >
> + rho_new = _sp.dot(rt,r);
>
> - // look if breakdown occured
> - if (std::abs(rho) <= EPSILON)
> + // look if breakdown occured
> + if (std::abs(rho) <= EPSILON)
> DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
> << rho << " <= EPSILON " << EPSILON
> << " after " << it << " iterations");
> @@ -701,109 +751,114 @@
> << " after " << it << " iterations");
>
>
> - if (it==0)
> - p = r;
> - else
> - {
> - beta = ( rho_new / rho ) * ( alpha / omega );
> - p.axpy(-omega,v); // p = r + beta (p - omega*v)
> - p *= beta;
> - p += r;
> - }
> + if (it==0)
> + p = r;
> + else
> + {
> + beta = ( rho_new / rho ) * ( alpha / omega );
> + p.axpy(-omega,v); // p = r + beta (p - omega*v)
> + p *= beta;
> + p += r;
> + }
>
> - // y = W^-1 * p
> - y = 0;
> - _prec.apply(y,p); // apply preconditioner
> + // y = W^-1 * p
> + y = 0;
> + _prec.apply(y,p); // apply preconditioner
>
> - // v = A * y
> - _op.apply(y,v);
> + // v = A * y
> + _op.apply(y,v);
>
> - // alpha = rho_new / < rt, v >
> - h = _sp.dot(rt,v);
> + // alpha = rho_new / < rt, v >
> + h = _sp.dot(rt,v);
>
> - if ( std::abs(h) < EPSILON )
> - DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
> -
> - alpha = rho_new / h;
> + if ( std::abs(h) < EPSILON )
> + DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
> +
> + alpha = rho_new / h;
>
> - // apply first correction to x
> - // x <- x + alpha y
> - x.axpy(alpha,y);
> + // apply first correction to x
> + // x <- x + alpha y
> + x.axpy(alpha,y);
>
> - // r = r - alpha*v
> - r.axpy(-alpha,v);
> + // r = r - alpha*v
> + r.axpy(-alpha,v);
>
> - //
> - // test stop criteria
> - //
> + //
> + // test stop criteria
> + //
>
> - it++;
> - norm = _sp.norm(r);
> + it++;
> + norm = _sp.norm(r);
>
> - if (_verbose>1) // print
> - std::cout << it << " " << norm << " " << norm/norm_old << std::endl;
> + if (_verbose>1) // print
> + {
> + this->printOutput(std::cout,it,norm,norm_old);
> + }
>
> - if ( norm < (_reduction * norm_0) )
> - {
> - res.converged = 1;
> - break;
> - }
> + if ( norm < (_reduction * norm_0) )
> + {
> + res.converged = 1;
> + break;
> + }
>
> - if (it >= _maxit)
> + if (it >= _maxit)
> break;
>
> - norm_old = norm;
> + norm_old = norm;
>
> - // y = W^-1 * r
> - y = 0;
> - _prec.apply(y,r);
> + // y = W^-1 * r
> + y = 0;
> + _prec.apply(y,r);
>
> - // t = A * y
> - _op.apply(y,t);
> -
> - // omega = < t, r > / < t, t >
> - omega = _sp.dot(t,r)/_sp.dot(t,t);
> + // t = A * y
> + _op.apply(y,t);
> +
> + // omega = < t, r > / < t, t >
> + omega = _sp.dot(t,r)/_sp.dot(t,t);
>
> - // apply second correction to x
> - // x <- x + omega y
> - x.axpy(omega,y);
> + // apply second correction to x
> + // x <- x + omega y
> + x.axpy(omega,y);
>
> - // r = s - omega*t (remember : r = s)
> - r.axpy(-omega,t);
> + // r = s - omega*t (remember : r = s)
> + r.axpy(-omega,t);
>
> - rho = rho_new;
> + rho = rho_new;
>
> - //
> - // test stop criteria
> - //
> + //
> + // test stop criteria
> + //
>
> - it++;
> + it++;
>
> - norm = _sp.norm(r);
> + norm = _sp.norm(r);
>
> - if (_verbose>1) // print
> - std::cout << it << " " << norm << " " << norm/norm_old << std::endl;
> + if (_verbose > 1) // print
> + {
> + this->printOutput(std::cout,it,norm,norm_old);
> + }
>
> - if ( norm < (_reduction * norm_0) || norm<1E-30)
> - {
> - res.converged = 1;
> - break;
> - }
> + if ( norm < (_reduction * norm_0) || norm<1E-30)
> + {
> + res.converged = 1;
> + break;
> + }
>
> - if (it >= _maxit)
> + if (it >= _maxit)
> break;
>
> - norm_old = norm;
> - }// while
> + norm_old = norm;
> + }// while
>
> - if (_verbose==1) // printing for non verbose
> - std::cout << it << " " << norm << std::endl;
> - _prec.post(x); // postprocess preconditioner
> - res.iterations = it; // fill statistics
> - res.reduction = norm/norm_0;
> - res.conv_rate = pow(res.reduction,1.0/it);
> - res.elapsed = watch.elapsed();
> - if (_verbose>0) // final print
> + if (_verbose==1) // printing for non verbose
> + this->printOutput(std::cout,it,norm);
> +
> + _prec.post(x); // postprocess preconditioner
> + res.iterations = it; // fill statistics
> + res.reduction = norm/norm_0;
> + res.conv_rate = pow(res.reduction,1.0/it);
> + res.elapsed = watch.elapsed();
> + if (_verbose>0) // final print
> std::cout << "=== rate=" << res.conv_rate
> << ", T=" << res.elapsed
> << ", TIT=" << res.elapsed/it
> @@ -822,13 +877,13 @@
> }
>
> private:
> - SeqScalarProduct<X> ssp;
> - LinearOperator<X,X>& _op;
> - Preconditioner<X,X>& _prec;
> - ScalarProduct<X>& _sp;
> - double _reduction;
> - int _maxit;
> - int _verbose;
> + SeqScalarProduct<X> ssp;
> + LinearOperator<X,X>& _op;
> + Preconditioner<X,X>& _prec;
> + ScalarProduct<X>& _sp;
> + double _reduction;
> + int _maxit;
> + int _verbose;
> };
>
>
>
>
> _______________________________________________
> Dune-Commit mailing list
> Dune-Commit at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-commit
>
--
************************************************************************
* Oliver Sander ** email: sander at mi.fu-berlin.de *
* Freie Universität Berlin ** phone: + 49 (30) 838 75217 *
* Institut für Mathematik II ** URL : page.mi.fu-berlin.de/~sander *
* Arnimallee 6 ** -------------------------------------*
* 14195 Berlin, Germany ** Member of MATHEON (www.matheon.de) *
************************************************************************
More information about the Dune
mailing list