dqed.f90

求解非线性方程组的一个高效算法

!
!    Input, integer N, the number of entries in the vector.
!
!    Input, double precision X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, double precision DAMAX, the maximum absolute value of an 
!    element of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  double precision damax
  double precision x(*)
!
  if ( n <= 0 ) then

    damax = 0.0D+00

  else if ( n == 1 ) then

    damax = abs ( x(1) )

  else if ( incx == 1 ) then

    damax = abs ( x(1) )

    do i = 2, n
      if ( abs ( x(i) ) > damax ) then
        damax = abs ( x(i) )
      end if
    end do

  else

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

    damax = abs ( x(ix) )
    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > damax ) then
        damax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function dasum ( n, x, incx )
!
!*******************************************************************************
!
!! DASUM sums the absolute values of the entries of a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, double precision X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, double precision DASUM, the sum of the absolute values of 
!    the vector.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer m
  integer n
  double precision dasum
  double precision stemp
  double precision x(*)
!
  stemp = 0.0D+00

  if ( n <= 0 ) then

  else if ( incx == 1 ) then

    m = mod ( n, 6 )

    do i = 1, m
      stemp = stemp + abs ( x(i) )
    end do

    do i = m+1, n, 6
      stemp = stemp + abs ( x(i)   ) + abs ( x(i+1) ) + abs ( x(i+2) ) &
                    + abs ( x(i+3) ) + abs ( x(i+4) ) + abs ( x(i+5) )
    end do

  else

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

    do i = 1, n
      stemp = stemp + abs ( x(ix) )
      ix = ix + incx
    end do

  end if

  dasum = stemp

  return
end
subroutine daxpy ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! DAXPY adds a constant times one vector to another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, double precision SA, the multiplier.
!
!    Input, double precision X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, double precision Y(*), the vector to which a multiple 
!    of X is to be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  double precision sa
  double precision x(*)
  double precision y(*)
!
  if ( n <= 0 ) then

  else if ( sa == 0.0D+00 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = y(1:n) + sa * x(1:n)

  else

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

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine dbocls ( w, mdw, mcon, mrows, ncols, bl, bu, ind, iopt, x, &
  rnormc, rnorm, mode, rw, iw )
!
!***********************************************************************
!
!! DBOCLS solves a bounded and constrained least squares problem.
!
!
!  Discussion:
!
!    DBOCLS solves the bounded and constrained least squares
!    problem consisting of solving the equation
!      E*X = F  (in the least squares sense)
!    subject to the linear constraints
!      C*X = Y.
!
!     This subprogram solves the bounded and constrained least squares
!     problem. The problem statement is:
!
!     Solve E*X = F (least squares sense), subject to constraints
!     C*X=Y.
!
!     In this formulation both X and Y are unknowns, and both may
!     have bounds on any of their components.  This formulation
!     of the problem allows the user to have equality and inequality
!     constraints as well as simple bounds on the solution components.
!
!     This constrained linear least squares subprogram solves E*X=F
!     subject to C*X=Y, where E is MROWS by NCOLS, C is MCON by NCOLS.
!
!      The user must have dimension statements of the form
!
!      DIMENSION W(MDW,NCOLS+MCON+1), BL(NCOLS+MCON), BU(NCOLS+MCON),
!     * X(2*(NCOLS+MCON)+2+NX), RW(6*NCOLS+5*MCON)
!       integer IND(NCOLS+MCON), IOPT(17+NI), IW(2*(NCOLS+MCON))
!
!     (here NX=number of extra locations required for the options; NX=0
!     if no options are in use. Also NI=number of extra locations
!     for options 1-9.
!
!  Author:
!
!    Richard Hanson, Sandia National Laboratory
!
!  Reference:
!
!    HANSON, R. J. 
!    LINEAR LEAST SQUARES WITH BOUNDS AND LINEAR CONSTRAINTS, 
!    SIAM J. SCI. STAT. COMPUT.,
!    VOL. 7, NO. 3, JULY, 1986.
!
!    INPUT
!    -----
!
!    -------------------------
!    W(MDW,*),MCON,MROWS,NCOLS
!    -------------------------
!     The array W contains the (possibly null) matrix [C:*] followed by
!     [E:F].  This must be placed in W as follows:
!          [C  :  *]
!     W  = [       ]
!          [E  :  F]
!     The (*) after C indicates that this data can be undefined. The
!     matrix [E:F] has MROWS rows and NCOLS+1 columns. The matrix C is
!     placed in the first MCON rows of W(*,*) while [E:F]
!     follows in rows MCON+1 through MCON+MROWS of W(*,*). The vector F
!     is placed in rows MCON+1 through MCON+MROWS, column NCOLS+1. The
!     values of MDW and NCOLS must be positive; the value of MCON must
!     be nonnegative. An exception to this occurs when using option 1
!     for accumulation of blocks of equations. In that case MROWS is an
!     OUTPUT variable only, and the matrix data for [E:F] is placed in
!     W(*,*), one block of rows at a time. See IOPT(*) contents, option
!     number 1, for further details. The row dimension, MDW, of the
!     array W(*,*) must satisfy the inequality:
!
!     If using option 1,
!                     MDW >= MCON + max(max. number of
!                     rows accumulated, NCOLS)
!
!     If using option 8, MDW >= MCON + MROWS.
!     Else, MDW >= MCON + max(MROWS, NCOLS).
!
!     Other values are errors, but this is checked only when using
!     option=2.  The value of MROWS is an output parameter when
!     using option number 1 for accumulating large blocks of least
!     squares equations before solving the problem.
!     See IOPT(*) contents for details about option 1.
!
!    ------------------
!    BL(*),BU(*),IND(*)
!    ------------------
!     These arrays contain the information about the bounds that the
!     solution values are to satisfy. The value of IND(J) tells the
!     type of bound and BL(J) and BU(J) give the explicit values for
!     the respective upper and lower bounds on the unknowns X and Y.
!     The first NVARS entries of IND(*), BL(*) and BU(*) specify
!     bounds on X; the next MCON entries specify bounds on Y.
!
!    1.    For IND(J)=1, require X(J) >= BL(J);
!          IF J > NCOLS,        Y(J-NCOLS) >= BL(J).
!          (the value of BU(J) is not used.)
!    2.    For IND(J)=2, require X(J) <= BU(J);
!          IF J > NCOLS,        Y(J-NCOLS) <= BU(J).
!          (the value of BL(J) is not used.)
!    3.    For IND(J)=3, require X(J) >= BL(J) and
!                                X(J) <= BU(J);
!          IF J > NCOLS,        Y(J-NCOLS) >= BL(J) and
!                                Y(J-NCOLS) <= BU(J).
!          (to impose equality constraints have BL(J)=BU(J)=
!          constraining value.)
!    4.    For IND(J)=4, no bounds on X(J) or Y(J-NCOLS) are required.
!          (the values of BL(J) and BU(J) are not used.)
!
!     Values other than 1,2,3 or 4 for IND(J) are errors. In the case
!     IND(J)=3 (upper and lower bounds) the condition BL(J)  >  BU(J)
!     is  an  error.   The values BL(J), BU(J), J  >  NCOLS, will be
!     changed.  Significant changes mean that the constraints are
!     infeasible.  (Users must make this decision themselves.)
!     The new values for BL(J), BU(J), J  >  NCOLS, define a
!     region such that the perturbed problem is feasible.  If users
!     know that their problem is feasible, this step can be skipped
!     by using option number 8 described below.
!
!    -------
!    IOPT(*)
!    -------
!     This is the array where the user can specify nonstandard options
!     for DBOCLS( ). Most of the time this feature can be ignored by
!     setting the input value IOPT(1)=99. Occasionally users may have
!     needs that require use of the following subprogram options. For
!     details about how to use the options see below: IOPT(*) CONTENTS.
!
!     Option Number   Brief Statement of Purpose
!     ------ ------   ----- --------- -- -------
!           1         Return to user for accumulation of blocks
!                     of least squares equations.  The values
!                     of IOPT(*) are changed with this option.
!                     The changes are updates to pointers for
!                     placing the rows of equations into position
!                     for processing.
!           2         Check lengths of all arrays used in the
!                     subprogram.
!           3         Column scaling of the data matrix, [C].
!                                                        [E]
!           4         User provides column scaling for matrix [C].
!                                                             [E]
!           5         Provide option array to the low-level
!                     subprogram DBOLS( ).
!                     {Provide option array to the low-level
!                     subprogram DBOLSM( ) by imbedding an
!                     option array within the option array to
!                     DBOLS(). Option 6 is now disabled.}
!           7         Move the IOPT(*) processing pointer.
!           8         Do not preprocess the constraints to
!                     resolve infeasibilities.
!           9         Do not pretriangularize the least squares matrix.
!          99         No more options to change.
!
!    ----
!    X(*)
!    ----
!     This array is used to pass data associated with options 4,5 and
!     6. Ignore this parameter (on input) if no options are used.
!     Otherwise see below: IOPT(*) CONTENTS.
!
!
!    OUTPUT
!    ------
!
!    -----------------
!    X(*),RNORMC,RNORM
!    -----------------
!     The array X(*) contains a solution (if MODE >=0 or  == -22) for
!     the constrained least squares problem. The value RNORMC is the
!     minimum residual vector length for the constraints C*X - Y = 0.
!     The value RNORM is the minimum residual vector length for the
!     least squares equations. Normally RNORMC=0, but in the case of
!     inconsistent constraints this value will be nonzero.
!     The values of X are returned in the first NVARS entries of X(*).
!     The values of Y are returned in the last MCON entries of X(*).
!
!    ----
!    MODE
!    ----
!     The sign of MODE determines whether the subprogram has completed
!     normally, or encountered an error condition or abnormal status. A
!     value of MODE >= 0 signifies that the subprogram has completed
!     normally. The value of mode (>= 0) is the number of variables
!     in an active status: not at a bound nor at the value zero, for
!     the case of free variables. A negative value of MODE will be one
!     of the cases (-57)-(-41), (-37)-(-22), (-19)-(-2). Values  <  -1
!     correspond to an abnormal completion of the subprogram. These
!     error messages are in groups for the subprograms DBOCLS(),
!     DBOLSM(), and DBOLS().  An approximate solution will be returned
!     to the user only when max. iterations is reached, MODE=-22.
!
!    -----------
!    RW(*),IW(*)
!    -----------
!     These are working arrays.  (normally the user can ignore the
!     contents of these arrays.)
!
!    IOPT(*) CONTENTS
!    ------- --------
!     The option array allows a user to modify some internal variables
!     in the subprogram without recompiling the source code. A central
!     goal of the initial software design was to do a good job for most
!     people. Thus the use of options will be restricted to a select
!     group of users. The processing of the option array proceeds as
!     follows: a pointer, here called LP, is initially set to the value
!     1. At the pointer position the option number is extracted and
!     used for locating other information that allows for options to be
!     changed. The portion of the array IOPT(*) that is used for each
!     option is fixed; the user and the subprogram both know how many
!     locations are needed for each option. The value of LP is updated
!     for each option based on the amount of storage in IOPT(*) that is
!     required. A great deal of error checking is done by the
!     subprogram on the contents of the option array. Nevertheless it
!     is still possible to give the subprogram optional input that is
!     meaningless. For example option 4 uses the locations
!     X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing scaling data.
!     The user must manage the allocation of these locations.
!
!   1
!   -
!     This option allows the user to solve problems with a large number
!     of rows compared to the number of variables. The idea is that the
!     subprogram returns to the user (perhaps many times) and receives
!     new least squares equations from the calling program unit.
!     Eventually the user signals "that's all" and a solution is then
!     computed. The value of MROWS is an output variable when this
!     option is used. Its value is always in the range 0 <= MROWS
!     <= NCOLS+1. It is the number of rows after the
!     triangularization of the entire set of equations. If LP is the
!     processing pointer for IOPT(*), the usage for the sequential
!     processing of blocks of equations is
!
!
!        IOPT(LP)=1
!         Move block of equations to W(*,*) starting at
!         the first row of W(*,*).
!        IOPT(LP+3)=# of rows in the block; user defined
!
!     The user now calls DBOCLS( ) in a loop. The value of IOPT(LP+1)