From e5db27afdc9ab86ce71a9a00a9e230a4267c6daa Mon Sep 17 00:00:00 2001 From: Laurent Saboret Date: Mon, 25 Aug 2008 15:26:38 +0000 Subject: [PATCH] Fixed bug in Bi-CGSTAB solver when omega==0, which may happen for small linear systems. Added traces. --- OpenNL/include/CGAL/OpenNL/bicgstab.h | 90 ++++++++++++++++++--------- 1 file changed, 62 insertions(+), 28 deletions(-) diff --git a/OpenNL/include/CGAL/OpenNL/bicgstab.h b/OpenNL/include/CGAL/OpenNL/bicgstab.h index 82e3f953044..79fa0c8adef 100644 --- a/OpenNL/include/CGAL/OpenNL/bicgstab.h +++ b/OpenNL/include/CGAL/OpenNL/bicgstab.h @@ -26,10 +26,10 @@ * URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics * } * - * Laurent Saboret 2005-2006: Changes for CGAL: + * Laurent Saboret 2005-2008: Changes for CGAL: * - Added OpenNL namespace * - solve() returns true on success - * - test divisions by zero with IsZero() method + * - check divisions by zero * - added comments and traces * - copied BICGSTAB algorithm WITH preconditioner from Graphite 1.9 code */ @@ -88,6 +88,9 @@ public: // Solve the sparse linear system "A*x = b". Return true on success. bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) { +#ifdef DEBUG_TRACE + std::cerr << " Call BICGSTAB" << std::endl; +#endif CGAL_assertion(A.dimension() > 0); unsigned int n = A.dimension() ; // (Square) matrix dimension @@ -106,21 +109,20 @@ public: CoeffType rTh, rTAd, rTr, alpha, beta, omega, st, tt; unsigned int its=0; // Loop counter CoeffType err=epsilon_*epsilon_*BLAS::dot(b,b); // Error to reach - Vector r(n) ; // Current residue r=A*x-b + Vector r(n) ; // residue r=A*x-b mult(A,x,r); BLAS::axpy(-1,b,r); BLAS::copy(r,d); BLAS::copy(d,h); BLAS::copy(h,rT); - //CGAL_assertion( ! IsZero( BLAS::dot(rT,rT) ) ); + //CGAL_assertion(BLAS::dot(rT,rT) == 0.0); // may happen for small systems rTh=BLAS::dot(rT,h); // rTh = (rT|h) - rTr=BLAS::dot(r,r); // Current error rTr = (r|r) + rTr=BLAS::dot(r,r); // error rTr = (r|r) while ( rTr>err && its < max_iter) { mult(A,d,Ad); rTAd=BLAS::dot(rT,Ad); - if (IsZero(rTAd)) - break; // stop if bad conditioning + CGAL_assertion(rTAd != 0.0); alpha=rTh/rTAd; BLAS::axpy(-alpha,Ad,r); BLAS::copy(h,s); @@ -130,24 +132,36 @@ public: BLAS::scal(alpha,u); st=BLAS::dot(s,t); tt=BLAS::dot(t,t); - if ( IsZero(st) || IsZero(tt) ) + if (st == 0.0 || tt == 0.0) omega = 0 ; else omega = st/tt; BLAS::axpy(-omega,t,r); BLAS::axpy(-alpha,d,x); BLAS::axpy(-omega,s,x); + rTr=BLAS::dot(r,r); BLAS::copy(s,h); BLAS::axpy(-omega,t,h); - if (IsZero(omega)) - break; // stop if bad conditioning - if (IsZero(rTh)) - break; // stop if bad conditioning - beta=(alpha/omega)/rTh; rTh=BLAS::dot(rT,h); beta*=rTh; + if (omega == 0.0) // stop if omega==0 (success) + { +#ifdef DEBUG_TRACE + std::cerr << " BICGSTAB: omega=0" << std::endl; +#endif + break; + } + if (rTh == 0.0) // stop if rTh==0 (failure) + { +#ifdef DEBUG_TRACE + std::cerr << " BICGSTAB: rTh=0" << std::endl; +#endif + break; + } + beta=(alpha/omega)/rTh; // beta = (rTh/"old rTh") * (alpha/omega) + rTh=BLAS::dot(rT,h); + beta*=rTh; BLAS::scal(beta,d); BLAS::axpy(1,h,d); BLAS::axpy(-beta*omega,Ad,d); - rTr=BLAS::dot(r,r); ++its; } @@ -215,6 +229,9 @@ public: // Solve the sparse linear system "A*x = b". Return true on success. bool solve(const MATRIX &A, const PC_MATRIX &C, const VECTOR& b, VECTOR& x) { +#ifdef DEBUG_TRACE + std::cerr << " Call BICGSTAB with preconditioner" << std::endl; +#endif CGAL_assertion(A.dimension() > 0); unsigned int n = A.dimension() ; // (Square) matrix dimension @@ -234,22 +251,27 @@ public: CoeffType rTh, rTSd, rTr, alpha, beta, omega, st, tt; unsigned int its=0; // Loop counter CoeffType err=epsilon_*epsilon_*BLAS::dot(b,b); // Error to reach - Vector r(n) ; // Current residue r=A*x-b + Vector r(n) ; // residue r=A*x-b mult(A,x,r); BLAS::axpy(-1,b,r); mult(C,r,d); BLAS::copy(d,h); BLAS::copy(h,rT); - //CGAL_assertion( ! IsZero( BLAS::dot(rT,rT) ) ); + //CGAL_assertion(BLAS::dot(rT,rT) == 0.0); // may happen for small systems rTh=BLAS::dot(rT,h); // rTh = (rT|h) - rTr=BLAS::dot(r,r); // Current error rTr = (r|r) + rTr=BLAS::dot(r,r); // error rTr = (r|r) - while ( rTr>err && its < max_iter) { + while (rTr>err && its < max_iter) { mult(A,d,aux); mult(C,aux,Sd); rTSd=BLAS::dot(rT,Sd); - if (IsZero(rTSd)) - break; // stop if bad conditioning + if (rTSd == 0.0) // stop if rTSd==0 (failure) + { +#ifdef DEBUG_TRACE + std::cerr << " BICGSTAB with preconditioner: rTSd=0" << std::endl; +#endif + break; + } alpha=rTh/rTSd; BLAS::axpy(-alpha,aux,r); BLAS::copy(h,s); @@ -260,24 +282,36 @@ public: BLAS::scal(alpha,u); st=BLAS::dot(s,t); tt=BLAS::dot(t,t); - if ( IsZero(st) || IsZero(tt) ) + if (st == 0.0 || tt == 0.0) omega = 0 ; else omega = st/tt; BLAS::axpy(-omega,aux,r); BLAS::axpy(-alpha,d,x); BLAS::axpy(-omega,s,x); + rTr=BLAS::dot(r,r); BLAS::copy(s,h); BLAS::axpy(-omega,t,h); - if (IsZero(omega)) - break; // stop if bad conditioning - if (IsZero(rTh)) - break; // stop if bad conditioning - beta=(alpha/omega)/rTh; rTh=BLAS::dot(rT,h); beta*=rTh; - BLAS::scal(beta,d); + if (omega == 0.0) // stop if omega==0 (success) + { +#ifdef DEBUG_TRACE + std::cerr << " BICGSTAB with preconditioner: omega=0" << std::endl; +#endif + break; + } + if (rTh == 0.0) // stop if rTh==0 (failure) + { +#ifdef DEBUG_TRACE + std::cerr << " BICGSTAB with preconditioner: rTh=0" << std::endl; +#endif + break; + } + beta=(alpha/omega)/rTh; // beta = (rTh/"old rTh") * (alpha/omega) + rTh=BLAS::dot(rT,h); + beta*=rTh; + BLAS::scal(beta,d); // d = h + beta * (d - omega * Sd); BLAS::axpy(1,h,d); BLAS::axpy(-beta*omega,Sd,d); - rTr=BLAS::dot(r,r); ++its; }