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 ---------- ---------- April 22, 2017 On p.25, right after Eq. 2.15, it says "where $a$ is the length of the $x$-axis of the ellipse" ==> "where $2 a$ is the length of the $x$-axis of the ellipse" reported by Greg Hammett ---------- Jack Dongarra on April 22, 2017 ----------