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

If the user requests it, and CG seems to have converged, recompute the

Browse Source 
residual from scratch (as b-Ax) and use that to redo the convergence test.
This commit is contained in:
painter 2007-05-09 21:36:12 +00:00
parent f9adfb23bc
commit b52bf2037c
4 changed files with 120 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
krylov/HYPRE_pcg.c
View File

@ -219,6 +219,24 @@ HYPRE_PCGGetRelChange( HYPRE_Solver solver,
return( hypre_PCGGetRelChange( (void *) solver, rel_change ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetRecomputeResidual, HYPRE_PCGGetRecomputeResidual
*--------------------------------------------------------------------------*/
int
HYPRE_PCGSetRecomputeResidual( HYPRE_Solver solver,
int recompute_residual )
{
return( hypre_PCGSetRecomputeResidual( (void *) solver, recompute_residual ) );
}
int
HYPRE_PCGGetRecomputeResidual( HYPRE_Solver solver,
int * recompute_residual )
{
return( hypre_PCGGetRecomputeResidual( (void *) solver, recompute_residual ) );
}
/*--------------------------------------------------------------------------
* HYPRE_PCGSetPrecond
*--------------------------------------------------------------------------*/

11
krylov/krylov.h
View File

@ -709,7 +709,7 @@ typedef struct
where the norm is an energy norm wrt preconditioner, |r|=sqrt(<Cr,r>).
- two_norm!=0 means: the norm is the L2 norm, |r|=sqrt(<r,r>)
- rel_change!=0 means: if pass the other stopping criteria, also check the
relative change in the solution x.
relative change in the solution x. Pass iff this relative change is small.
- stop_crit!=0 means: pure absolute error tolerance rather than a pure relative
error tolerance on the residual. Never applies if rel_change!=0 or atolf!=0.
- atolf = absolute error tolerance factor to be used _together_ with the
@ -717,6 +717,10 @@ typedef struct
- tol = relative error tolerance, as above
- cf_tol = convergence factor tolerance; if >0 used for special test
for slow convergence
- recompute_residual means: when the iteration seems to be converged, recompute the
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.
*/
typedef struct
@ -727,6 +731,7 @@ typedef struct
int max_iter;
int two_norm;
int rel_change;
int recompute_residual;
int stop_crit;
int converged;
@ -951,6 +956,8 @@ int HYPRE_PCGSetTwoNorm ( HYPRE_Solver solver , int two_norm );
int HYPRE_PCGGetTwoNorm ( HYPRE_Solver solver , int *two_norm );
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_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 );
@ -981,6 +988,8 @@ int hypre_PCGSetTwoNorm ( void *pcg_vdata , int two_norm );
int hypre_PCGGetTwoNorm ( void *pcg_vdata , int *two_norm );
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_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 );

103
krylov/pcg.c
View File

@ -113,6 +113,7 @@ hypre_PCGCreate( hypre_PCGFunctions *pcg_functions )
(pcg_data -> max_iter) = 1000;
(pcg_data -> two_norm) = 0;
(pcg_data -> rel_change) = 0;
(pcg_data -> recompute_residual) = 0;
(pcg_data -> stop_crit) = 0;
(pcg_data -> converged) = 0;
(pcg_data -> owns_matvec_data ) = 1;
@ -288,6 +289,7 @@ hypre_PCGSolve( void *pcg_vdata,
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 stop_crit = (pcg_data -> stop_crit);
/*
int converged = (pcg_data -> converged);
@ -315,7 +317,8 @@ hypre_PCGSolve( void *pcg_vdata,
double weight;
double ratio;
double guard_zero_residual, sdotp;
double guard_zero_residual, sdotp;
int tentatively_converged = 0;
int i = 0;
int ierr = 0;
@ -447,7 +450,7 @@ hypre_PCGSolve( void *pcg_vdata,
if ( logging>0 || print_level>0 ) norms[0] = sqrt(i_prod_0);
}
if ( print_level > 1 && my_id==0 ) /* formerly for par_csr only */
if ( print_level > 1 && my_id==0 )
{
printf("\n\n");
if (two_norm)
@ -470,6 +473,9 @@ hypre_PCGSolve( void *pcg_vdata,
while ((i+1) <= max_iter)
{
/*--------------------------------------------------------------------
* the core CG calculations...
*--------------------------------------------------------------------*/
i++;
/* s = A*p */
@ -506,6 +512,9 @@ hypre_PCGSolve( void *pcg_vdata,
else
i_prod = gamma;
/*--------------------------------------------------------------------
* optional output
*--------------------------------------------------------------------*/
#if 0
if (two_norm)
printf("Iter (%d): ||r||_2 = %e, ||r||_2/||b||_2 = %e\n",
@ -543,25 +552,49 @@ hypre_PCGSolve( void *pcg_vdata,
}
/* check for convergence */
if (i_prod / bi_prod < eps)
/*--------------------------------------------------------------------
* check for convergence
*--------------------------------------------------------------------*/
if (i_prod / bi_prod < eps) /* the basic convergence test */
tentatively_converged = 1;
if ( tentatively_converged && recompute_residual )
/* At user request, don't trust the convergence test until we've recomputed
the residual from scratch. This is expensive in the usual case where an
the norm is the energy norm.
This calculation is coded on the assumption that r's accuracy is only a
concern for problems where CG takes many iterations. */
{
if (rel_change && i_prod > guard_zero_residual)
/* r = b - Ax */
(*(pcg_functions->CopyVector))(b, r);
(*(pcg_functions->Matvec))(matvec_data, -1.0, A, x, 1.0, r);
/* set i_prod for convergence test */
if (two_norm)
i_prod = (*(pcg_functions->InnerProd))(r,r);
else
{
/* s = C*r */
(*(pcg_functions->ClearVector))(s);
precond(precond_data, A, r, s);
/* iprod = gamma = <r,s> */
i_prod = (*(pcg_functions->InnerProd))(r, s);
}
if (i_prod / bi_prod >= eps) tentatively_converged = 0;
}
if ( tentatively_converged && rel_change && (i_prod > guard_zero_residual ))
/* At user request, don't treat this as converged unless x didn't change
much in the last iteration. */
{
pi_prod = (*(pcg_functions->InnerProd))(p,p);
xi_prod = (*(pcg_functions->InnerProd))(x,x);
ratio = alpha*alpha*pi_prod/xi_prod;
if (ratio < eps)
{
(pcg_data -> converged) = 1;
break;
}
}
else
{
(pcg_data -> converged) = 1;
break;
}
if (ratio >= eps) tentatively_converged = 0;
}
if ( tentatively_converged )
/* we've passed all the convergence tests, it's for real */
{
(pcg_data -> converged) = 1;
break;
}
if ( (gamma<1.0e-292) && ((-gamma)<1.0e-292) ) {
@ -605,6 +638,10 @@ hypre_PCGSolve( void *pcg_vdata,
if (weight * cf_ave_1 > cf_tol) break;
}
/*--------------------------------------------------------------------
* back to the core CG calculations
*--------------------------------------------------------------------*/
/* beta = gamma / gamma_old */
beta = gamma / gamma_old;
@ -613,7 +650,11 @@ hypre_PCGSolve( void *pcg_vdata,
(*(pcg_functions->Axpy))(1.0, s, p);
}
if ( print_level > 1 && my_id==0 ) /* formerly for par_csr only */
/*--------------------------------------------------------------------
* Finish up with some outputs.
*--------------------------------------------------------------------*/
if ( print_level > 1 && my_id==0 )
printf("\n\n");
(pcg_data -> num_iterations) = i;
@ -793,6 +834,34 @@ hypre_PCGGetRelChange( void *pcg_vdata,
return ierr;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetRecomputeResidual, hypre_PCGGetRecomputeResidual
*--------------------------------------------------------------------------*/
int
hypre_PCGSetRecomputeResidual( void *pcg_vdata,
int recompute_residual )
{
hypre_PCGData *pcg_data = pcg_vdata;
int ierr = 0;
(pcg_data -> recompute_residual) = recompute_residual;
return ierr;
}
int
hypre_PCGGetRecomputeResidual( void *pcg_vdata,
int * recompute_residual )
{
hypre_PCGData *pcg_data = pcg_vdata;
int ierr = 0;
*recompute_residual = (pcg_data -> recompute_residual);
return ierr;
}
/*--------------------------------------------------------------------------
* hypre_PCGSetStopCrit, hypre_PCGGetStopCrit
*--------------------------------------------------------------------------*/

7
krylov/pcg.h
View File

@ -100,7 +100,7 @@ typedef struct
where the norm is an energy norm wrt preconditioner, |r|=sqrt(<Cr,r>).
- two_norm!=0 means: the norm is the L2 norm, |r|=sqrt(<r,r>)
- rel_change!=0 means: if pass the other stopping criteria, also check the
relative change in the solution x.
relative change in the solution x. Pass iff this relative change is small.
- stop_crit!=0 means: pure absolute error tolerance rather than a pure relative
error tolerance on the residual. Never applies if rel_change!=0 or atolf!=0.
- atolf = absolute error tolerance factor to be used _together_ with the
@ -108,6 +108,10 @@ typedef struct
- tol = relative error tolerance, as above
- cf_tol = convergence factor tolerance; if >0 used for special test
for slow convergence
- recompute_residual means: when the iteration seems to be converged, recompute the
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.
*/
typedef struct
@ -118,6 +122,7 @@ typedef struct
int max_iter;
int two_norm;
int rel_change;
int recompute_residual;
int stop_crit;
int converged;