Commit b3572ba7 by Praetorius, Simon

### Convert tabs to spaces

parent ff4baac0
 ... ... @@ -31,46 +31,46 @@ class basic_iteration template basic_iteration(const Vector& r0, int max_iter_, Real t, Real a = Real(0)) : error(0), i(0), my_norm_r0(std::abs(two_norm(r0))), max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) { } max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) { } /// Constructor basic_iteration(Real nb, int max_iter_, Real t, Real a = Real(0)) : error(0), i(0), my_norm_r0(nb), max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) {} my_quite(false), my_suppress(false) {} virtual ~basic_iteration() {} bool check_max() { if (i >= max_iter) error= 1, is_finished= true, err_msg= "Too many iterations."; return is_finished; if (i >= max_iter) error= 1, is_finished= true, err_msg= "Too many iterations."; return is_finished; } /// Iteration finished according to the norm of r template bool finished(const Vector& r) { if (converged(two_norm(r))) return is_finished= true; return check_max(); if (converged(two_norm(r))) return is_finished= true; return check_max(); } /// Iteration finished according to residual value r bool finished(const Real& r) { if (converged(r)) return is_finished= true; return check_max(); if (converged(r)) return is_finished= true; return check_max(); } /// Iteration finished according to complex residual value r template bool finished(const std::complex& r) { if (converged(std::abs(r))) return is_finished= true; return check_max(); if (converged(std::abs(r))) return is_finished= true; return check_max(); } /// Iteration finished according to last provided residual ... ... @@ -84,12 +84,12 @@ class basic_iteration bool converged() const { // std::cout << "abs: " << resid_ << " <= " << atol_ << '\n'; // std::cout << "rel: " << resid_ << " <= " << rtol_ << " * " << my_norm_r0 // << " = " << rtol_ * my_norm_r0 << '\n'; if (my_norm_r0 == 0) return resid_ <= atol_; // ignore relative tolerance if |r0| is zero return resid_ <= rtol_ * my_norm_r0 || resid_ <= atol_; // std::cout << "abs: " << resid_ << " <= " << atol_ << '\n'; // std::cout << "rel: " << resid_ << " <= " << rtol_ << " * " << my_norm_r0 // << " = " << rtol_ * my_norm_r0 << '\n'; if (my_norm_r0 == 0) return resid_ <= atol_; // ignore relative tolerance if |r0| is zero return resid_ <= rtol_ * my_norm_r0 || resid_ <= atol_; } public: ... ... @@ -144,14 +144,14 @@ class basic_iteration void update_progress(const basic_iteration& that) { i= that.i; resid_= that.resid_; if (that.error > 1) { // copy error except too many iterations error= that.error; err_msg= that.err_msg; is_finished= true; } else finished(resid_); i= that.i; resid_= that.resid_; if (that.error > 1) { // copy error except too many iterations error= that.error; err_msg= that.err_msg; is_finished= true; } else finished(resid_); } protected: ... ...
 ... ... @@ -28,15 +28,15 @@ namespace itl { void print_resid() { if (!this->my_quite && this->i % cycle == 0) if (multi_print || this->i != last_print) { // Avoid multiple print-outs in same iteration out << "iteration " << this->i << ": resid " << this->resid() if (!this->my_quite && this->i % cycle == 0) if (multi_print || this->i != last_print) { // Avoid multiple print-outs in same iteration out << "iteration " << this->i << ": resid " << this->resid() # ifdef MTL_VERBOSE_ITERATION << " / " << this->norm_r0 << " = " << this->resid() / this->norm_r0 << " (rel. error)" << " / " << this->norm_r0 << " = " << this->resid() / this->norm_r0 << " (rel. error)" # endif << std::endl; last_print= this->i; } << std::endl; last_print= this->i; } } public: ... ... @@ -44,14 +44,14 @@ namespace itl { /// Constructor template cyclic_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100, OStream& out = std::cout) : super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out) OStream& out = std::cout) : super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out) {} /// Constructor cyclic_iteration(Real r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100, OStream& out = std::cout) : super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out) OStream& out = std::cout) : super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out) {} ... ... @@ -62,9 +62,9 @@ namespace itl { template bool finished(const T& r) { bool ret= super::finished(r); print_resid(); return ret; bool ret= super::finished(r); print_resid(); return ret; } inline self& operator++() { ++this->i; return *this; } ///< Increment counter ... ... @@ -84,14 +84,14 @@ namespace itl { /// Error code with final resume int error_code() const { if (!this->my_suppress) out << "finished! error code = " << this->error << '\n' << this->iterations() << " iterations\n" << this->resid() << " is actual final residual. \n" << this->relresid() << " is actual relative tolerance achieved. \n" << "Relative tol: " << this->rtol_ << " Absolute tol: " << this->atol_ << '\n' << "Convergence: " << pow(this->relresid(), 1.0 / double(this->iterations())) << std::endl; return this->error; if (!this->my_suppress) out << "finished! error code = " << this->error << '\n' << this->iterations() << " iterations\n" << this->resid() << " is actual final residual. \n" << this->relresid() << " is actual relative tolerance achieved. \n" << "Relative tol: " << this->rtol_ << " Absolute tol: " << this->atol_ << '\n' << "Convergence: " << pow(this->relresid(), 1.0 / double(this->iterations())) << std::endl; return this->error; } protected: int cycle, last_print; ... ...
 ... ... @@ -24,8 +24,8 @@ namespace itl { public: template noisy_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0), OStream& out = std::cout) : cyclic_iteration(r0, max_iter_, tol_, atol_, 1, out) OStream& out = std::cout) : cyclic_iteration(r0, max_iter_, tol_, atol_, 1, out) {} }; ... ...
 ... ... @@ -25,61 +25,61 @@ namespace itl { namespace pc { template struct solver; template class identity; // template Vector solve(const identity&, const Vector& x); // template Vector adjoint_solve(const identity&, const Vector& x); template struct solver; template class identity; // template Vector solve(const identity&, const Vector& x); // template Vector adjoint_solve(const identity&, const Vector& x); template class diagonal; // template Vector solve(const diagonal& P, const Vector& x); // template Vector adjoint_solve(const diagonal& P, const Vector& x); template class diagonal; // template Vector solve(const diagonal& P, const Vector& x); // template Vector adjoint_solve(const diagonal& P, const Vector& x); template class ilu; template class ilu_0; // Maybe we should declare the default here??? template class ilut; // Maybe we should declare the default here??? // template Vector solve(const ilu_0& P, const Vector& x); // template Vector adjoint_solve(const ilu_0& P, const Vector& x); template class ilu; template class ilu_0; // Maybe we should declare the default here??? template class ilut; // Maybe we should declare the default here??? // template Vector solve(const ilu_0& P, const Vector& x); // template Vector adjoint_solve(const ilu_0& P, const Vector& x); template class ic_0; // Maybe we should declare the default here??? // template Vector solve(const ic_0& P, const Vector& x); // template Vector adjoint_solve(const ic_0& P, const Vector& x); template class ic_0; // Maybe we should declare the default here??? // template Vector solve(const ic_0& P, const Vector& x); // template Vector adjoint_solve(const ic_0& P, const Vector& x); } // namespace pc template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB, typename Preconditioner, typename Iteration > typename Preconditioner, typename Iteration > int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, const Preconditioner& M, Iteration& iter); const Preconditioner& M, Iteration& iter); template < typename LinearOperator, typename Preconditioner= pc::identity, typename RightPreconditioner= pc::identity > typename RightPreconditioner= pc::identity > class cg_solver; template < typename LinearOperator, typename Vector, typename Preconditioner, typename Iteration > typename Preconditioner, typename Iteration > int bicg(const LinearOperator &A, Vector &x, const Vector &b, const Preconditioner &M, Iteration& iter); const Preconditioner &M, Iteration& iter); template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner, class Iteration > class Preconditioner, class Iteration > int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, const Preconditioner& M, Iteration& iter); const Preconditioner& M, Iteration& iter); template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner, class Iteration > class Preconditioner, class Iteration > int bicgstab_2(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, const Preconditioner& M, Iteration& iter); const Preconditioner& M, Iteration& iter); template < typename LinearOperator, typename Vector, typename LeftPreconditioner, typename RightPreconditioner, typename Iteration > typename LeftPreconditioner, typename RightPreconditioner, typename Iteration > int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b, const LeftPreconditioner &L, const RightPreconditioner &R, Iteration& iter, size_t l); const LeftPreconditioner &L, const RightPreconditioner &R, Iteration& iter, size_t l); template < typename LinearOperator, typename Preconditioner= pc::identity, typename RightPreconditioner= pc::identity > typename RightPreconditioner= pc::identity > class bicgstab_ell_solver; template ... ...
 ... ... @@ -35,37 +35,37 @@ struct base_solver template int step(HilbertSpaceX& x, const HilbertSpaceB& b) const { itl::basic_iteration iter(b, 1, 0, 0); return static_cast(this)->solve(x, b, iter); itl::basic_iteration iter(b, 1, 0, 0); return static_cast(this)->solve(x, b, iter); } /// Solve using the iteration (resets initial residuum) template int operator()(HilbertSpaceX& x, const HilbertSpaceB& b) { my_iteration.restart(); if (two_norm(x) == 0) my_iteration.set_norm_r0(two_norm(b)); else { HilbertSpaceB r(b); r-= A * x; my_iteration.set_norm_r0(two_norm(r)); } return static_cast(this)->solve(x, b, my_iteration); my_iteration.restart(); if (two_norm(x) == 0) my_iteration.set_norm_r0(two_norm(b)); else { HilbertSpaceB r(b); r-= A * x; my_iteration.set_norm_r0(two_norm(r)); } return static_cast(this)->solve(x, b, my_iteration); } /// Set log level to \p level /** Concrete meaning of log level depends on specific solver. The general understanding is: - Level 0: no logging at all - Level 1: print final log (convergence and alike) - Level 2: print log on iterations - Level 3 or higher: treated as 2 **/ The general understanding is: - Level 0: no logging at all - Level 1: print final log (convergence and alike) - Level 2: print log on iterations - Level 3 or higher: treated as 2 **/ void set_log_level(unsigned level) { this->level= level; my_iteration.set_quite(level < 2); my_iteration.suppress_resume(level == 0); this->level= level; my_iteration.set_quite(level < 2); my_iteration.suppress_resume(level == 0); } /// Get log level ... ...
 ... ... @@ -26,42 +26,42 @@ namespace itl { /// Bi-Conjugate Gradient template < typename LinearOperator, typename Vector, typename Preconditioner, typename Iteration > typename Preconditioner, typename Iteration > int bicg(const LinearOperator &A, Vector &x, const Vector &b, const Preconditioner &M, Iteration& iter) const Preconditioner &M, Iteration& iter) { mtl::vampir_trace<7003> tracer; using mtl::conj; typedef typename mtl::Collection::value_type Scalar; Scalar rho_1(0), rho_2(0), alpha(0), beta(0); Vector r(b - A * x), z(resource(x)), p(resource(x)), q(resource(x)), r_tilde(r), z_tilde(resource(x)), p_tilde(resource(x)), q_tilde(resource(x)); r_tilde(r), z_tilde(resource(x)), p_tilde(resource(x)), q_tilde(resource(x)); while ( ! iter.finished(r)) { ++iter; z= solve(M, r); z_tilde= adjoint_solve(M, r_tilde); rho_1= dot(z_tilde, z); if (rho_1 == 0.) return iter.fail(2, "bicg breakdown"); if (iter.first()) { p= z; p_tilde= z_tilde; } else { beta= rho_1 / rho_2; p= z + beta * p; p_tilde= z_tilde + conj(beta) * p_tilde; } q= A * p; q_tilde= adjoint(A) * p_tilde; alpha= rho_1 / dot(p_tilde, q); x+= alpha * p; r-= alpha * q; r_tilde-= conj(alpha) * q_tilde; rho_2= rho_1; ++iter; z= solve(M, r); z_tilde= adjoint_solve(M, r_tilde); rho_1= dot(z_tilde, z); if (rho_1 == 0.) return iter.fail(2, "bicg breakdown"); if (iter.first()) { p= z; p_tilde= z_tilde; } else { beta= rho_1 / rho_2; p= z + beta * p; p_tilde= z_tilde + conj(beta) * p_tilde; } q= A * p; q_tilde= adjoint(A) * p_tilde; alpha= rho_1 / dot(p_tilde, q); x+= alpha * p; r-= alpha * q; r_tilde-= conj(alpha) * q_tilde; rho_2= rho_1; } return iter; } ... ... @@ -69,7 +69,7 @@ int bicg(const LinearOperator &A, Vector &x, const Vector &b, /// Solver class for BiCG method; right preconditioner ignored (prints warning if not identity) /** Methods inherited from \ref base_solver. **/ template < typename LinearOperator, typename Preconditioner= pc::identity, typename RightPreconditioner= pc::identity > typename RightPreconditioner= pc::identity > class bicg_solver : public base_solver< bicg_solver, LinearOperator > { ... ... @@ -78,22 +78,22 @@ class bicg_solver /// Construct solver from a linear operator; generate (left) preconditioner from it explicit bicg_solver(const LinearOperator& A) : base(A), L(A) { if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; } /// Construct solver from a linear operator and (left) preconditioner bicg_solver(const LinearOperator& A, const Preconditioner& L) : base(A), L(L) { if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; } /// Solve linear system approximately as specified by \p iter template < typename HilbertSpaceX, typename HilbertSpaceB, typename Iteration > int solve(HilbertSpaceX& x, const HilbertSpaceB& b, Iteration& iter) const { return bicg(this->A, x, b, L, iter); return bicg(this->A, x, b, L, iter); } private: ... ...
 ... ... @@ -25,9 +25,9 @@ namespace itl { /// Bi-Conjugate Gradient Stabilized template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner, class Iteration > class Preconditioner, class Iteration > int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, const Preconditioner& M, Iteration& iter) const Preconditioner& M, Iteration& iter) { typedef typename mtl::Collection::value_type Scalar; typedef HilbertSpaceX Vector; ... ... @@ -80,7 +80,7 @@ int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, /// Solver class for BiCGStab method; right preconditioner ignored (prints warning if not identity) /** Methods inherited from \ref base_solver. **/ template < typename LinearOperator, typename Preconditioner= pc::identity, typename RightPreconditioner= pc::identity > typename RightPreconditioner= pc::identity > class bicgstab_solver : public base_solver< bicgstab_solver, LinearOperator > { ... ... @@ -89,22 +89,22 @@ class bicgstab_solver /// Construct solver from a linear operator; generate (left) preconditioner from it explicit bicgstab_solver(const LinearOperator& A) : base(A), L(A) { if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; } /// Construct solver from a linear operator and (left) preconditioner bicgstab_solver(const LinearOperator& A, const Preconditioner& L) : base(A), L(L) { if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; if (!pc::static_is_identity::value) std::cerr << "Right Preconditioner ignored!" << std::endl; } /// Solve linear system approximately as specified by \p iter template < typename HilbertSpaceX, typename HilbertSpaceB, typename Iteration > int solve(HilbertSpaceX& x, const HilbertSpaceB& b, Iteration& iter) const { return bicgstab(this->A, x, b, L, iter); return bicgstab(this->A, x, b, L, iter); } private: ... ...
 ... ... @@ -33,9 +33,9 @@ namespace itl { /// Bi-Conjugate Gradient Stabilized(2) template < typename LinearOperator, typename Vector, typename Preconditioner, typename Iteration > typename Preconditioner, typename Iteration > int bicgstab_2(const LinearOperator &A, Vector &x, const Vector &b, const Preconditioner &L, Iteration& iter) const Preconditioner &L, Iteration& iter) { mtl::vampir_trace<7005> tracer; using mtl::size; using mtl::irange; using mtl::imax; using mtl::mat::strict_upper; using mtl::lazy; ... ... @@ -53,9 +53,9 @@ int bicgstab_2(const LinearOperator &A, Vector &x, const Vector &b, x0= zero; r_hat[0]= b; if (two_norm(x) != zero) { r_hat[0]-= A * x; x0= x; x= zero; r_hat[0]-= A * x; x0= x; x= zero; } Vector r0_tilde(r_hat[0]/two_norm(r_hat[0])); ... ... @@ -68,62 +68,62 @@ int bicgstab_2(const LinearOperator &A, Vector &x, const Vector &b, mtl::dense_vector sigma(l+1), gamma(l+1), gamma_a(l+1), gamma_aa(l+1); while (! iter.finished(r_hat[0])) { ++iter; rho_0= -omega * rho_0; ++iter; rho_0= -omega * rho_0; for (Size j= 0; j < 2; ++j) { rho_1= dot(r0_tilde, r_hat[j]); beta= alpha * rho_1/rho_0; rho_0= rho_1; for (Size j= 0; j < 2; ++j) { rho_1= dot(r0_tilde, r_hat[j]); beta= alpha * rho_1/rho_0; rho_0= rho_1; for (Size i= 0; i <= j; ++i) u_hat[i]= r_hat[i] - beta * u_hat[i]; for (Size i= 0; i <= j; ++i) u_hat[i]= r_hat[i] - beta * u_hat[i]; y= A * u_hat[j]; u_hat[j+1]= solve(L, y); Gamma= dot(r0_tilde, u_hat[j+1]); alpha= rho_0 / Gamma; y= A * u_hat[j]; u_hat[j+1]= solve(L, y); Gamma= dot(r0_tilde, u_hat[j+1]); alpha= rho_0 / Gamma; for (Size i= 0; i <= j; ++i) r_hat[i]-= alpha * u_hat[i+1]; for (Size i= 0; i <= j; ++i) r_hat[i]-= alpha * u_hat[i+1]; if (iter.finished(r_hat[j])) { x+= x0; return iter; } y= A * r_hat[j]; r_hat[j+1]= solve(L, y); x+= alpha * u_hat[0]; } // mod GS (MR part) irange i1m(1, imax); mtl::dense_vector r_hat_tail(r_hat[i1m]); tau[i1m][i1m]= orthogonalize_factors(r_hat_tail); for (Size j= 1; j <= l; ++j) gamma_a[j]= dot(r_hat[j], r_hat[0]) / tau[j][j];