[CONTACT]

[ABOUT]

[POLICY]

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 ----------

.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]