# FORTRAN Corrected initial residual c

## Found at: ftp.icm.edu.pl:70/packages/netlib/templates/templates_code.errata

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"

---------- Jack Dongarra on April 22, 2017 ----------

.