mirror of https://github.com/CGAL/cgal
Fixed bug in Bi-CGSTAB solver when omega==0, which may happen for small linear systems.
Added traces.
This commit is contained in:
Laurent Saboret
parent
f524fb618e
commit
e5db27afdc
|
|
@ -26,10 +26,10 @@
|
||||||
* URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics
|
* 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
|
* - Added OpenNL namespace
|
||||||
* - solve() returns true on success
|
* - solve() returns true on success
|
||||||
* - test divisions by zero with IsZero() method
|
* - check divisions by zero
|
||||||
* - added comments and traces
|
* - added comments and traces
|
||||||
* - copied BICGSTAB algorithm WITH preconditioner from Graphite 1.9 code
|
* - 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.
|
// Solve the sparse linear system "A*x = b". Return true on success.
|
||||||
bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x)
|
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);
|
CGAL_assertion(A.dimension() > 0);
|
||||||
unsigned int n = A.dimension() ; // (Square) matrix dimension
|
unsigned int n = A.dimension() ; // (Square) matrix dimension
|
||||||
|
|
||||||
|
|
@ -106,21 +109,20 @@ public:
|
||||||
CoeffType rTh, rTAd, rTr, alpha, beta, omega, st, tt;
|
CoeffType rTh, rTAd, rTr, alpha, beta, omega, st, tt;
|
||||||
unsigned int its=0; // Loop counter
|
unsigned int its=0; // Loop counter
|
||||||
CoeffType err=epsilon_*epsilon_*BLAS<Vector>::dot(b,b); // Error to reach
|
CoeffType err=epsilon_*epsilon_*BLAS<Vector>::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);
|
mult(A,x,r);
|
||||||
BLAS<Vector>::axpy(-1,b,r);
|
BLAS<Vector>::axpy(-1,b,r);
|
||||||
BLAS<Vector>::copy(r,d);
|
BLAS<Vector>::copy(r,d);
|
||||||
BLAS<Vector>::copy(d,h);
|
BLAS<Vector>::copy(d,h);
|
||||||
BLAS<Vector>::copy(h,rT);
|
BLAS<Vector>::copy(h,rT);
|
||||||
//CGAL_assertion( ! IsZero( BLAS<Vector>::dot(rT,rT) ) );
|
//CGAL_assertion(BLAS<Vector>::dot(rT,rT) == 0.0); // may happen for small systems
|
||||||
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
|
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
|
||||||
rTr=BLAS<Vector>::dot(r,r); // Current error rTr = (r|r)
|
rTr=BLAS<Vector>::dot(r,r); // error rTr = (r|r)
|
||||||
|
|
||||||
while ( rTr>err && its < max_iter) {
|
while ( rTr>err && its < max_iter) {
|
||||||
mult(A,d,Ad);
|
mult(A,d,Ad);
|
||||||
rTAd=BLAS<Vector>::dot(rT,Ad);
|
rTAd=BLAS<Vector>::dot(rT,Ad);
|
||||||
if (IsZero(rTAd))
|
CGAL_assertion(rTAd != 0.0);
|
||||||
break; // stop if bad conditioning
|
|
||||||
alpha=rTh/rTAd;
|
alpha=rTh/rTAd;
|
||||||
BLAS<Vector>::axpy(-alpha,Ad,r);
|
BLAS<Vector>::axpy(-alpha,Ad,r);
|
||||||
BLAS<Vector>::copy(h,s);
|
BLAS<Vector>::copy(h,s);
|
||||||
|
|
@ -130,24 +132,36 @@ public:
|
||||||
BLAS<Vector>::scal(alpha,u);
|
BLAS<Vector>::scal(alpha,u);
|
||||||
st=BLAS<Vector>::dot(s,t);
|
st=BLAS<Vector>::dot(s,t);
|
||||||
tt=BLAS<Vector>::dot(t,t);
|
tt=BLAS<Vector>::dot(t,t);
|
||||||
if ( IsZero(st) || IsZero(tt) )
|
if (st == 0.0 || tt == 0.0)
|
||||||
omega = 0 ;
|
omega = 0 ;
|
||||||
else
|
else
|
||||||
omega = st/tt;
|
omega = st/tt;
|
||||||
BLAS<Vector>::axpy(-omega,t,r);
|
BLAS<Vector>::axpy(-omega,t,r);
|
||||||
BLAS<Vector>::axpy(-alpha,d,x);
|
BLAS<Vector>::axpy(-alpha,d,x);
|
||||||
BLAS<Vector>::axpy(-omega,s,x);
|
BLAS<Vector>::axpy(-omega,s,x);
|
||||||
|
rTr=BLAS<Vector>::dot(r,r);
|
||||||
BLAS<Vector>::copy(s,h);
|
BLAS<Vector>::copy(s,h);
|
||||||
BLAS<Vector>::axpy(-omega,t,h);
|
BLAS<Vector>::axpy(-omega,t,h);
|
||||||
if (IsZero(omega))
|
if (omega == 0.0) // stop if omega==0 (success)
|
||||||
break; // stop if bad conditioning
|
{
|
||||||
if (IsZero(rTh))
|
#ifdef DEBUG_TRACE
|
||||||
break; // stop if bad conditioning
|
std::cerr << " BICGSTAB: omega=0" << std::endl;
|
||||||
beta=(alpha/omega)/rTh; rTh=BLAS<Vector>::dot(rT,h); beta*=rTh;
|
#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<Vector>::dot(rT,h);
|
||||||
|
beta*=rTh;
|
||||||
BLAS<Vector>::scal(beta,d);
|
BLAS<Vector>::scal(beta,d);
|
||||||
BLAS<Vector>::axpy(1,h,d);
|
BLAS<Vector>::axpy(1,h,d);
|
||||||
BLAS<Vector>::axpy(-beta*omega,Ad,d);
|
BLAS<Vector>::axpy(-beta*omega,Ad,d);
|
||||||
rTr=BLAS<Vector>::dot(r,r);
|
|
||||||
++its;
|
++its;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
@ -215,6 +229,9 @@ public:
|
||||||
// Solve the sparse linear system "A*x = b". Return true on success.
|
// 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)
|
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);
|
CGAL_assertion(A.dimension() > 0);
|
||||||
unsigned int n = A.dimension() ; // (Square) matrix dimension
|
unsigned int n = A.dimension() ; // (Square) matrix dimension
|
||||||
|
|
||||||
|
|
@ -234,22 +251,27 @@ public:
|
||||||
CoeffType rTh, rTSd, rTr, alpha, beta, omega, st, tt;
|
CoeffType rTh, rTSd, rTr, alpha, beta, omega, st, tt;
|
||||||
unsigned int its=0; // Loop counter
|
unsigned int its=0; // Loop counter
|
||||||
CoeffType err=epsilon_*epsilon_*BLAS<Vector>::dot(b,b); // Error to reach
|
CoeffType err=epsilon_*epsilon_*BLAS<Vector>::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);
|
mult(A,x,r);
|
||||||
BLAS<Vector>::axpy(-1,b,r);
|
BLAS<Vector>::axpy(-1,b,r);
|
||||||
mult(C,r,d);
|
mult(C,r,d);
|
||||||
BLAS<Vector>::copy(d,h);
|
BLAS<Vector>::copy(d,h);
|
||||||
BLAS<Vector>::copy(h,rT);
|
BLAS<Vector>::copy(h,rT);
|
||||||
//CGAL_assertion( ! IsZero( BLAS<Vector>::dot(rT,rT) ) );
|
//CGAL_assertion(BLAS<Vector>::dot(rT,rT) == 0.0); // may happen for small systems
|
||||||
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
|
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
|
||||||
rTr=BLAS<Vector>::dot(r,r); // Current error rTr = (r|r)
|
rTr=BLAS<Vector>::dot(r,r); // error rTr = (r|r)
|
||||||
|
|
||||||
while ( rTr>err && its < max_iter) {
|
while (rTr>err && its < max_iter) {
|
||||||
mult(A,d,aux);
|
mult(A,d,aux);
|
||||||
mult(C,aux,Sd);
|
mult(C,aux,Sd);
|
||||||
rTSd=BLAS<Vector>::dot(rT,Sd);
|
rTSd=BLAS<Vector>::dot(rT,Sd);
|
||||||
if (IsZero(rTSd))
|
if (rTSd == 0.0) // stop if rTSd==0 (failure)
|
||||||
break; // stop if bad conditioning
|
{
|
||||||
|
#ifdef DEBUG_TRACE
|
||||||
|
std::cerr << " BICGSTAB with preconditioner: rTSd=0" << std::endl;
|
||||||
|
#endif
|
||||||
|
break;
|
||||||
|
}
|
||||||
alpha=rTh/rTSd;
|
alpha=rTh/rTSd;
|
||||||
BLAS<Vector>::axpy(-alpha,aux,r);
|
BLAS<Vector>::axpy(-alpha,aux,r);
|
||||||
BLAS<Vector>::copy(h,s);
|
BLAS<Vector>::copy(h,s);
|
||||||
|
|
@ -260,24 +282,36 @@ public:
|
||||||
BLAS<Vector>::scal(alpha,u);
|
BLAS<Vector>::scal(alpha,u);
|
||||||
st=BLAS<Vector>::dot(s,t);
|
st=BLAS<Vector>::dot(s,t);
|
||||||
tt=BLAS<Vector>::dot(t,t);
|
tt=BLAS<Vector>::dot(t,t);
|
||||||
if ( IsZero(st) || IsZero(tt) )
|
if (st == 0.0 || tt == 0.0)
|
||||||
omega = 0 ;
|
omega = 0 ;
|
||||||
else
|
else
|
||||||
omega = st/tt;
|
omega = st/tt;
|
||||||
BLAS<Vector>::axpy(-omega,aux,r);
|
BLAS<Vector>::axpy(-omega,aux,r);
|
||||||
BLAS<Vector>::axpy(-alpha,d,x);
|
BLAS<Vector>::axpy(-alpha,d,x);
|
||||||
BLAS<Vector>::axpy(-omega,s,x);
|
BLAS<Vector>::axpy(-omega,s,x);
|
||||||
|
rTr=BLAS<Vector>::dot(r,r);
|
||||||
BLAS<Vector>::copy(s,h);
|
BLAS<Vector>::copy(s,h);
|
||||||
BLAS<Vector>::axpy(-omega,t,h);
|
BLAS<Vector>::axpy(-omega,t,h);
|
||||||
if (IsZero(omega))
|
if (omega == 0.0) // stop if omega==0 (success)
|
||||||
break; // stop if bad conditioning
|
{
|
||||||
if (IsZero(rTh))
|
#ifdef DEBUG_TRACE
|
||||||
break; // stop if bad conditioning
|
std::cerr << " BICGSTAB with preconditioner: omega=0" << std::endl;
|
||||||
beta=(alpha/omega)/rTh; rTh=BLAS<Vector>::dot(rT,h); beta*=rTh;
|
#endif
|
||||||
BLAS<Vector>::scal(beta,d);
|
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<Vector>::dot(rT,h);
|
||||||
|
beta*=rTh;
|
||||||
|
BLAS<Vector>::scal(beta,d); // d = h + beta * (d - omega * Sd);
|
||||||
BLAS<Vector>::axpy(1,h,d);
|
BLAS<Vector>::axpy(1,h,d);
|
||||||
BLAS<Vector>::axpy(-beta*omega,Sd,d);
|
BLAS<Vector>::axpy(-beta*omega,Sd,d);
|
||||||
rTr=BLAS<Vector>::dot(r,r);
|
|
||||||
++its;
|
++its;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue