7.1.1 

Inaccuracy.f90:

! Don't worry about this - it is described in the next lecture.  All it
! does is declare an interface for the functions used here.

MODULE declarations
    INTERFACE

    FUNCTION copy (arg)
        IMPLICIT NONE
        REAL :: copy, arg
    END FUNCTION copy

    FUNCTION series ()
        IMPLICIT NONE
        REAL :: series
    END FUNCTION series

    FUNCTION decode (arg)
        IMPLICIT NONE
        REAL :: decode
        CHARACTER(LEN=*) :: arg
    END FUNCTION decode

    END INTERFACE
END MODULE declarations



FUNCTION copy (arg)
    IMPLICIT NONE
    REAL :: copy, arg
    copy = arg
END FUNCTION copy

FUNCTION decode (arg)
!
! This is not truly portable, but you would be unlucky to have problems
! nowadays, and ISO C actually requires '0' to '9' to be contiguous in
! numeric value.  Note that it is not the best way to code this,
! numerically, but is intended to show the issues.
!
    IMPLICIT NONE
    REAL :: decode, t, s
    CHARACTER(LEN=*) :: arg
    INTEGER :: i
    CHARACTER(LEN=*), PARAMETER :: format = &
        "('Invalid format of number: ''',A,'''')"
    t = 0.0
    s = 0.0
    DO i = 1,LEN(arg)
        IF (arg(i:i) .EQ. ' ') THEN
            CONTINUE
        ELSE IF (arg(i:i) .EQ. '.' .AND. s .EQ. 0.0) THEN
            s = 1.0
        ELSE IF (arg(i:i) .GE. '0' .AND. arg(i:i) .LE. '9') THEN
            t = t*10.0+(ICHAR(arg(i:i))-ICHAR('0'))
            s = s*0.1
        ELSE
            WRITE (*, format) arg
            STOP
        END IF
    END DO
    decode = t*s
END FUNCTION decode

FUNCTION series ()
!
! This is Leibnitz's series for pi, with pairs of terms collected
! together and simplified.  Please note that it is not a particularly
! poorly (or well) behaved series, and is NOT being summed very cleverly,
! but this is done deliberately to demonstrate the effect of limited
! precision.
!
    IMPLICIT NONE
    INTEGER :: n
    REAL :: series, t, y
    CHARACTER(LEN=*), PARAMETER :: format = &
        "('The series converged after ',I9,' iterations')"
    n = 0
    t = 0.0
    DO
       y = t
       t = t+2.0/((4.0*n+1.0)*(4.0*n+3.0))
       n = n+1
       IF (y .EQ. t) EXIT
    END DO
    WRITE (*, format) n
    series = 4.0*t
END FUNCTION series

PROGRAM verysimple
    USE declarations
    IMPLICIT NONE
    REAL, PARAMETER :: a  = 3.14159265358979323846, &
        b  = 314159.265358979323846, p = 1.0
    REAL :: c, d, e, f, g
    CHARACTER(LEN=100) :: buf1, buf2
    c = 3.14159265358979323846
    d = 314159.265358979323846
    OPEN (UNIT=1, FILE='inaccuracy.data', ACTION='READ')
    READ (1, '(A)') buf1
    WRITE (*, '(3X,A)') buf1(1:25)
    READ (buf1, *) e
    READ (1, '(A)') buf2
    WRITE (*, '(3X,A)') buf2(1:25)
    READ (buf2, *) f
    WRITE (*, '(1X)')
!
! Display pi and 10^5*pi calculated directly in several ways, and note
! the discrepancies.
!
    WRITE (*, '(F25.20)') a
    WRITE (*, '(F25.15)') b
    WRITE (*, '(F25.15)') 1.0e5*a
    WRITE (*, '(F25.20)') c
    WRITE (*, '(F25.15)') d
    WRITE (*, '(F25.15)') 1.0e5*c
    WRITE (*, '(F25.20)') e
    WRITE (*, '(F25.15)') f
    WRITE (*, '(F25.15)') 1.0e5*e
    WRITE (*, '(F25.20)') 4.0*atan(1.0)
    WRITE (*, '(F25.20)') 4.0*atan(p)
    WRITE (*, '(F25.15)') 4.0e5*atan(1.0)
!
! For arcane historical reasons, Fortran does not allow I/O in functions
! called from I/O lists; in THIS case, it wouldn't be a good idea even
! if in compilers that do allow it, as both are writing to the same unit.
!
    g = copy(3.14159265358979323846)
    WRITE (*, '(F25.20)') g
    g = copy(314159.265358979323846)
    WRITE (*, '(F25.15)') g
    g = copy(a)
    WRITE (*, '(F25.20)') g
    g = copy(c)
    WRITE (*, '(F25.20)') g
    WRITE (*, *)
!
! Note the loss of accuracy, especially in the last one.
!
    g = decode(buf1)
    WRITE (*, '(F25.20)') g
    g = decode(buf2)
    WRITE (*, '(F25.15)') g
    g = series()
    WRITE (*, '(F25.20)') g
END PROGRAM verysimple

inaccuracy.data:
3.14159265358979323846
314159.265358979323846

output:
   3.14159265358979323846   
   314159.265358979323846   

   3.14159274101257324219
   314159.250000000000000
   314159.281250000000000
   3.14159274101257324219
   314159.250000000000000
   314159.281250000000000
   3.14159274101257324219
   314159.250000000000000
   314159.281250000000000
   3.14159274101257324219
   3.14159274101257324219
   314159.281250000000000
   3.14159274101257324219
   314159.250000000000000
   3.14159274101257324219
   3.14159274101257324219

   3.14159464836120605469
   314159.375000000000000
The series converged after      2049 iterations
   3.14138388633728027344

7.1.2 modified version of inaccuracy.f90 to include double.f90:

! Don't worry about this - it is described in the next lecture.  All it
! does is declare an interface for the functions used here.

MODULE declarations
    INTERFACE

    FUNCTION copy (arg)
        USE double
        IMPLICIT NONE
        REAL(KIND=dp) :: copy, arg
    END FUNCTION copy

    FUNCTION series ()
        USE double
        IMPLICIT NONE
        REAL(KIND=dp) :: series
    END FUNCTION series

    FUNCTION decode (arg)
        USE double
        IMPLICIT NONE
        REAL(KIND=dp) :: decode
        CHARACTER(LEN=*) :: arg
    END FUNCTION decode

    END INTERFACE
END MODULE declarations



FUNCTION copy (arg)
    USE double
    IMPLICIT NONE
    REAL(KIND=dp) :: copy, arg
    copy = arg
END FUNCTION copy

FUNCTION decode (arg)
!
! This is not truly portable, but you would be unlucky to have problems
! nowadays, and ISO C actually requires '0' to '9' to be contiguous in
! numeric value.  Note that it is not the best way to code this,
! numerically, but is intended to show the issues.
!
    USE double
    IMPLICIT NONE
    REAL(KIND=dp) :: decode, t, s
    CHARACTER(LEN=*) :: arg
    INTEGER :: i
    CHARACTER(LEN=*), PARAMETER :: format = &
        "('Invalid format of number: ''',A,'''')"
    t = 0.0_dp
    s = 0.0_dp
    DO i = 1,LEN(arg)
        IF (arg(i:i) .EQ. ' ') THEN
            CONTINUE
        ELSE IF (arg(i:i) .EQ. '.' .AND. s .EQ. 0.0_dp) THEN
            s = 1.0_dp
        ELSE IF (arg(i:i) .GE. '0' .AND. arg(i:i) .LE. '9') THEN
            t = t*10.0_dp+(ICHAR(arg(i:i))-ICHAR('0'))
            s = s*0.1_dp
        ELSE
            WRITE (*, format) arg
            STOP
        END IF
    END DO
    decode = t*s
END FUNCTION decode

FUNCTION series ()
!
! This is Leibnitz's series for pi, with pairs of terms collected
! together and simplified.  Please note that it is not a particularly
! poorly (or well) behaved series, and is NOT being summed very cleverly,
! but this is done deliberately to demonstrate the effect of limited
! precision.
!
    USE double
    IMPLICIT NONE
    INTEGER :: n
    REAL(KIND=dp) :: series, t, y
    CHARACTER(LEN=*), PARAMETER :: format = &
        "('The series converged after ',I9,' iterations')"
    n = 0
    t = 0.0_dp
    DO
       y = t
       t = t+2.0_dp/((4.0_dp*n+1.0_dp)*(4.0_dp*n+3.0_dp))
       n = n+1
       IF (y .EQ. t) EXIT
    END DO
    WRITE (*, format) n
    series = 4.0_dp*t
END FUNCTION series

PROGRAM verysimple
    USE declarations
    USE double
    IMPLICIT NONE
    REAL(KIND=dp), PARAMETER :: a  = 3.14159265358979323846_dp, &
        b  = 314159.265358979323846_dp, p = 1.0_dp
    REAL(KIND=dp) :: c, d, e, f, g
    CHARACTER(LEN=100) :: buf1, buf2
    c = 3.14159265358979323846_dp
    d = 314159.265358979323846_dp
    OPEN (UNIT=1, FILE='inaccuracy.data', ACTION='READ')
    READ (1, '(A)') buf1
    WRITE (*, '(3X,A)') buf1(1:25)
    READ (buf1, *) e
    READ (1, '(A)') buf2
    WRITE (*, '(3X,A)') buf2(1:25)
    READ (buf2, *) f
    WRITE (*, '(1X)')
!
! Display pi and 10^5*pi calculated directly in several ways, and note
! the discrepancies.
!
    WRITE (*, '(F25.20)') a
    WRITE (*, '(F25.15)') b
    WRITE (*, '(F25.15)') 1.0e5_dp*a
    WRITE (*, '(F25.20)') c
    WRITE (*, '(F25.15)') d
    WRITE (*, '(F25.15)') 1.0e5_dp*c
    WRITE (*, '(F25.20)') e
    WRITE (*, '(F25.15)') f
    WRITE (*, '(F25.15)') 1.0e5_dp*e
    WRITE (*, '(F25.20)') 4.0_dp*atan(1.0_dp)
    WRITE (*, '(F25.20)') 4.0_dp*atan(p)
    WRITE (*, '(F25.15)') 4.0e5_dp*atan(1.0_dp)
!
! For arcane historical reasons, Fortran does not allow I/O in functions
! called from I/O lists; in THIS case, it wouldn't be a good idea even
! if in compilers that do allow it, as both are writing to the same unit.
!
    g = copy(3.14159265358979323846_dp)
    WRITE (*, '(F25.20)') g
    g = copy(314159.265358979323846_dp)
    WRITE (*, '(F25.15)') g
    g = copy(a)
    WRITE (*, '(F25.20)') g
    g = copy(c)
    WRITE (*, '(F25.20)') g
    WRITE (*, *)
!
! Note the loss of accuracy, especially in the last one.
!
    g = decode(buf1)
    WRITE (*, '(F25.20)') g
    g = decode(buf2)
    WRITE (*, '(F25.15)') g
    g = series()
    WRITE (*, '(F25.20)') g
END PROGRAM verysimple

output:

   3.14159265358979323846   
   314159.265358979323846   
 
   3.14159265358979311600
   314159.265358979348093
   314159.265358979289886
   3.14159265358979311600
   314159.265358979348093
   314159.265358979289886
   3.14159265358979311600
   314159.265358979348093
   314159.265358979289886
   3.14159265358979311600
   3.14159265358979311600
   314159.265358979289886
   3.14159265358979311600
   314159.265358979348093
   3.14159265358979311600
   3.14159265358979311600

   3.14159265358979711280
   314159.265358979580924
The series converged after  47441554 iterations
   3.14159264457183695640

7.2.1 cuberoot.f90:

PROGRAM Newton
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: target, current
    REAL(KIND=DP), DIMENSION(5) :: previous
    INTEGER :: i, k
mainloop: DO
        PRINT *, 'Type in a real number'
        READ (*, *, IOSTAT=k) target
        IF (k < 0) THEN
            STOP
        ELSE IF (k > 0) THEN
            PRINT *, 'Some sort of horrible I/O error'
            STOP
        END IF 
        IF (target .EQ. 0.0_DP) THEN
            PRINT *, 'The cube root of zero is, er, zero'
            CYCLE
        END IF
!
! This is cheating, but the hardest part of Newton-Raphson solution of
! Nth roots is getting the starting value, and doing so properly would
! merely be confusing.  So use a horrible hack to get an approximation.
!
        current = 1.1*CMPLX(target,KIND=KIND(0.0))**0.3
        DO i = 1,5
            previous(i) = 0.0_DP
        END DO
loop:   DO
            current = current - &
                value(current,target)/derivative(current)
            PRINT *, current
            DO i = 1,5
                if (current .EQ. previous(i)) EXIT loop
            END DO
            DO i = 1,4
                previous(i+1) = previous(i)
            END DO
            previous(1) = current
        END DO loop
        PRINT *, current**3, current**3.0_DP
    END DO mainloop

CONTAINS

FUNCTION value (arg, targ)
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: value, arg, targ
    value = arg**3-targ
END FUNCTION value

FUNCTION derivative (arg)
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: derivative, arg
    derivative = 3*arg**2
END FUNCTION derivative

END PROGRAM Newton

output:
 Type in a real number
1.23
   1.0795850679102725
   1.0715025475716751
   1.0714412731951111
   1.0714412696907731
   1.0714412696907731
   1.2300000000000000   1.2300000000000000
 Type in a real number
4.56
   1.6615392639846105
   1.6582753012424118
   1.6582688683948941
   1.6582688683699394
   1.6582688683699394
   4.5599999999999996   4.5599999999999996
 Type in a real number
-3
  -0.6380672516545566
  -2.8815971836389327
  -2.0414944484631032
  -1.6009368177339121
  -1.4574591827210368
  -1.4424077411895482
  -1.4422495876514039
  -1.4422495703074085
  -1.4422495703074083
  -1.4422495703074083
  -3.0000000000000000  -3.0000000000000000

7.2.2 modified cuberoot.f90:

PROGRAM Newton
    USE double
    IMPLICIT NONE
    COMPLEX(KIND=DP) :: target, current
    COMPLEX(KIND=DP), DIMENSION(5) :: previous
    INTEGER :: i, k
mainloop: DO
        PRINT *, 'Type in a complex number'
        READ (*, *, IOSTAT=k) target
        IF (k < 0) THEN
            STOP
        ELSE IF (k > 0) THEN
            PRINT *, 'Some sort of horrible I/O error'
            STOP
        END IF 
        IF (target .EQ. 0.0_DP) THEN
            PRINT *, 'The cube root of zero is, er, zero'
            CYCLE
        END IF
!
! This is cheating, but the hardest part of Newton-Raphson solution of
! Nth roots is getting the starting value, and doing so properly would
! merely be confusing.  So use a horrible hack to get an approximation.
!
        current = 1.1*CMPLX(target,KIND=KIND(0.0))**0.3
        DO i = 1,5
            previous(i) = 0.0_DP
        END DO
loop:   DO
            current = current - &
                value(current,target)/derivative(current)
            PRINT *, current
            DO i = 1,5
                if (current .EQ. previous(i)) EXIT loop
            END DO
            DO i = 1,4
                previous(i+1) = previous(i)
            END DO
            previous(1) = current
        END DO loop
        PRINT *, current**3, current**3.0_DP
    END DO mainloop

CONTAINS

FUNCTION value (arg, targ)
    USE double
    IMPLICIT NONE
    COMPLEX(KIND=DP) :: value, arg, targ
    value = arg**3-targ
END FUNCTION value

FUNCTION derivative (arg)
    USE double
    IMPLICIT NONE
    COMPLEX(KIND=DP) :: derivative, arg
    derivative = 3*arg**2
END FUNCTION derivative

END PROGRAM Newton

output:
 Type in a complex number
(1.23,0.0)
 (1.0795850679102725,0.0000000000000000)
 (1.0715025475716748,0.0000000000000000)
 (1.0714412731951111,0.0000000000000000)
 (1.0714412696907731,0.0000000000000000)
 (1.0714412696907731,0.0000000000000000)
 (1.2300000000000000,0.0000000000000000) (1.2300000000000000,0.0000000000000000)
 Type in a complex number
(4.56,0.0)
 (1.6615392639846105,0.0000000000000000)
 (1.6582753012424118,0.0000000000000000)
 (1.6582688683948941,0.0000000000000000)
 (1.6582688683699394,0.0000000000000000)
 (1.6582688683699394,0.0000000000000000)
 (4.5599999999999996,0.0000000000000000) (4.5599999999999996,0.0000000000000000)
 Type in a complex number
(-3,0.0)
 (0.7314233808603878,1.2314709447793017)
 (0.7208350727638332,1.2490159447439546)
 (0.7211248172875359,1.2490247178887515)
 (0.7211247851537018,1.2490247664834062)
 (0.7211247851537042,1.2490247664834064)
 (0.7211247851537042,1.2490247664834064)
 (-3.0000000000000000,0.0000000000000000) (-3.0000000000000000,0.0000000000000000)
 Type in a complex number
(1.23,4.56)
 (1.5237177133208530,0.7025657833679676)
 (1.5209703771166918,0.7082105561740765)
 (1.5209913463373150,0.7082211456217790)
 (1.5209913460485927,0.7082211454642361)
 (1.5209913460485927,0.7082211454642361)
 (1.2299999999999998,4.5599999999999996) (1.2299999999999998,4.5599999999999996)

7.3.1

The argument to function CPY is of the wrong type, and this is diagnosed.

NAG Fortran Compiler Release 5.2(668)
Error: ex_07.1.2.f90, line 135: Incorrect data type REAL (expected DOUBLE PRECISION) for argument ARG (no. 1) of COPY
Error: ex_07.1.2.f90, line 137: Incorrect data type REAL (expected DOUBLE PRECISION) for argument ARG (no. 1) of COPY
[NAG Fortran Compiler error termination, 2 errors]

7.3.2

The program runs, but the results are wrong, because the constants
are stored to only 6 digits accuracy.