CFD/hypre
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added two new features to PCG:

Browse Source 
1) A new, residual-based stopping criteria ||r_old - r_new|| < rtol ||b|| which
   is activated with HYPRE_PCGSetResidualTol() and is designed to be used in
   high-accuracy computations with rtol <= 1e-16.

2) HYPRE_PCGSetRecomputeResidual(), which periodically recomputes the residual
   from scratch (r=b-Ax). This, again, can be useful in high-accuracy
   computations, but is potentially expensive and can lead to degraded
   convergence!
This commit is contained in:
kolev1 2009-07-20 00:16:39 +00:00
parent e1b1b3170b
commit 4502255f79
5 changed files with 176 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
krylov/HYPRE_krylov.h
View File

@ -122,8 +122,15 @@ int HYPRE_PCGSetTol(HYPRE_Solver solver,
int HYPRE_PCGSetAbsoluteTol(HYPRE_Solver solver,
double a_tol);
/*
/**
* (Optional) Set a residual-based convergence tolerance which checks if
* $\|r_old-r_new\| < rtol \|b\|$. This is useful when trying to converge to
* very low relative and/or absolute tolerances, in order to bail-out before
* roundoff errors affect the approximation.
**/
int HYPRE_PCGSetResidualTol(HYPRE_Solver solver,
double rtol);
/*
* RE-VISIT
**/
int HYPRE_PCGSetAbsoluteTolFactor(HYPRE_Solver solver, double abstolf);
@ -157,6 +164,18 @@ int HYPRE_PCGSetTwoNorm(HYPRE_Solver solver,
int HYPRE_PCGSetRelChange(HYPRE_Solver solver,
int rel_change);
/**
* (Optional) Recompute the residual at the end to double-check convergence.
**/
int HYPRE_PCGSetRecomputeResidual( HYPRE_Solver solver,
int recompute_residual );
/**
* (Optional) Periodically recompute the residual while iterating.
**/
int HYPRE_PCGSetRecomputeResidualP( HYPRE_Solver solver,
int recompute_residual_p );
/**
* (Optional) Set the preconditioner to use.
**/
@ -199,6 +218,10 @@ int HYPRE_PCGGetResidual(HYPRE_Solver solver,
**/
int HYPRE_PCGGetTol(HYPRE_Solver solver, double *tol);
/**
**/
int HYPRE_PCGGetResidualTol(HYPRE_Solver solver, double *rtol);
/*
* RE-VISIT
**/

36
krylov/HYPRE_pcg.c
View File

@ -113,6 +113,24 @@ HYPRE_PCGGetAbsoluteTol( HYPRE_Solver solver,
return( hypre_PCGGetAbsoluteTol( (void *) solver, a_tol ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetResidualTol, HYPRE_PCGGetResidualTol
*--------------------------------------------------------------------------*/
int
HYPRE_PCGSetResidualTol( HYPRE_Solver solver,
double rtol )
{
return( hypre_PCGSetResidualTol( (void *) solver, rtol ) );
}
int
HYPRE_PCGGetResidualTol( HYPRE_Solver solver,
double * rtol )
{
return( hypre_PCGGetResidualTol( (void *) solver, rtol ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetAbsoluteTolFactor, HYPRE_PCGGetAbsoluteTolFactor
*--------------------------------------------------------------------------*/
@ -241,6 +259,24 @@ HYPRE_PCGGetRecomputeResidual( HYPRE_Solver solver,
return( hypre_PCGGetRecomputeResidual( (void *) solver, recompute_residual ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetRecomputeResidualP, HYPRE_PCGGetRecomputeResidualP
*--------------------------------------------------------------------------*/
int
HYPRE_PCGSetRecomputeResidualP( HYPRE_Solver solver,
int recompute_residual_p )
{
return( hypre_PCGSetRecomputeResidualP( (void *) solver, recompute_residual_p ) );
}
int
HYPRE_PCGGetRecomputeResidualP( HYPRE_Solver solver,
int * recompute_residual_p )
{
return( hypre_PCGGetRecomputeResidualP( (void *) solver, recompute_residual_p ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetPrecond
*--------------------------------------------------------------------------*/

16
krylov/krylov.h
View File

@ -1010,7 +1010,11 @@ typedef struct
residual from scratch (r=b-Ax) and use this new residual to repeat the convergence test.
This can be expensive, use this only if you have seen a problem with the regular
residual computation.
*/
- recompute_residual_p means: recompute the residual from scratch (r=b-Ax)
every "recompute_residual_p" iterations. This can be expensive and degrade the
convergence. Use it only if you have seen a problem with the regular residual
computation.
*/
typedef struct
{
@ -1018,10 +1022,12 @@ typedef struct
double atolf;
double cf_tol;
double a_tol;
double rtol;
int max_iter;
int two_norm;
int rel_change;
int recompute_residual;
int recompute_residual_p;
int stop_crit;
int converged;
@ -1365,6 +1371,8 @@ int HYPRE_PCGSetAbsoluteTol ( HYPRE_Solver solver , double a_tol );
int HYPRE_PCGGetAbsoluteTol ( HYPRE_Solver solver , double *a_tol );
int HYPRE_PCGSetAbsoluteTolFactor ( HYPRE_Solver solver , double abstolf );
int HYPRE_PCGGetAbsoluteTolFactor ( HYPRE_Solver solver , double *abstolf );
int HYPRE_PCGSetResidualTol ( HYPRE_Solver solver , double rtol );
int HYPRE_PCGGetResidualTol ( HYPRE_Solver solver , double *rtol );
int HYPRE_PCGSetConvergenceFactorTol ( HYPRE_Solver solver , double cf_tol );
int HYPRE_PCGGetConvergenceFactorTol ( HYPRE_Solver solver , double *cf_tol );
int HYPRE_PCGSetMaxIter ( HYPRE_Solver solver , int max_iter );
@ -1377,6 +1385,8 @@ int HYPRE_PCGSetRelChange ( HYPRE_Solver solver , int rel_change );
int HYPRE_PCGGetRelChange ( HYPRE_Solver solver , int *rel_change );
int HYPRE_PCGSetRecomputeResidual ( HYPRE_Solver solver , int recompute_residual );
int HYPRE_PCGGetRecomputeResidual ( HYPRE_Solver solver , int *recompute_residual );
int HYPRE_PCGSetRecomputeResidualP ( HYPRE_Solver solver , int recompute_residual_p );
int HYPRE_PCGGetRecomputeResidualP ( HYPRE_Solver solver , int *recompute_residual_p );
int HYPRE_PCGSetPrecond ( HYPRE_Solver solver , HYPRE_PtrToSolverFcn precond , HYPRE_PtrToSolverFcn precond_setup , HYPRE_Solver precond_solver );
int HYPRE_PCGGetPrecond ( HYPRE_Solver solver , HYPRE_Solver *precond_data_ptr );
int HYPRE_PCGSetLogging ( HYPRE_Solver solver , int level );
@ -1401,6 +1411,8 @@ int hypre_PCGSetAbsoluteTol ( void *pcg_vdata , double a_tol );
int hypre_PCGGetAbsoluteTol ( void *pcg_vdata , double *a_tol );
int hypre_PCGSetAbsoluteTolFactor ( void *pcg_vdata , double atolf );
int hypre_PCGGetAbsoluteTolFactor ( void *pcg_vdata , double *atolf );
int hypre_PCGSetResidualTol ( void *pcg_vdata , double rtol );
int hypre_PCGGetResidualTol ( void *pcg_vdata , double *rtol );
int hypre_PCGSetConvergenceFactorTol ( void *pcg_vdata , double cf_tol );
int hypre_PCGGetConvergenceFactorTol ( void *pcg_vdata , double *cf_tol );
int hypre_PCGSetMaxIter ( void *pcg_vdata , int max_iter );
@ -1411,6 +1423,8 @@ int hypre_PCGSetRelChange ( void *pcg_vdata , int rel_change );
int hypre_PCGGetRelChange ( void *pcg_vdata , int *rel_change );
int hypre_PCGSetRecomputeResidual ( void *pcg_vdata , int recompute_residual );
int hypre_PCGGetRecomputeResidual ( void *pcg_vdata , int *recompute_residual );
int hypre_PCGSetRecomputeResidualP ( void *pcg_vdata , int recompute_residual_p );
int hypre_PCGGetRecomputeResidualP ( void *pcg_vdata , int *recompute_residual_p );
int hypre_PCGSetStopCrit ( void *pcg_vdata , int stop_crit );
int hypre_PCGGetStopCrit ( void *pcg_vdata , int *stop_crit );
int hypre_PCGGetPrecond ( void *pcg_vdata , HYPRE_Solver *precond_data_ptr );

96
krylov/pcg.c
View File

@ -98,10 +98,12 @@ hypre_PCGCreate( hypre_PCGFunctions *pcg_functions )
(pcg_data -> atolf) = 0.0;
(pcg_data -> cf_tol) = 0.0;
(pcg_data -> a_tol) = 0.0;
(pcg_data -> rtol) = 0.0;
(pcg_data -> max_iter) = 1000;
(pcg_data -> two_norm) = 0;
(pcg_data -> rel_change) = 0;
(pcg_data -> recompute_residual) = 0;
(pcg_data -> recompute_residual_p) = 0;
(pcg_data -> stop_crit) = 0;
(pcg_data -> converged) = 0;
(pcg_data -> owns_matvec_data ) = 1;
@ -274,10 +276,12 @@ hypre_PCGSolve( void *pcg_vdata,
double a_tol = (pcg_data -> a_tol);
double atolf = (pcg_data -> atolf);
double cf_tol = (pcg_data -> cf_tol);
double rtol = (pcg_data -> rtol);
int max_iter = (pcg_data -> max_iter);
int two_norm = (pcg_data -> two_norm);
int rel_change = (pcg_data -> rel_change);
int recompute_residual = (pcg_data -> recompute_residual);
int recompute_residual_p = (pcg_data -> recompute_residual_p);
int stop_crit = (pcg_data -> stop_crit);
/*
int converged = (pcg_data -> converged);
@ -494,8 +498,32 @@ hypre_PCGSolve( void *pcg_vdata,
(*(pcg_functions->Axpy))(alpha, p, x);
/* r = r - alpha*s */
(*(pcg_functions->Axpy))(-alpha, s, r);
if ( !recompute_residual_p || ((i+1)%(recompute_residual_p)) )
{
(*(pcg_functions->Axpy))(-alpha, s, r);
}
else
{
/* At user request, periodically recompute the residual from the formula
r = b - A x (instead of the above recursive definition). Note that this
is potentially expensive and can lead to degraded convergence! */
printf("Recomputing the residual...\n");
(*(pcg_functions->CopyVector))(b, r);
(*(pcg_functions->Matvec))(matvec_data, -1.0, A, x, 1.0, r);
}
/* residual-based stopping criteria: ||r_new-r_old|| < rtol ||b|| */
if (rtol && two_norm)
{
/* use that r_new-r_old = alpha * s */
double drob2 = alpha*alpha*(*(pcg_functions->InnerProd))(s,s)/bi_prod;
if ( drob2 < rtol*rtol )
{
printf("\n\n||r_old-r_new||/||b||: %e\n", sqrt(drob2));
break;
}
}
/* s = C*r */
(*(pcg_functions->ClearVector))(s);
precond(precond_data, A, r, s);
@ -503,6 +531,18 @@ hypre_PCGSolve( void *pcg_vdata,
/* gamma = <r,s> */
gamma = (*(pcg_functions->InnerProd))(r, s);
/* residual-based stopping criteria: ||r_new-r_old||_C < rtol ||b||_C */
if (rtol && !two_norm)
{
/* use that ||r_new-r_old||_C^2 = (r_new ,C r_new) + (r_old, C r_old) */
double r2ob2 = (gamma + gamma_old)/bi_prod;
if ( r2ob2 < rtol*rtol)
{
printf("\n\n||r_old-r_new||_C/||b||_C: %e\n", sqrt(r2ob2));
break;
}
}
/* set i_prod for convergence test */
if (two_norm)
i_prod = (*(pcg_functions->InnerProd))(r,r);
@ -744,6 +784,32 @@ hypre_PCGGetAbsoluteTolFactor( void *pcg_vdata,
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetResidualTol, hypre_PCGGetResidualTol
*--------------------------------------------------------------------------*/
int
hypre_PCGSetResidualTol( void *pcg_vdata,
double rtol )
{
hypre_PCGData *pcg_data = pcg_vdata;
(pcg_data -> rtol) = rtol;
return hypre_error_flag;
}
int
hypre_PCGGetResidualTol( void *pcg_vdata,
double * rtol )
{
hypre_PCGData *pcg_data = pcg_vdata;
*rtol = (pcg_data -> rtol);
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetConvergenceFactorTol, hypre_PCGGetConvergenceFactorTol
*--------------------------------------------------------------------------*/
@ -881,6 +947,32 @@ hypre_PCGGetRecomputeResidual( void *pcg_vdata,
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetRecomputeResidualP, hypre_PCGGetRecomputeResidualP
*--------------------------------------------------------------------------*/
int
hypre_PCGSetRecomputeResidualP( void *pcg_vdata,
int recompute_residual_p )
{
hypre_PCGData *pcg_data = pcg_vdata;
(pcg_data -> recompute_residual_p) = recompute_residual_p;
return hypre_error_flag;
}
int
hypre_PCGGetRecomputeResidualP( void *pcg_vdata,
int * recompute_residual_p )
{
hypre_PCGData *pcg_data = pcg_vdata;
*recompute_residual_p = (pcg_data -> recompute_residual_p);
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetStopCrit, hypre_PCGGetStopCrit
*--------------------------------------------------------------------------*/

6
krylov/pcg.h
View File

@ -106,6 +106,10 @@ typedef struct
residual from scratch (r=b-Ax) and use this new residual to repeat the convergence test.
This can be expensive, use this only if you have seen a problem with the regular
residual computation.
- recompute_residual_p means: recompute the residual from scratch (r=b-Ax)
every "recompute_residual_p" iterations. This can be expensive and degrade the
convergence. Use it only if you have seen a problem with the regular residual
computation.
*/
typedef struct
@ -114,10 +118,12 @@ typedef struct
double atolf;
double cf_tol;
double a_tol;
double rtol;
int max_iter;
int two_norm;
int rel_change;
int recompute_residual;
int recompute_residual_p;
int stop_crit;
int converged;