errors_prb.f90

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

    else if ( i == 2 ) then
      h = 1.0E-05
    else if ( i == 3 ) then
      h = 1.0E-08
    end if

    r = ( fx15 ( x + h ) - 2.0E+00 * fx15 ( x ) + fx15 ( x - h ) ) / h**2

    write ( *, '(2g14.6)' ) h, r

  end do

  return
end
subroutine test16
!
!*******************************************************************************
!
!! TEST16
!
  implicit none
!
  real scale
  real x
  real y
  real z
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 16                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  P(X,Y) =                                   |'
  write ( *, '(a)' ) '|     83521 *        y**8                     |'
  write ( *, '(a)' ) '|    +  578 * x**2 * y**4                     |'
  write ( *, '(a)' ) '|    -    2 * x**4                            |'
  write ( *, '(a)' ) '|    +    2 * x**6                            |'
  write ( *, '(a)' ) '|    -        x**8                            |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Evaluate P(X,Y) for                        |' 
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    X = 9478657                              |'
  write ( *, '(a)' ) '|    Y = 2298912                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P(X,Y) = - 179,689,877,047,297           |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'

  scale = 1000000.0E+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The unscaled calculation will overflow.'
  write ( *, '(a,g14.6)' ) '  We use a scale factor on X and Y of ', scale

  x = 9478657.0E+00 / scale
  y = 2298912.0E+00 / scale

  z =     83521.0E+00 *        y**8 &
        +   578.0E+00 * x**2 * y**4 / scale**2 &
        -     2.0E+00 * x**4 / scale**4 &
        +     2.0E+00 * x**6 * scale**2 &
        -           x**8
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Scaled value of P(X,Y) = ', z
  write ( *, '(a,g14.6,a,g14.6,a)' ) '  Computed value of P(X,Y) = ', z , &
    ' * ', scale, ' ** 8'

  return
end
subroutine test17
!
!*******************************************************************************
!
!! TEST17
!
  implicit none
!
  integer, parameter :: n = 5
!
  real, dimension ( n ) :: a = (/ 2.718281828E+00, -3.141592654E+00, &
    1.414213562E+00, 0.5772156649E+00, 0.3010299957E+00 /)
  real b(n)
  double precision, dimension ( n ) :: da = (/ 2.718281828D+00, &
    -3.141592654D+00, 1.414213562D+00, 0.5772156649D+00, &
    0.3010299957D+00 /)
  double precision db(n)
  double precision dr
  integer i
  real r
  real sdot
  real sdsdot
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 17                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  A =                                        |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|     2.718281828                             |'
  write ( *, '(a)' ) '|    -3.141592654                             |'
  write ( *, '(a)' ) '|     1.414213562                             |'
  write ( *, '(a)' ) '|     0.5772156649                            |'
  write ( *, '(a)' ) '|     0.3010299957                            |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  B =                                        |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|        1486.2497                            |'
  write ( *, '(a)' ) '|      878366.9879                            |'
  write ( *, '(a)' ) '|        - 22.37492                           |'
  write ( *, '(a)' ) '|     4773714.647                             |'
  write ( *, '(a)' ) '|           0.000185049                       |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Compute the scalar product:                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = A dot B                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    R = - 0.100657107E-10                    |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '

  b(1) =    1486.2497E+00
  b(2) =  878366.9879E+00
  b(3) =    - 22.37492E+00
  b(4) = 4773714.647E+00
  b(5) =       0.000185049E+00

  r = 0.0E+00
  do i = 1, n
    r = r + a(i) * b(i)
  end do

  write ( *, '(a,g14.6)' ) '  Standard loop     R = ', r

  r = sdot ( n, a, 1, b, 1 )
  write ( *, '(a,g14.6)' ) '  SDOT:             R = ', r

  r = sdsdot ( n, a, 1, b, 1 )
  write ( *, '(a,g14.6)' ) '  SDSDOT:           R = ', r

  r = dot_product ( a, b )
  write ( *, '(a,g14.6)' ) '  DOT_PRODUCT:      R = ', r

  db(1) =    1486.2497D+00
  db(2) =  878366.9879D+00
  db(3) =    - 22.37492D+00
  db(4) = 4773714.647D+00
  db(5) =       0.000185049D+00

  dr = 0.0E+00
  do i = 1, n
    dr = dr + da(i) * db(i)
  end do

  write ( *, '(a,g14.6)' ) '  Standard loop    DR = ', dr
  dr = dot_product ( da, db )

  write ( *, '(a,g14.6)' ) '  DOT_PRODUCT:     DR = ', dr

  return
end
subroutine test18
!
!*******************************************************************************
!
!! TEST18
!
  implicit none
!
  real scale
  real x
  real y
  real z
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 18                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  X = 192119201                              |'
  write ( *, '(a)' ) '|  Y = 35675640                               |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Z(X,Y) =                                   |'
  write ( *, '(a)' ) '|    (                                        |'
  write ( *, '(a)' ) '|         1682 * x    * y**4                  |'
  write ( *, '(a)' ) '|       +    3 * x**3                         |'
  write ( *, '(a)' ) '|       +   29 * x    * y**2                  |'
  write ( *, '(a)' ) '|       -    2 * x**5                         |'
  write ( *, '(a)' ) '|       +  832                                |'
  write ( *, '(a)' ) '|      ) / 107751                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct Z:                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    1783.0                                   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'

  scale = 1000000.0E+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The unscaled calculation will overflow.'
  write ( *, '(a,g14.6)' ) '  We use a scale factor on X and Y of ', scale

  x = 192119201.0E+00 / scale
  y = 35675640.0E+00 / scale

  z = (      1682.0E+00 * x    * y**4 &
           +    3.0E+00 * x**3 / scale**2 &
           +   29.0E+00 * x    * y**2 / scale**2 &
           -    2.0E+00 * x**5 &
           +  832.0E+00 / scale**5 &
         ) / 107751.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Scaled value of Z = ', z
  write ( *, '(a,g14.6,a,g14.6,a,i5)' ) '  Computed value of Z = ', z, &
    ' * ', scale, ' ** ', 5

  return
end
subroutine test19
!
!*******************************************************************************
!
!! TEST19
!
  implicit none
!
  real fx1
  real fx2
  real fx3
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 19                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    F1(X) = 1 - Cos(X)                       |'
  write ( *, '(a)' ) '|    F2(X) = Sin**2(X) / ( 1 + Cos(X) )       |'
  write ( *, '(a)' ) '|    F3(X) = Taylor series for 1 - Cos(X)     |'
  write ( *, '(a)' ) '|      around X = 0, up to X**6 terms.        |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  For all X                                  |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    F1(X) = F2(X)                            |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  For X near 0,                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    F3(X) should approximate F1(X)           |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       X          1-Cos        Sin**2/(1+Cos)  Taylor'
  write ( *, '(a)' ) ' '

  x = 0.5E+00

  do

    fx1 = 1.0E+00 - cos ( x )
    fx2 = ( sin ( x ) )**2 / ( 1.0E+00 + cos ( x ) )
    fx3 = 0.5E+00 * x**2 * ( 1.0E+00 - ( x**2 / 12.0E+00 ) + ( x**4 / 360.0E+00 ) )

    write ( *, '(4g14.6)' ) x, fx1, fx2, fx3

    x = 0.25E+00 * x
 
    if ( x <= 1.0E-08 ) then
      exit
    end if

  end do

  return
end
subroutine test20
!
!*******************************************************************************
!
!! TEST20
!
  implicit none
!
  integer, parameter :: n = 2
!
  double precision, dimension (0:n) :: dp = &
    (/ 0.005D+00, 100.0D+00, 0.00005D+00 /)
  double precision dpval(2)
  complex ( kind = 16 ) dr(2)
  integer ierror
  real, dimension (0:n) :: p = (/ 0.005E+00, 100.0E+00, 0.00005E+00 /)
  real pval(2)
  complex r(2)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 20                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  P(X) =                                     |'
  write ( *, '(a)' ) '|        0.00005 * x**2                       |'
  write ( *, '(a)' ) '|    + 100        * x                         |'
  write ( *, '(a)' ) '|    +   0.005                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Evaluate the standard quadratic formula:   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    X1 = 0.5 * ( -B + sqrt(B**2-4*A*C) ) / A |'
  write ( *, '(a)' ) '|    X2 = 0.5 * ( -B - sqrt(B**2-4*A*C) ) / A |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Evaluate an alternate quadratic formula:   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    X3 = 2 * C / ( - B + sqrt(B**2-4*A*C) )  |'
  write ( *, '(a)' ) '|    X4 = 2 * C / ( - B - sqrt(B**2-4*A*C) )  |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct results:                           |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    Root1 =       -0.00005...                |'
  write ( *, '(a)' ) '|    Root2 = -1999999.99995...                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'

  call rpoly2_roots ( p, r )

  call rpoly_val_horner ( n, p, real ( r(1) ), pval(1) )
  call rpoly_val_horner ( n, p, real ( r(2) ), pval(2) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) '  X1 = ', real ( r(1) ), &
    ' with P(X1) = ', pval(1)
  write ( *, '(a,g14.6,a,g14.6)' ) '  X2 = ', real ( r(2) ), &
    ' with P(X2) = ', pval(2)

  call rpoly2_roots2 ( p, r, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  RPOLY2_ROOTS2 returned an error flag.'
  else
    call rpoly_val_horner ( n, p, real ( r(1) ), pval(1) )
    call rpoly_val_horner ( n, p, real ( r(2) ), pval(2) )
    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6,a,g14.6)' ) '  X3 = ', real ( r(1) ), &
      ' with P(X3) = ', pval(1)
    write ( *, '(a,g14.6,a,g14.6)' ) '  X4 = ', real ( r(2) ), &
      ' with P(X4) = ', pval(2)
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Repeat calculations in double precision:'
  write ( *, '(a)' ) ' '

  call dpoly2_roots ( dp, dr )

  call dpoly_val_horner ( n, dp, dble ( dr(1) ), dpval(1) )
  call dpoly_val_horner ( n, dp, dble ( dr(2) ), dpval(2) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g20.12,a,g14.6)' ) '  X1 = ', dble ( dr(1) ), &
    ' with P(X1) = ', dpval(1)
  write ( *, '(a,g20.12,a,g14.6)' ) '  X2 = ', dble ( dr(2) ), &
    ' with P(X2) = ', dpval(2)

  call dpoly2_roots2 ( dp, dr, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  DPOLY2_ROOTS2 returned an error flag.'
  else
    call dpoly_val_horner ( n, dp, dble ( dr(1) ), dpval(1) )
    call dpoly_val_horner ( n, dp, dble ( dr(2) ), dpval(2) )
    write ( *, '(a)'