From 456741cb265ceb8ca3cf92559faf96f4f08e94b8 Mon Sep 17 00:00:00 2001 From: Laurent Saboret Date: Mon, 4 Dec 2006 08:29:51 +0000 Subject: [PATCH] Added comments --- OpenNL/include/OpenNL/bicgstab.h | 29 +++++++-------- OpenNL/include/OpenNL/conjugate_gradient.h | 41 +++++++++++----------- 2 files changed, 35 insertions(+), 35 deletions(-) diff --git a/OpenNL/include/OpenNL/bicgstab.h b/OpenNL/include/OpenNL/bicgstab.h index d12649e023a..0b284fb3abd 100644 --- a/OpenNL/include/OpenNL/bicgstab.h +++ b/OpenNL/include/OpenNL/bicgstab.h @@ -26,7 +26,7 @@ * URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics * } * - * Laurent Saboret 01/2005-04/2005: Changes for CGAL: + * Laurent Saboret 2005-2006: Changes for CGAL: * - Added OpenNL namespace * - solve() returns true on success * - test divisions by zero with IsZero() method @@ -63,44 +63,45 @@ namespace OpenNL { * by Christian Badura, available from: * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ * - * @param A generic matrix, a function + * @param A generic square matrix; a function * mult(const MATRIX& M, const double* x, double* y) - * needs to be defined. + * and a member function + * int dimension() const + * must to be defined. * @param b right hand side of the system. * @param x initial value. * @param eps threshold for the residual. * @param max_iter maximum number of iterations. */ -template < class MATRIX, class VECTOR> class Solver_BICGSTAB { +template class Solver_BICGSTAB { public: typedef MATRIX Matrix ; typedef VECTOR Vector ; typedef typename Vector::CoeffType CoeffType ; +public: Solver_BICGSTAB() { epsilon_ = 1e-6 ; max_iter_ = 0 ; } + // Default copy constructor, operator =() and destructor are fine + void set_epsilon(CoeffType eps) { epsilon_ = eps ; } void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } // Solve the sparse linear system "A*x = b". Return true on success. - // - // Preconditions: - // * A.dimension() == b.dimension() - // * A.dimension() == x.dimension() - bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) { - assert(A.dimension() == b.dimension()) ; - assert(A.dimension() == x.dimension()) ; + bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) + { assert (A.dimension() > 0); - unsigned int n = A.dimension() ; // (Square) matrix dimension + unsigned int max_iter = max_iter_ ; // Max number of iterations if(max_iter == 0) { max_iter = 5 * n ; } + Vector rT(n) ; // Initial residue rT=Ax-b Vector d(n) ; Vector h(n) ; @@ -154,7 +155,7 @@ public: BLAS::axpy(1,h,d); BLAS::axpy(-beta*omega,Ad,d); rTr=BLAS::dot(r,r); - its++ ; + ++its; } bool success = (rTr <= err); @@ -171,7 +172,7 @@ public: private: // Test if a floating point number is (close to) 0.0 - static inline bool IsZero(CoeffType a) + static bool IsZero(CoeffType a) { return (std::fabs(a) < 10.0 * (std::numeric_limits::min)()); } diff --git a/OpenNL/include/OpenNL/conjugate_gradient.h b/OpenNL/include/OpenNL/conjugate_gradient.h index 66b70194625..d1383c1209c 100644 --- a/OpenNL/include/OpenNL/conjugate_gradient.h +++ b/OpenNL/include/OpenNL/conjugate_gradient.h @@ -26,11 +26,11 @@ * URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics * } * - * Laurent Saboret 01/2005-04/2005: Changes for CGAL: + * Laurent Saboret 2005-2006: Changes for CGAL: * - Added OpenNL namespace * - solve() returns true on success * - test divisions by zero with IsZero() method - * - added comments and traces + * - added comments */ #ifndef __OPENNL_CONJUGATE_GRADIENT__ @@ -63,46 +63,46 @@ namespace OpenNL { * by Christian Badura, available from: * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ * - * @param A generic matrix, a function + * @param A generic square matrix; a function * mult(const MATRIX& M, const double* x, double* y) - * needs to be defined. + * and a member function + * int dimension() const + * must to be defined. * @param b right hand side of the system. * @param x initial value. * @param eps threshold for the residual. * @param max_iter maximum number of iterations. */ -template< class MATRIX, class VECTOR > class Solver_CG { +template class Solver_CG { public: typedef MATRIX Matrix ; typedef VECTOR Vector ; typedef typename Vector::CoeffType CoeffType ; +public: Solver_CG() { epsilon_ = 1e-6 ; max_iter_ = 0 ; } + // Default copy constructor, operator =() and destructor are fine + void set_epsilon(CoeffType eps) { epsilon_ = eps ; } void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } - // Solve the sparse linear system "A*x = b" for A SYMMETRIC POSITIVE + // Solve the sparse linear system "A*x = b" for A symmetric positive definite // Return true on success - // - // Preconditions: - // * A.dimension() == b.dimension() - // * A.dimension() == x.dimension() - bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) { - assert(A.dimension() == b.dimension()) ; - assert(A.dimension() == x.dimension()) ; - assert (A.dimension() > 0); - + bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) + { + assert(A.dimension() > 0); unsigned int n = A.dimension() ; // (Square) matrix dimension unsigned int max_iter = max_iter_ ; // Max number of iterations if(max_iter == 0) { max_iter = 5 * n ; } + Vector g(n) ; Vector r(n) ; Vector p(n) ; @@ -116,10 +116,9 @@ public: BLAS::scal(-1,g); // Initially, r=g=b-A*x BLAS::copy(g,r); // r = g - CoeffType gg=BLAS::dot(g,g); // Current error gg = (g|g) - while ( gg>err && its < max_iter) - { + CoeffType gg=BLAS::dot(g,g); // Current error gg = (g|g) + while ( gg>err && its < max_iter) { mult(A,r,p); rho=BLAS::dot(p,p); sig=BLAS::dot(r,p); @@ -134,8 +133,8 @@ public: gam=(t*t*rho-tau)/tau; BLAS::scal(gam,r); BLAS::axpy(1,g,r); - gg=BLAS::dot(g,g); // Current error gg = (g|g) - its++; + gg=BLAS::dot(g,g); // Update error gg = (g|g) + ++its; } bool success = (gg <= err); @@ -152,7 +151,7 @@ public: private: // Test if a floating point number is (close to) 0.0 - static inline bool IsZero(CoeffType a) + static bool IsZero(CoeffType a) { return (fabs(a) < 10.0 * (std::numeric_limits::min)()); }