Jan 2 94
FORTRAN code:
Corrected initial residual check in BiCG and QMR
( was R = A*X; i.e changed ZERO to -ONE )
Feb 11 94 Upon failure to converge, set INFO = ITER; was INFO = 1
Corrected in FORTRAN and C versions (single and double)
June 17, 1996
In the Fortran code for GMRESREVCOM ITER is not set to zero before
starting the iteration.
Before 10 continue in the code iter should be set to 0.
*
ITER = 0
10 CONTINUE
*
ITER = ITER + 1
----------
August 1, 1996: Corrected in Fortran versions of GMRES and GMRESREVCOM.
I. Problem I
The problem reported arises from inconsistency between allocation and
use of workspace for V.
In GMRES() routine, two working spaces, WORK(LDW,5) and
WORK2(LDW2,2*RESTRT+2) are initially assumed. Here, it is required
that LDW =< max(1,n), LDW2 =< max(1,RESTRT), and RESTRT =< n, where
n is the problem size.
In GMRESREVCOM(),
1. WORK is used as temporary storages for storing
5 intermediate vectors.
2. WORK2 is used for storing the upper Hessenberg matrix
H(m+1,m) and the vectors V(n,m) for restarted GMRES.
3. Extra two columns of WORK2, c(1:m) and s(1:m) are used
for storing Given rotations.
Here, m = RESTRT for restarted GMRES.
Thus, the initial layout of workspaces is assumed as below.
WORK(LDW,5) WORK2(LDW2,2*RESTRT+2)
_ _ _ _ _ _________________________
| | | | | | | | | | |
| | | | | | | | | | |
| | | | | | | | | | |
| | | | | | | | | | |
|R|S|W|Y|A| | H(m+1,m) | V(n,m) |c|s| m = RESTRT
| | | | |V| | | | | |
| | | | | | | | |---|
| | | | | | ----------| | |
| | | | | | | | | |
| | | | | | - - - - - |- - - - - |- - LDW2
| | | | | | | |
LDW - - - - - ----------
As shown, in their uses in GMRESREVCOM(), the vectors V have length n
which may not be as small as LDW2. In this case, incorrect access
of memory location would occur. Moreover, V consists of m+1 vectors
rather than m.
::::::::::
Solution
::::::::::
Allocate space for V in WORK instead of WORK2. The corrected layout
looks like the following.
WORK(LDW,6+RESTRT) WORK2(LDW2,RESTRT+2)
_ _ _ _ _ ____________ ______________
| | | | | | | | | | |
| | | | | | | | | | |
| | | | | | | | | | |
| | | | | | | | | | |
|R|S|W|Y|A| V(n,m+1) | | H(m+1,m) |c|s|
| | | | |V| | | | | |
| | | | | | | | |---|
| | | | | | | - - - - - - - LDW2 = m+1
| | | | | | |
| | | | | | |
| | | | | | |
- - - - - - - - - - - LDW = n
So, allocate the amount of WORK(LDW,6+RESTRT) and WORK2(LDW2,RESTRT+2).
This requires the users to allocate additional working space for V
in WORK in their routine. The routines QMRES() and QMRESREVCOM() need
to be modified (e.g. WORK2(1,V) --> WORK(1,V)).
::::::::::::::
Corrections
::::::::::::::
The following corrections were made in Fortran versions of GMRES and
GMRESREVCOM to work correctly. The corresponding corrections were
also made in both versions (double and single).
In this errata note, the corrections are explained only for the single
precision Fortran version.
1. GMRES.f
line 175: Now use the workspace of V allocated in array WORK.
CALL MATVEC(SCLR1, WORK2(NDX1), SCLR2, WORK(NDX2))
==> CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))
2. GMRESREVCOM.f
lines 152,153: correct the pointer to V in array WORK
AV = 5
V = 6
line 156: correct the pointer to GIV
GIV = H + RESTRT
line 172: correct the pointer
NEED1 = ((V - 1) * LDW) + 1
line 200: correct the pointer
NEED2 = ((V - 1) * LDW) + 1
line 241: added the following statement (was the previous errata)
ITER = 0
line 262,263: now use the workspace of V in array WORK
RNORM = SNRM2( N, WORK( 1,V ), 1 )
CALL SSCAL( N, RNORM, WORK( 1,V ), 1 )
line 267: now use the workspace of V in array WORK
******CALL MATVEC( ONE, WORK( 1,V+I-1 ), ZERO, WORK( 1,AV ) )
line 269,270: the correct leading dimension for V and AV is LDW.
NDX1 = ((V+I-1 - 1) * LDW) + 1
NDX2 = ((AV - 1) * LDW) + 1
line 298: now use the workspace of V in array WORK
CALL ORTHOH( I, N, WORK2( 1,I+H-1 ), WORK( 1,V ), LDW,
line 317: now use the workspace of V in array WORK
$ WORK(1,Y), WORK(1,S), WORK( 1,V ), LDW)
line 325: now use the workspace of V in array WORK
$ WORK(1,Y), WORK( 1,S ), WORK( 1,V ), LDW )
3. sspdchk.f
line 281: the workspace for V is allocated in array WORK now.
CALL GMRES( N, B, X(1,7), RESTRT, WORK, LDW,
$ WORK(6*LDW+1), LDW, ITER(7), RESID(7),
$ MATVEC, PSOLVE, INFO(7) )
==> CALL GMRES( N, B, X(1,7), RESTRT, WORK, LDW,
$ WORK((6+RESTRT)*LDW+1), LDW, ITER(7),
$ RESID(7), MATVEC, PSOLVE, INFO(7) )
4. snsychk.f
line 242: the workspace for V is allocated in array WORK now.
CALL GMRES( N, B, X(1,7), RESTRT, WORK, LDW,
$ WORK(6*LDW+1), LDW, ITER(7), RESID(7),
$ MATVEC, PSOLVE, INFO(7))
==> CALL GMRES( N, B, X(1,7), RESTRT, WORK, LDW,
$ WORK((RESTRT+6)*LDW+1), LDW, ITER(7),
$ RESID(7), MATVEC, PSOLVE, INFO(7) )
II. Other Problems
1. GMRESREVCOM.f
Line 313: Pass the correct pointer to WORK(1,S) (was WORK(I,S))
RESID = APPROXRES( I, WORK2( 1,I+H-1 ), WORK( I,S ),
==>
RESID = APPROXRES( I, WORK2( 1,I+H-1 ), WORK( 1,S ),
2. GMRESREVCOM.f
line 475: Need the absolute value for residual.
APPROXRES = S( I+1 )
==>
APPROXRES = ABS( S( I+1 ) )
line 466: declare the intrinsic function of ABS()
INTRINSIC ABS
III. The above problems don't occur in C versions.
---------- Done by Youngbae Kim on August 1, 1996 ----------
----------
March 7, 2006
Conjugate Gradient Method (2.11) correction:
p(i)=r(i)+beta(i-1)*p(i-1)
==> p(i)=r(i-1)+beta(i-1)*p(i-1)
URL location: node20.html#SECTION00731000000000000000
Book location: page 14
----------