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

straightforward computation of residual sum-of-squares if more than one

Browse Source 
pointset (when the other method is inaccurate)
This commit is contained in:
painter 2006-12-12 22:15:57 +00:00
parent 446e6345d7
commit 057632667b
1 changed files with 47 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
struct_ls/point_relax.c
View File

@ -377,6 +377,7 @@ hypre_PointRelax( void *relax_vdata,
double *bp;
double *xp;
double *tp;
void *matvec_data = NULL;
int Ai;
int bi;
@ -435,6 +436,11 @@ hypre_PointRelax( void *relax_vdata,
p = 0;
iter = 0;
if ( tol>0.0 && num_pointsets>1 )
{
matvec_data = hypre_StructMatvecCreate();
hypre_StructMatvecSetup( matvec_data, A, x );
}
if (zero_guess)
{
@ -522,8 +528,8 @@ hypre_PointRelax( void *relax_vdata,
}
}
if ( tol>0.0 )
{
if ( tol>0.0 && num_pointsets==1 )
{ /* accumulate residual sum-squares for convergence test */
ierr += hypre_sumsqdiff_pointset( relax_data, pointset, x, b, A, &ss );
rsumsq += ss;
}
@ -542,6 +548,13 @@ hypre_PointRelax( void *relax_vdata,
tol>0.0 means to do a convergence test, using tol.
The test is simply ||r||/||b||<tol, where r=residual, b=r.h.s., unweighted L2 norm */
{
if ( num_pointsets>1 )
{ /* do residual sum-squares the expensive way */
hypre_StructCopy( b, t ); /* t = b */
hypre_StructMatvecCompute( matvec_data,
-1.0, A, x, 1.0, t ); /* t = - A x + t = - A x + b */
rsumsq = hypre_StructInnerProd( t, t ); /* <t,t> */
}
if ( rsumsq/bsumsq<tol2 ) max_iter = iter; /* converged; reset max_iter to prevent more iterations */
}
}
@ -646,8 +659,8 @@ hypre_PointRelax( void *relax_vdata,
}
if ( tol>0.0 )
{
if ( tol>0.0 && num_pointsets==1 )
{ /* accumulate residual sum-squares for convergence test */
ierr += hypre_sumsqdiff_pointset( relax_data, pointset, t, x, A, &ss );
rsumsq += ss;
}
@ -673,10 +686,22 @@ hypre_PointRelax( void *relax_vdata,
tol>0.0 means to do a convergence test, using tol.
The test is simply ||r||/||b||<tol, where r=residual, b=r.h.s., unweighted L2 norm */
{
if ( num_pointsets>1 )
{ /* do residual sum-squares the expensive way */
hypre_StructCopy( b, t ); /* t = b */
hypre_StructMatvecCompute( matvec_data,
-1.0, A, x, 1.0, t ); /* t = - A x + t = - A x + b */
rsumsq = hypre_StructInnerProd( t, t ); /* <t,t> */
}
if ( rsumsq/bsumsq<tol2 ) break;
}
}
if ( tol>0.0 && num_pointsets>1 )
{
hypre_StructMatvecDestroy( matvec_data );
}
if ( tol>0.0 ) (relax_data -> rresnorm) = sqrt( rsumsq/bsumsq );
(relax_data -> num_iterations) = iter;
@ -1876,7 +1901,17 @@ int hypre_relax_copy( void *relax_vdata, int pointset,
}
/*--------------------------------------------------------------------------
* hypre_sumsqdiff_pointset
* Computes the sum of squares of the residual but only on the specified pointset.
* Computes the sum of squares of the residual for a Jacobi iteration.
* This is a cheap computation, but the residual sum-of-squares actually belongs
* to the previous iteration. But note that a straightforward computation of
* a current value would cost as much as doing another iteration.
* Note also that PointRelax is doing a real Jacobi iteration only in the
* case where there is a single pointset. This is because if there are two or
* more pointsets, the iteration for the second pointset uses updated values
* from the first one. This function will not compute a true correct residual
* sum-of-squares in the multiple-pointset case, although what it computes may
* be useful in early stages of convergence testing.
*
* The residual is computed from the old x iterate, x, the new x iterate, t,
* and the diagonal of the matrix A as follows. We assume t was computed by
* a Jacobi iteration on the pointset.
@ -1950,14 +1985,14 @@ int hypre_sumsqdiff_pointset( void *relax_vdata, int pointset,
tp = hypre_StructVectorBoxData(t, i);
Ap = hypre_StructMatrixBoxData(A, i, diag_rank);
compute_box = hypre_BoxArrayBox(compute_box_a, j);
start = hypre_BoxIMin(compute_box);
hypre_BoxGetStrideSize(compute_box, stride, loop_size);
if ( constant_coefficient==1 ) /* constant diagonal */
{
hypre_ForBoxI(j, compute_box_a)
{
compute_box = hypre_BoxArrayBox(compute_box_a, j);
start = hypre_BoxIMin(compute_box);
hypre_BoxGetStrideSize(compute_box, stride, loop_size);
Ai = hypre_CCBoxIndexRank( A_data_box, start );
Api = Ap[Ai];
@ -1978,6 +2013,9 @@ int hypre_sumsqdiff_pointset( void *relax_vdata, int pointset,
{
hypre_ForBoxI(j, compute_box_a)
{
compute_box = hypre_BoxArrayBox(compute_box_a, j);
start = hypre_BoxIMin(compute_box);
hypre_BoxGetStrideSize(compute_box, stride, loop_size);
hypre_BoxLoop3Begin(loop_size,
A_data_box, start, stride, Ai,
x_data_box, start, stride, xi,