errors_prb.f90

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

    call sgesl ( a, lda, n, ipivot, r, job )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  SGEFA/SGESL solution (Gauss elimination):'
    write ( *, '(a)' ) ' '

    do i = 1, n
      write ( *, '(g14.6)' ) r(i)
    end do

  end if

  a(1,1) =  -367296.0E+00
  a(1,2) =   -43199.0E+00
  a(1,3) =   519436.0E+00
  a(1,4) =  -954302.0E+00
  a(2,1) =   259718.0E+00
  a(2,2) =  -477151.0E+00
  a(2,3) =  -367295.0E+00
  a(2,4) = -1043199.0E+00
  a(3,1) =   886731.0E+00
  a(3,2) =    88897.0E+00
  a(3,3) = -1254026.0E+00
  a(3,4) = -1132096.0E+00
  a(4,1) =   627013.0E+00
  a(4,2) =   566048.0E+00
  a(4,3) =  -886732.0E+00
  a(4,4) =   911103.0E+00

  r(1:4) = (/ 1.0E+00, 1.0E+00, 1.0E+00, 0.0E+00 /)

  job = 0

  call sqrdc ( a, lda, n, n, qraux, ipivot, work, job )

  job = 110

  call sqrsl ( a, lda, n, n, qraux, r, qy, qty, r, rsd, ab, job, info )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  SQRDC/SQRSL solution (QR factorization):'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(g14.6)' ) r(i)
  end do

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08
!
  implicit none
!
  real a
  real b
  real f1
  real f2
  real f3
  real f4
  real fmin
  real, external :: p8horn
  real, external :: p8reg
  real r1
  real r2
  real r3
  real r4
  real, parameter :: tol = 1.0E-14
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 8                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P(X) =                                   |'
  write ( *, '(a)' ) '|        170.4  * x**3                        |'
  write ( *, '(a)' ) '|      - 356.41 * x**2                        |'
  write ( *, '(a)' ) '|      + 168.97 * x                           |'
  write ( *, '(a)' ) '|      +  18.601                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Seek the minimizer X* of P(X)              |'
  write ( *, '(a)' ) '|  over the interval 0 <= X <= 2.             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  The true minimizer lies in the interval    |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  [1.091607978, 1.091607981 ]                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'

  a = 0.0E+00
  b = 2.0E+00

  r1 = fmin ( a, b, p8reg, tol )
  f1 = p8reg ( r1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Evaluating P(X) in the usual way'
  write ( *, '(a,g14.6)' ) '  FMIN computes the minimizer at ', r1
  write ( *, '(a,g14.6)' ) '  with function value ', f1

  r2 = fmin ( a, b, p8horn, tol )
  f2 = p8horn ( r2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Evaluating P(X) using Horner''s method,'
  write ( *, '(a,g14.6)' ) '  FMIN computes the minimizer at ',r2
  write ( *, '(a,g14.6)' ) '  with function value ', f2

  r3 = 1.091607978E+00
  r4 = 1.091607981E+00
  f3 = p8horn ( r1 )
  f4 = p8horn ( r2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) '  The true value lies between ', r3, &
    ' and ', r4
  write ( *, '(a,2g14.6)' ) '  with function value below ', f3, f4

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09
!
  implicit none
!
  integer, parameter :: n = 3
!
  real, dimension (0:n) :: c = (/ -945804881.0E+00, 1753426039.0E+00, &   
    -1083557822.0E+00, 223200658.0E+00 /)
  double precision, dimension (0:n) :: dc = (/ -945804881.0D+00, &
    1753426039.0D+00, -1083557822.0D+00, 223200658.0D+00 /)
  double precision dphorn
  double precision dpreg
  double precision dx
  integer i
  real phorn
  real preg
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 9                                 |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    P(X) =                                   |'
  write ( *, '(a)' ) '|         223200658 * X**3                    |'
  write ( *, '(a)' ) '|      - 1083557822 * X**2                    |'
  write ( *, '(a)' ) '|      + 1753426039 * X                       |'
  write ( *, '(a)' ) '|       - 945804881                           |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Evaluate at 11 equally spaced values       |'
  write ( *, '(a)' ) '|  from 1.61801916 to 1.61801917              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct table:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  ===== X =====    ========= P(X) ========   |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    1.618019160  - 0.17081105112320e-11      |'
  write ( *, '(a)' ) '|    1.618019161  - 0.89804011575510e-12      |'
  write ( *, '(a)' ) '|    1.618019162  - 0.34596536943057e-12      |'
  write ( *, '(a)' ) '|    1.618019163  - 0.51884933054474e-13      |'
  write ( *, '(a)' ) '|    1.618019164  - 0.15797467422848e-13      |'
  write ( *, '(a)' ) '|    1.618019165  - 0.23770163333175e-12      |'
  write ( *, '(a)' ) '|    1.618019166  - 0.71759609157723e-12      |'
  write ( *, '(a)' ) '|    1.618019167  - 0.14554795029553e-11      |'
  write ( *, '(a)' ) '|    1.618019168  - 0.24513505282621e-11      |'
  write ( *, '(a)' ) '|    1.618019169  - 0.37052078282936e-11      |'
  write ( *, '(a)' ) '|    1.618019170  - 0.52170500638460e-11      |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  R1: compute P(X) the regular way.'
  write ( *, '(a)' ) '  R2: compute P(X) using Horner''s method.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Single precision:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X, R1, R2'
  write ( *, '(a)' ) ' '

  do i = 0, 10

    x = ( real ( 10 - i ) * 1.61801916E+00 &
        + real (      i ) * 1.61801917E+00 ) / 10.0E+00

    call rpoly_val ( n, c, x, preg )

    call rpoly_val_horner ( n, c, x, phorn )

    write ( *, '(f14.10,3g14.6)' ) x, preg, phorn

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Double precision:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X, DR1, DR2'
  write ( *, '(a)' ) ' '
  do i = 0, 10

    dx = ( dble ( 10 - i ) * 1.61801916D+00 &
      +  dble ( i ) * 1.61801917D+00 ) / 10.0D+00

    dpreg = dc(3) * dx**3 + dc(2) * dx**2 + dc(1) * dx + dc(0)

    call dpoly_val_horner ( n, dc, dx, dphorn )

    write ( *, '(f14.10,3g14.6)' ) dx, dpreg, dphorn

  end do

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 seeks the nearest straight line to given data.
!
  implicit none
!
  double precision drb
  double precision drm
  double precision drm_bot
  double precision drm_top
  double precision drv
  double precision dx
  double precision dx1
  double precision dx2
  double precision dx3
  double precision dy1
  double precision dy2
  double precision dy3
  real rb
  real rm
  real rm_bot
  real rm_top
  real rv
  real x
  real x1
  real x2
  real x3
  real y1
  real y2
  real y3
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 10                                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  For the data:                              |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|       X       Y                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    5201477   99999                          |'
  write ( *, '(a)' ) '|    5201478  100000                          |'
  write ( *, '(a)' ) '|    5201479  100001                          |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Find the "nearest" straight line:          |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    Y(X) = A * X + B                         |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  in the least squares sense.                |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct values:                            |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    A = 1.0                                  |'
  write ( *, '(a)' ) '|    B = - 5101478                            |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Evaluate:                                  |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    Y(5201480)                               |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Correct value:                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|    Y(5201480) = 100002.0                    |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '                    Slope A    Intercept B    Y(5201480) '
  write ( *, '(a)' ) ' '

  x1 = 5201477.0E+00
  x2 = 5201478.0E+00
  x3 = 5201479.0E+00

  y1 =  99999.0E+00
  y2 = 100000.0E+00
  y3 = 100001.0E+00

  rm_top = ( x1 * y1 + x2 * y2 + x3 * y3 &
         - ( x1 + x2 + x3 ) * ( y1 + y2 + y3 ) / 3.0E+00 )
  rm_bot = ( x1**2 + x2**2 + x3**2 - ( x1 + x2 + x3 )**2 / 3.0E+00 )

  if ( rm_bot == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Single precision slope calculation fails.'
    write ( *, '(a)' ) '  The divisor was computed as 0.'
  else

    rm = rm_top / rm_bot

    rb = ( y1 + y2 + y3 ) / 3.0E+00 - rm * ( x1 + x2 + x3 ) / 3.0E+00

    x = 5201480.0E+00
    rv = rm * x + rb

    write ( *, '(a,3g14.6)' ) '  Single precision ', rm, rb, rv

  end if

  dx1 = 5201477.0D+00
  dx2 = 5201478.0D+00
  dx3 = 5201479.0D+00

  dy1 = 99999.0D+00
  dy2 = 100000.0D+00
  dy3 = 100001.0D+00

  drm_top = ( dx1 * dy1 + dx2 * dy2 + dx3 * dy3 &
          - ( dx1 + dx2 + dx3 ) * ( dy1 + dy2 + dy3 ) / 3.0D+00 )
  drm_bot = ( dx1**2 + dx2**2 + dx3**2 - ( dx1 + dx2 + dx3 )**2 / 3.0D+00 )

  if ( drm_bot == 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Double precision slope calculation fails.'
    write ( *, '(a)' ) '  The divisor was computed as 0.'
  else
    drm = drm_top / drm_bot
    drb = ( dy1 + dy2 + dy3 ) / 3.0D+00 - drm * ( dx1 + dx2 + dx3 ) / 3.0D+00

    dx = 5201480.0D+00
    drv = drm * dx + drb

    write ( *, '(a,3g14.6)' ) '  Double precision ', drm, drb, drv
  end if

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11
!
  implicit none
!
  real a
  real b
  real f1
  real f2
  real, external :: fhorn11
  real fmin
  real, external :: freg11
  real r1
  real r2
  real :: tol = 1.0E-14
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '+=============================================+'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  Exercise 11                                |' 
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  P(X) =                                     |'
  write ( *, '(a)' ) '|      2124476931 * X**4                      |'
  write ( *, '(a)' ) '|    - 1226567328 * X**3                      |'
  write ( *, '(a)' ) '|    -  708158977 * X**2                      |'
  write ( *, '(a)' ) '|    +  408855776 * X                         |'
  write ( *, '(a)' ) '|    +    1.0E-27                             |'
  write ( *, '(a)' ) '|                                             |'
  write ( *, '(a)' ) '|  P(X) is positive for X >= 0.               |'