errors_prb.f90

数值计算常用的出错处理！可以看看！学习一下：）

program errors_prb
!
!*******************************************************************************
!
!! ERRORS_PRB runs the ERRORS examples.
!
!
!  Discussion:
!
!    This program demonstrates how unreliable certain
!    calculations can be when performed with a finite accuracy
!    on a computer.  Standard software is used, where possible,
!    for the solution of these problems.
!
!  Quoting from the reference:
!
!    Despite the increasing spread of computers and calculators,
!    many users still suffer from a certain computational naivete.
!    Whatever the computer produces as a result is regarded as
!    almost mathematically proved.  It is often quite astonishing
!    for a user to discover that even in simple calculations with
!    just a few operations it is possible to produce completely
!    incorrect results, and that this is a necessary consequence
!    of current computer techniques (not just the user's method,
!    but the arithmetic unit of the computer as well).
!
!    The following examples are offered to show how few operations
!    it takes to cause a computer to return meaningless results.
!    This fact applies not just to hand calculators
!    but also to the biggest supercomputers.  The particular
!    problem involved in these effects is that one never knows
!    whether, when or where these sorts of errors
!    will enter into a long, involved calculation and destroy the
!    results.
!
!    Only in the last few years have methods and solution
!    processes been developed which can mathematically describe
!    these effects, and numerically control them.
!
!  Modified:
!
!    08 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    U Kulisch, C Ullrich, Editors,
!    Wissenschaftliches Rechnen und Programmiersprachen,
!    (Scientific Computing and Programming Languages),
!    Berichte des German Chapter of the ACM,
!    (Reports of the German Chapter of the ACM)
!    Volume 10, Teubner Verlag, 1982.
!
!    Yves Nievergelt,
!    Numerical Linear Algebra on the HP-28, or How to Lie with Supercalculators,
!    The American Mathematical Monthly,
!    Volume 98, Number 6, June-July 1991, pages 539-544.
!
!    S M Rump,
!    Wie Zuverlaessig Sind die Ergebnisse Unserer Rechenanlagen?
!    (How Reliable are the Results of our Computations?)
!    Jahrbuch Ueberblicke Mathematik 1983, pages 163-168.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ERRORS_PRB'
  write ( *, '(a)' ) '  Carry out the ERROR tests.'

  call header

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test065
  call test07
  call test08
  call test09
  call test10

  call test11
  call test12
  call test13
  call test14
  call test15
  call test16
  call test17
  call test18
  call test19
  call test20

  call test21
  call test215
  call test22
  call test23
  call test24
  call test25
  call test26

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ERRORS_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine header
!
!*******************************************************************************
!
!! HEADER prints out a header.
!
  implicit none
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HOW RELIABLE ARE OUR CALCULATIONS?'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'This program demonstrates how unreliable certain'
  write ( *, '(a)' ) 'calculations are when performed with a finite'
  write ( *, '(a)' ) 'accuracy on a computer. '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Standard software is used where possible for the '
  write ( *, '(a)' ) 'solution of these problems.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Paraphrasing the reference:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Despite the spread of computers and calculators,'
  write ( *, '(a)' ) 'many users suffer a certain computational '
  write ( *, '(a)' ) 'gullibility.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Whatever the computer produces is regarded as'
  write ( *, '(a)' ) 'almost mathematically proven.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'It is often quite astonishing to discover that'
  write ( *, '(a)' ) 'even in simple calculations with just a few '
  write ( *, '(a)' ) 'operations it is possible to produce completely '
  write ( *, '(a)' ) 'incorrect results, and that this is a necessary '
  write ( *, '(a)' ) 'consequence of current computer techniques - not '
  write ( *, '(a)' ) 'just the user''s method, but the arithmetic '
  write ( *, '(a)' ) 'processes of the computer as well.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'The following examples show how few operations '
  write ( *, '(a)' ) 'are needed for a computer to return meaningless '
  write ( *, '(a)' ) 'results.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'This fact applies not just to hand calculators'
  write ( *, '(a)' ) 'but also to the biggest supercomputers.  The '
  write ( *, '(a)' ) 'worst aspect of these effects is that one never'
  write ( *, '(a)' ) 'knows whether, when or where these sorts of '
  write ( *, '(a)' ) 'errors will enter into a long, involved '
  write ( *, '(a)' ) 'calculation and destroy the results.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Only in the last few years have methods and '
  write ( *, '(a)' ) 'solution processes been developed to '
  write ( *, '(a)' ) 'mathematically describe these effects, and '
  write ( *, '(a)' ) 'numerically control them.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '

  return
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01
!
  implicit none
!
  double precision, parameter :: dp = 665857.0D+00
  double precision, parameter :: dq = 470832.0D+00
  double precision dr
  real, parameter :: p = 665857.0E+00
  real, parameter :: q = 470832.0E+00
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 1                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P = 665857                               |' 
  write ( *, '(a)' ) '|    Q = 470832                               |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Compute:                                   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = P**2 - 2 * Q**2                      |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = 1                                    |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '

  r = p**2 - 2.0E+00 * q**2

  write ( *, '(a,g14.6)' ) '  Single precision R = ', r

  dr = dp**2 - 2.0D+00 * dq**2

  write ( *, '(a,g14.6)' ) '  Double precision R = ', dr

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02
!
  implicit none
!
  double precision, parameter :: dp = 10864.0D+00
  double precision, parameter :: dq = 18817.0D+00
  double precision dr
  real, parameter :: p = 10864.0E+00
  real, parameter :: q = 18817.0E+00
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 2                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P = 10864                                |'
  write ( *, '(a)' ) '|    Q = 18817                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Compute:                                   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = 9 * P**4 - Q**4 + 2 * Q**2           |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = 1                                    |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '

  r = 9.0E+00 * p**4 - q**4 + 2.0E+00 * q**2

  write ( *, '(a,g14.6)' ) '  Single precision R = ', r

  dr = 9.0D+00 * dp**4 - dq**4 + 2.0D+00 * dq**2

  write ( *, '(a,g14.6)' ) '  Double precision R = ', dr

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03
!
  implicit none
!
  double precision dp
  double precision dq
  double precision dr
  real p
  real q
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 3                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P = 10**34                               |'
  write ( *, '(a)' ) '|    Q = - 2                                  |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Compute:                                   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = ( P + Q ) - P                        |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = - 2                                  |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '

  p = 1.0E+34
  q = - 2.0E+00

  r = ( p + q ) - p

  write ( *, '(a,g14.6)' ) '  Single precision R = ', r

  dp = 1.0D+34
  dq = - 2.0D+00

  dr = ( dp + dq ) - dp

  write ( *, '(a,g14.6)' ) '  Double precision R = ', dr

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04
!
  implicit none
!
  integer, parameter :: n = 3
!
  real c(0:n)
  double precision dc(0:n)
  double precision dq
  double precision dr
  real q
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 4                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    Q = 1.091608                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Compute:                                   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R =                                      |'
  write ( *, '(a)' ) '|        170.4  * Q**3                        |'
  write ( *, '(a)' ) '|      - 356.41 * Q**2                        |'
  write ( *, '(a)' ) '|      + 168.97 * Q                           |'
  write ( *, '(a)' ) '|      +  18.601                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = 0.821248E-13                         |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '

  q = 1.091608E+00

  c(0) =    18.601E+00
  c(1) =   168.97E+00
  c(2) = - 356.41E+00
  c(3) =   170.4E+00

  call rpoly_val ( n, c, q, r )

  write ( *, '(a,g14.6)' ) '  Single precision R (direct evaluation) = ', r

  call rpoly_val_horner ( n, c, q, r )

  write ( *, '(a,g14.6)' ) '  Single precision R (Horner''s method) =  ', r

  dq = 1.091608D+00

  dc(0) =    18.601D+00
  dc(1) =   168.97D+00
  dc(2) = - 356.41D+00
  dc(3) =   170.4D+00

  call dpoly_val ( n, dc, dq, dr )

  write ( *, '(a,g14.6)' ) '  Double precision R (direct evaluation) = ', dr

  call dpoly_val_horner ( n, dc, dq, dr )

  write ( *, '(a,g14.6)' ) '  Double precision R (Horner''s method) =  ', dr