dqed.f90

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

!     This array is used to pass data associated with option 4. Ignore
!     this parameter if this option is not used. Otherwise see below:
!     IOPT(*) CONTENTS.
!
!    OUTPUT
!    ------
!
!    ----------
!    X(*),RNORM
!    ----------
!     The array X(*) contains a solution (if MODE >=0 or  == -22) for
!     the constrained least squares problem. The value RNORM is the
!     minimum residual vector length.
!
!    ----
!    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 (.GE. 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 -37,-36,...,-22, or -17,...,-2. Values  <  -1
!     correspond to an abnormal completion of the subprogram. To
!     understand the abnormal completion codes see below: ERROR
!     MESSAGES for DBOLS( ). AN approximate solution will be returned
!     to the user only when max. iterations is reached, MODE=-22.
!     Values for MODE=-37,...,-22 come from the low-level subprogram
!     DBOLSM(). See the section ERROR MESSAGES for DBOLSM() in the
!     documentation for DBOLSM().
!
!    -----------
!    RW(*),IW(*)
!    -----------
!     These are working arrays with 5*NCOLS and 2*NCOLS entries.
!     (normally the user can ignore the contents of these arrays,
!     but they must be dimensioned properly.)
!
!    IOPT(*) CONTENTS
!    ------- --------
!     The option array allows a user to modify 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. This value is updated as each option is processed. 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. 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 then computes the
!     solution with one final call to subprogram DBOLS( ). The value of
!     MROWS is an OUTPUT variable when this option is used. Its value
!     is always in the range 0 <= MROWS .le. NCOLS+1. It is equal to
!     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 DBOLS( ) in a loop. The value of IOPT(LP+1)
!     directs the user's action. The value of IOPT(LP+2) points to
!     where the subsequent rows are to be placed in W(*,*).
!
!      .<LOOP
!      . CALL DBOLS()
!      . IF(IOPT(LP+1) .EQ. 1) THEN
!      .    IOPT(LP+3)=# OF ROWS IN THE NEW BLOCK; USER DEFINED
!      .    PLACE NEW BLOCK OF IOPT(LP+3) ROWS IN
!      .    W(*,*) STARTING AT ROW IOPT(LP+2).
!      .
!      .    IF( THIS IS THE LAST BLOCK OF EQUATIONS ) THEN
!      .       IOPT(LP+1)=2
!      .<------CYCLE LOOP
!      .    ELSE IF (IOPT(LP+1) .EQ. 2) THEN
!      <-------EXIT LOOP SOLUTION COMPUTED IF MODE .GE. 0
!      . ELSE
!      . ERROR CONDITION; SHOULD NOT HAPPEN.
!      .<END LOOP
!
!     Use of this option adds 4 to the required length of IOPT(*).
!
!
!   2
!   -
!     This option is useful for checking the lengths of all arrays used
!     by DBOLS() against their actual requirements for this problem.
!     The idea is simple: the user's program unit passes the declared
!     dimension information of the arrays. These values are compared
!     against the problem-dependent needs within the subprogram. If any
!     of the dimensions are too small an error message is printed and a
!     negative value of MODE is returned, -11 to -17. The printed error
!     message tells how long the dimension should be. If LP is the
!     processing pointer for IOPT(*),
!
!        IOPT(LP)=2
!        IOPT(LP+1)=Row dimension of W(*,*)
!        IOPT(LP+2)=Col. dimension of W(*,*)
!        IOPT(LP+3)=Dimensions of BL(*),BU(*),IND(*)
!        IOPT(LP+4)=Dimension of X(*)
!        IOPT(LP+5)=Dimension of RW(*)
!        IOPT(LP+6)=Dimension of IW(*)
!        IOPT(LP+7)=Dimension of IOPT(*)
!         .
!        CALL DBOLS()
!
!     Use of this option adds 8 to the required length of IOPT(*).
!
!   3
!   -
!     This option changes the type of scaling for the data matrix E.
!     Nominally each nonzero column of E is scaled so that the
!     magnitude of its largest entry is equal to the value ONE. If LP
!     is the processing pointer for IOPT(*),
!
!        IOPT(LP)=3
!        IOPT(LP+1)=1,2 or 3
!            1= Nominal scaling as noted;
!            2= Each nonzero column scaled to have length ONE;
!            3= Identity scaling; scaling effectively suppressed.
!         .
!        CALL DBOLS()
!
!     Use of this option adds 2 to the required length of IOPT(*).
!
!   4
!   -
!     This option allows the user to provide arbitrary (positive)
!     column scaling for the matrix E. If LP is the processing pointer
!     for IOPT(*),
!
!        IOPT(LP)=4
!        IOPT(LP+1)=IOFF
!        X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1)
!        = Positive scale factors for cols. of E.
!         .
!        CALL DBOLS()
!
!     Use of this option adds 2 to the required length of IOPT(*) and
!     NCOLS to the required length of X(*).
!
!   5
!   -
!     This option allows the user to provide an option array to the
!     low-level subprogram DBOLSM(). If LP is the processing pointer
!     for IOPT(*),
!
!        IOPT(LP)=5
!        IOPT(LP+1)= Position in IOPT(*) where option array
!                    data for DBOLSM() begins.
!         .
!        CALL DBOLS()
!
!     Use of this option adds 2 to the required length of IOPT(*).
!
!   6
!   -
!     Move the processing pointer (either forward or backward) to the
!     location IOPT(LP+1). The processing point is moved to entry
!     LP+2 of IOPT(*) if the option is left with -6 in IOPT(LP).  For
!     example to skip over locations 3,...,NCOLS+2 of IOPT(*),
!
!       IOPT(1)=6
!       IOPT(2)=NCOLS+3
!       (IOPT(I), I=3,...,NCOLS+2 are not defined here.)
!       IOPT(NCOLS+3)=99
!       CALL DBOLS()
!
!     CAUTION: Misuse of this option can yield some very hard
!     -to-find bugs.  Use it with care.
!
!   7
!   -
!     If the user is calling DBOLS() directly, use this option.
!     (This is necessary because DBOCLS() uses DBOLS() as a
!     low-level subprogram.  Due to weighting required within
!     DBOCLS(), the two cases must be known.) For example,
!
!       IOPT(1)=7
!       IOPT(1)=99
!
!   99
!   --
!     There are no more options to change.
!
!     Only option numbers -99, -7,-6,-5,...,-1, 1,2,...,7, and 99 are
!     permitted. Other values are errors. Options -99,-1,...,-7 mean
!     that the repective options 99,1,...,7 are left at their default
!     values. An example is the option to modify the (rank) tolerance:
!
!       IOPT(1)=-3 Option is recognized but not changed
!       IOPT(2)=2  Scale nonzero cols. to have length ONE
!       IOPT(3)=99
!
!***REFERENCES  HANSON, R. J. LINEAR LEAST SQUARES WITH BOUNDS AND
!                 LINEAR CONSTRAINTS, SIAM J. SCI. STAT. COMPUT., VOL. 7,
!                 NO. 3, JULY, 1986.
!***ROUTINES CALLED  IDAMAX,DBOLSM,DCOPY,DNRM2,DROT,DROTG,XERRWV
!***COMMON BLOCKS    DBOCOM
!***END PROLOGUE  DBOLS
!
!     SOLVE LINEAR LEAST SQUARES SYSTEM WITH BOUNDS ON
!     SELECTED VARIABLES.
!
  implicit none
!
  integer mdw
!
  double precision bl(*)
  double precision bu(*)
  logical checkl
  double precision dnrm2
  integer i
  integer ibig
  integer idamax
  integer idope
  integer idum
  integer, save :: igo = 0
  integer ind(*)
  integer inrows
  integer iopt(*)
  integer ip
  integer iscale
  integer iw(*)
  integer j
  integer jp
  integer lds
  integer lenx
  integer level
  integer liopt
  integer llb
  integer lliw
  integer llrw
  integer llx
  integer lmdw
  integer lndw
  integer locacc
  integer locdim
  integer lopt
  integer lp
  integer mnew
  integer mode
  integer mrows
  integer ncols
  integer nerr
  double precision, parameter :: one = 1.0D+00
  double precision rdum
  double precision rnorm
  double precision rw(*)
  double precision sc
  double precision ss
  double precision w(mdw,*)
  double precision x(*)
!
  common /dbocom/ idope(5)
!
  save iscale
  save locacc
  save lopt
!
  idum = 0
  rdum = 0.0D+00
  nerr = 0
  mode = 0
  level = 1
  if ( igo == 0 ) then
!
!  CHECK VALIDITY OF INPUT DATA
!
!  SEE THAT MDW IS .GT.0. GROSS CHECK ONLY.
!
      if ( mdw<=0) then
          nerr = 2
          call xerrwv('dbols(). mdw=(i1) must be positive.', &
            nerr,level,1,mdw,idum,0,rdum,rdum)
          go to 190
      end if
!
!  SEE THAT NUMBER OF UNKNOWNS IS POSITIVE.
!
      if ( ncols<=0) then
          nerr = 3
          call xerrwv('dbols(). ncols=(i1) the no. of variables must ' // &
            'be positive.',nerr,level,1,ncols,idum,0,rdum,rdum)
          go to 190
      end if
!
!  SEE THAT CONSTRAINT INDICATORS ARE ALL WELL-DEFINED.
!
      do 10 j = 1,ncols
          if ( ind(j) < 1 .or. ind(j) > 4) then
              nerr = 4
              call xerrwv( &
                     'dbols(). for j=(i1), ind(j)=(i2) must be 1-4.' &
                          ,nerr,level,2,j,ind(j),0,rdum,rdum)
              go to 190
          end if
   10     continue
!
!  SEE THAT BOUNDS ARE CONSISTENT.
!
      do j = 1,ncols
          if ( ind(j) == 3) then
              if ( bl(j) > bu(j)) then
                  nerr = 5
                  call xerrwv( &
         'dbols(). for j=(i1), bound bl(j)=(r1) is  >  bu(j)=(r2).' &
                   ,nerr,level,1,j,idum,2,bl(j),bu(j))
                  go to 190
              end if
          end if
      end do
!
!  PROCESS OPTION ARRAY
!
      checkl = .false.
      lenx = ncols
      iscale = idope(4)
      igo = 2
      lopt = 0
      lp = 0
      lds = 0

30    continue

      lp = lp + lds
      ip = iopt(lp+1)
      jp = abs(ip)
!
!  TEST FOR NO MORE OPTIONS.
!
      if ( ip == 99) then
          if ( lopt == 0) lopt = lp + 1
          go to 50
      else if ( jp == 99) then
          lds = 1
          go to 30
      else if ( jp == 1) then
          if ( ip > 0) then
!
!  SET UP DIRECTION FLAG, ROW STACKING POINTER
!  LOCATION, AND LOCATION FOR NUMBER OF NEW ROWS.
!
              locacc = lp + 2
!
!                  IOPT(LOCACC-1)=OPTION NUMBER FOR SEQ. ACCUMULATION.
!     CONTENTS..   IOPT(LOCACC  )=USER DIRECTION FLAG, 1 OR 2.
!                  IOPT(LOCACC+1)=ROW STACKING POINTER.
!                  IOPT(LOCACC+2)=NUMBER OF NEW ROWS TO PROCESS.
!     USER ACTION WITH THIS OPTION..
!      (SET UP OPTION DATA FOR SEQ. ACCUMULATION IN IOPT(*).
!      MUST ALSO START PROCESS WITH IOPT(LOCACC)=1.)
!      (MOVE BLOCK OF EQUATIONS INTO W(*,*)  STARTING AT FIRST
!       ROW OF W(*,*).  SET IOPT(LOCACC+2)=NO. OF ROWS IN BLOCK.)
!              LOOP
!              CALL DBOLS()
!
!                  IF(IOPT(LOCACC) .EQ. 1) THEN
!                      STACK EQUAS., STARTING AT ROW IOPT(LOCACC+1),
!                       INTO W(*,*).
!                       SET IOPT(LOCACC+2)=NO. OF EQUAS.
!                      IF LAST BLOCK OF EQUAS., SET IOPT(LOCACC)=2.
!                  ELSE IF IOPT(LOCACC) .EQ. 2) THEN
!                      (PROCESS IS OVER. EXIT LOOP.)
!                  ELSE
!                      (ERROR CONDITION. SHOULD NOT HAPPEN.)
!                  END IF
!              END LOOP
!              SET IOPT(LOCACC-1)=-OPTION NUMBER FOR SEQ. ACCUMULATION.
!              CALL DBOLS( )
              iopt(locacc+1) = 1
              igo = 1
          end if
          lds = 4
          go to 30
      else if ( jp == 2) then
          if ( ip > 0) then
!
!  GET ACTUAL LENGTHS OF ARRAYS FOR CHECKING AGAINST NEEDS.
!
              locdim = lp + 2
!
!  LMDW.GE.MROWS
!  LNDW.GE.NCOLS+1
!  LLB .GE.NCOLS
!  LLX .GE.NCOLS+EXTRA REQD. IN OPTIONS.
!  LLRW.GE.5*NCOLS
!  LLIW.GE.2*NCOLS
!  LIOP.GE. AMOUNT REQD. FOR IOPTION ARRAY.
!
              lmdw = iopt(locdim)
              lndw = iopt(locdim+1)
              llb = iopt(locdim+2)
              llx = iopt(locdim+3)
              llrw = iopt(locdim+4)
              lliw = iopt(locdim+5)
              liopt = iopt(locdim+6)
              checkl = .true.
          end if
          lds = 8
          go to 30
!
!  OPTION TO MODIFY THE COLUMN SCALING.
!
      else if ( jp == 3) then
          if ( ip > 0) then
              iscale = iopt(lp+2)
!
!  SEE THAT ISCALE IS 1 THRU 3.
!
              if ( iscale < 1 .or. iscale > 3) then
                  nerr = 7
                  call xerrwv( 'dbols(). iscale option=(i1) must be 1-3.' &
                    ,nerr,level,1,iscale,idum,0,rdum,rdum)
                  go to 190
              end if
          end if
          lds = 2
          go to 30
!
!  IN THIS OPTION THE USER HAS PROVIDED SCALING.  THE
!  SCALE FACTORS FOR THE COLUMNS BEGIN IN X(NCOLS+IOPT(LP+2)).
!
      else if ( jp == 4) then
          if ( ip > 0) then
              iscale = 4
              if ( iopt(lp+2)<=0) then
                  nerr = 8
                  call xerrwv('dbols(). offset past x(ncols) (i1) ' // &
                    'for user-provided column scaling must be positive.', &
                    nerr,level,1,iopt(lp+2),idum,0,rdum,rdum)
                  go to 190
              end if
              call dcopy(ncols,x(ncols+iopt(lp+2)),1,rw,1)
              lenx = lenx + ncols
              do 40 j = 1,ncols
                  if ( rw(j)<= 0.