errors.f90

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

      dsdot = dsdot + dble ( x(ix) ) * dble ( y(iy) )
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  sdsdot = sngl ( dsdot )

  return
end
subroutine sgedi ( a, lda, n, ipvt, det, work, job )

!*****************************************************************************80
!
!! SGEDI computes the determinant and inverse of a matrix factored by SGECO or SGEFA.
!
!  Discussion:
!
!    A division by zero will occur if the input factor contains
!    a zero on the diagonal and the inverse is requested.
!    It will not occur if the subroutines are called correctly
!    and if SGECO has set RCOND > 0.0E+00 or SGEFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real A(LDA,N), on input, the N by N factored matrix.
!    as output by SGECO or SGEFA.  On output, contains the inverse
!    matrix if requested.
!
!    Input, integer LDA, the leading dimension of the array A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer IPVT(N), the pivot vector from SGECO or SGEFA.
!
!    Workspace, real WORK(N).
!
!    Output, real DET(2), the determinant of original matrix if requested.
!    determinant = DET(1) * 10.0**DET(2)
!    with  1.0E+00 <= abs ( DET(1) ) < 10.0E+00
!    or DET(1) == 0.0E+00.
!
!    Input, integer JOB, specifies what is to be computed.
!    11, both determinant and inverse.
!    01, inverse only.
!    10, determinant only.
!
  implicit none

  integer lda
  integer n

  real a(lda,n)
  real det(2)
  integer i
  integer ipvt(n)
  integer j
  integer job
  integer k
  integer kp1
  integer l
  real t
  real, parameter :: ten = 10.0E+00
  real work(n)
!
!  Compute the determinant.
!
  if ( job / 10 /= 0 ) then

     det(1) = 1.0E+00
     det(2) = 0.0E+00

     do i = 1, n

        if ( ipvt(i) /= i ) then
          det(1) = - det(1)
        end if

        det(1) = a(i,i) * det(1)

        if ( det(1) == 0.0E+00 ) then
          exit
        end if

        do while ( abs ( det(1) ) < 1.0E+00 )
          det(1) = ten * det(1)
          det(2) = det(2) - 1.0E+00
        end do

        do while ( abs ( det(1) ) >= ten )
          det(1) = det(1) / ten
          det(2) = det(2) + 1.0E+00
        end do

    end do

  end if
!
!  Compute inverse(U).
!
  if ( mod ( job, 10 ) /= 0 ) then

     do k = 1, n

        a(k,k) = 1.0E+00 / a(k,k)
        t = - a(k,k)
        call sscal ( k-1, t, a(1,k), 1 )

        do j = k+1, n
           t = a(k,j)
           a(k,j) = 0.0E+00
           call saxpy ( k, t, a(1,k), 1, a(1,j), 1 )
        end do

     end do
!
!  Form inverse(U) * inverse(L).
!
     do k = n-1, 1, -1

        do i = k+1, n
           work(i) = a(i,k)
           a(i,k) = 0.0E+00
        end do

        do j = k+1, n
           t = work(j)
           call saxpy ( n, t, a(1,j), 1, a(1,k), 1 )
        end do

        l = ipvt(k)
        if ( l /= k ) then
          call sswap ( n, a(1,k), 1, a(1,l), 1 )
        end if

     end do

  end if

  return
end
subroutine sgefa ( a, lda, n, ipvt, info )

!*****************************************************************************80
!
!! SGEFA factors a real matrix.
!
!  Modified:
!
!    07 March 2001
!
!  Parameters:
!
!    Input/output, real A(LDA,N).
!    On intput, the matrix to be factored.
!    On output, an upper triangular matrix and the multipliers used to 
!    obtain it.  The factorization can be written A=L*U, where L is a product 
!    of permutation and unit lower triangular matrices, and U is upper
!    triangular.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, integer INFO, singularity indicator.
!    0, normal value.
!    K, if U(K,K) == 0.  This is not an error condition for this subroutine,
!    but it does indicate that SGESL or SGEDI will divide by zero if called.
!    Use RCOND in SGECO for a reliable indication of singularity.
!
  implicit none

  integer lda
  integer n

  real a(lda,n)
  integer info
  integer ipvt(n)
  integer isamax
  integer j
  integer k
  integer l
  real t
!
!  Gaussian elimination with partial pivoting.
!
  info = 0

  do k = 1, n - 1
!
!  Find L = pivot index.
!
    l = isamax ( n-k+1, a(k,k), 1 ) + k - 1
    ipvt(k) = l
!
!  Zero pivot implies this column already triangularized.
!
    if ( a(l,k) == 0.0E+00 ) then
      info = k
      cycle
    end if
!
!  Interchange if necessary.
!
    if ( l /= k ) then
      t = a(l,k)
      a(l,k) = a(k,k)
      a(k,k) = t
    end if
!
!  Compute multipliers.
!
    t = -1.0E+00 / a(k,k)
    call sscal ( n-k, t, a(k+1,k), 1 )
!
!  Row elimination with column indexing.
!
    do j = k+1, n
      t = a(l,j)
      if ( l /= k ) then
        a(l,j) = a(k,j)
        a(k,j) = t
      end if
      call saxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 )
    end do

  end do

  ipvt(n) = n

  if ( a(n,n) == 0.0E+00 ) then
    info = n
  end if

  return
end
subroutine sgesl ( a, lda, n, ipvt, b, job )

!*****************************************************************************80
!
!! SGESL solves a real general linear system A * X = B.
!
!  Discussion:
!
!    SGESL can solve either of the systems A * X = B or transpose ( A ) * X = B.
!
!    The system matrix must have been factored by SGECO or SGEFA.
!
!    A division by zero will occur if the input factor contains a
!    zero on the diagonal.  Technically this indicates singularity
!    but it is often caused by improper arguments or improper
!    setting of LDA.  It will not occur if the subroutines are
!    called correctly and if SGECO has set RCOND > 0.0E+00
!    or SGEFA has set INFO == 0.
!
!  Modified:
!
!    07 March 2001
!
!  Parameters:
!
!    Input, real A(LDA,N), the output from SGECO or SGEFA.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer IPVT(N), the pivot vector from SGECO or SGEFA.
!
!    Input/output, real B(N).
!    On input, the right hand side vector.
!    On output, the solution vector.
!
!    Input, integer JOB.
!    0, solve A * X = B;
!    nonzero, solve A' * X = B.
!
  implicit none

  integer lda
  integer n

  real a(lda,n)
  real b(n)
  integer ipvt(n)
  integer job
  integer k
  integer l
  real sdot
  real t
!
!  Solve A * X = B.
!
  if ( job == 0 ) then

    do k = 1, n-1

      l = ipvt(k)
      t = b(l)

      if ( l /= k ) then
        b(l) = b(k)
        b(k) = t
      end if

      call saxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 )

    end do

    do k = n, 1, -1
      b(k) = b(k) / a(k,k)
      t = -b(k)
      call saxpy ( k-1, t, a(1,k), 1, b(1), 1 )
    end do

  else
!
!  Solve A' * X = B.
!
    do k = 1, n
      t = sdot ( k-1, a(1,k), 1, b(1), 1 )
      b(k) = ( b(k) - t ) / a(k,k)
    end do

    do k = n-1, 1, -1

      b(k) = b(k) + sdot ( n-k, a(k+1,k), 1, b(k+1), 1 )
      l = ipvt(k)

      if ( l /= k ) then
        t = b(l)
        b(l) = b(k)
        b(k) = t
      end if

    end do

  end if

  return
end
function rvec_norm2 ( n, x, incx )

!*****************************************************************************80
!
!! RVEC_NORM2 computes the Euclidean norm of a vector.
!
!  Discussion:
!
!    The original SNRM2 algorithm is accurate but written in a bizarre,
!    unreadable and obsolete format.  This version goes for clarity.
!
!  Modified:
!
!    01 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector whose norm is to be computed.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SNRM2, the Euclidean norm of X.
!
  implicit none

  integer i
  integer incx
  integer ix
  integer n
  real samax
  real rvec_norm2
  real stemp
  real x(*)
  real xmax

  if ( n <= 0 ) then

    rvec_norm2 = 0.0E+00

  else

    xmax = samax ( n, x, incx )

    if ( xmax == 0.0E+00 ) then

      rvec_norm2 = 0.0E+00

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = ( - n + 1 ) * incx + 1
      end if

      stemp = 0.0E+00
      do i = 1, n
        stemp = stemp + ( x(ix) / xmax )**2
        ix = ix + incx
      end do

      rvec_norm2 = xmax * sqrt ( stemp )

    end if

  end if

  return
end
subroutine sqrdc ( a, lda, n, p, qraux, jpvt, work, job )

!*****************************************************************************80
!
!! SQRDC computes the QR factorization of a real rectangular matrix.
!
!  Discussion:
!
!    SQRDC uses Householder transformations.
!
!    Column pivoting based on the 2-norms of the reduced columns may be
!    performed at the user's option.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real A(LDA,P).  On input, the N by P matrix
!    whose decomposition is to be computed.  On output, A contains in
!    its upper triangle the upper triangular matrix R of the QR
!    factorization.  Below its diagonal A contains information from
!    which the orthogonal part of the decomposition can be recovered.
!    Note that if pivoting has been requested, the decomposition is not that
!    of the original matrix A but that of A with its columns permuted
!    as described by JPVT.
!
!    Input, integer LDA, the leading dimension of the array A.  LDA must
!    be at least N.
!
!    Input, integer N, the number of rows of the matrix A.
!
!    Input, integer P, the number of columns of the matrix A.
!
!    Output, real QRAUX(P), contains further information required to recover
!    the orthogonal part of the decomposition.
!
!    Input/output, integer JPVT(P).  On input, JPVT contains integers that
!    control the selection of the pivot columns.  The K-th column A(*,K) of A
!    is placed in one of three classes according to the value of JPVT(K).
!      > 0, then A(K) is an initial column.
!      = 0, then A(K) is a free column.
!      < 0, then A(K) is a final column.
!    Before the decomposition is computed, initial columns are moved to
!    the beginning of the array A and final columns to the end.  Both
!    initial and final columns are frozen in place during the computation
!    and only free columns are moved.  At the K-th stage of the
!    reduction, if A(*,K) is occupied by a free column it is interchanged
!    with the free column of largest reduced norm.  JPVT is not referenced
!    if JOB == 0.  On output, JPVT(K) contains the index of the column of the
!    original matrix that has been interchanged into the K-th column, if
!    pivoting was requested.
!
!    Workspace, real WORK(P).  WORK is not referenced if JOB == 0.
!
!    Input, integer JOB, initiates column pivoting.
!    0, no pivoting is done.
!    nonzero, pivoting is done.
!
  implicit none

  integer lda
  integer n
  integer p

  real a(lda,p)
  integer jpvt(p)
  real qraux(p)
  real work(p)
  integer j
  integer jj
  integer job
  integer jp
  integer l
  integer lup
  integer maxj
  real maxnrm
  real nrmxl
  integer pl
  integer pu
  real sdot
  real rvec_norm2
  logical swapj
  real t
  real tt

  pl = 1
  pu = 0
!
!  If pivoting is requested, rearrange the columns.
!
  if ( job /= 0 ) then

    do j = 1, p

      swapj = jpvt(j) > 0

      if ( jpvt(j) < 0 ) then
        jpvt(j) = - j
      else
        jpvt(j) = j
      end if

      if ( swapj ) then

        if ( j /= pl ) then
          call sswap ( n, a(1,pl), 1, a(1,j), 1 )
        end if

        jpvt(j) = jpvt(pl)
        jpvt(pl) = j
        pl = pl + 1

      end if

    end do

    pu = p

    do j = p, 1, -1

      if ( jpvt(j) < 0 ) then

        jpvt(j) = - jpvt(j)

        if ( j /= pu ) then
          call sswap ( n, a(1,pu), 1, a(1,j), 1 )
          jp = jpvt(pu)
          jpvt(pu) = jpvt(j)
          jpvt(j) = jp
        end if

        pu = pu - 1

      end if

    end do

  end if
!
!  Compute the norms of the free columns.
!
  do j = pl, pu
    qraux(j) = rvec_norm2 ( n, a(1,j), 1 )
  end do

  work(pl:pu) = qraux(pl:pu)
!
!  Perform the Householder reduction of A.
!
  lup = min ( n, p )

  do l = 1, lup
!
!  Bring the column of largest norm into the pivot position.
!
    if ( l >= pl .and. l < pu ) then

      maxnrm = 0.0E+00
      maxj = l
      do j = l, pu
        if ( qraux(j) > maxnrm ) then
          maxnrm = qraux(j)
          maxj = j
        end if
      end do

      if ( maxj /= l ) then
        call sswap ( n, a(1,l), 1, a(1,maxj), 1 )
        qraux(maxj) = qraux(l)
        work(maxj) = work(l)
        jp = jpvt(maxj)
        jpvt(maxj) = jpvt(l)
        jpvt(l) = jp
      end if

    end if
!
!  Compute the Householder transformation for column L.
!
    qraux(l) = 0.0E+00

    if ( l /= n ) then