dqed.f90

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

!     directs the user's action. The value of IOPT(LP+2) points to
!     where the subsequent rows are to be placed in W(*,*). Both of
!     these values are first defined in the subprogram. The user
!     changes the value of IOPT(LP+1) (to 2) as a signal that all of
!     the rows have been processed.
!
!
!      .<LOOP
!      . CALL DBOCLS( )
!      . 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 MCON + 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 DBOCLS( ) 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, -41 to -47. 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 DBOCLS( )
!
!     Use of this option adds 8 to the required length of IOPT(*).
!
!   3
!   -
!     This option can change the type of scaling for the data matrix.
!     Nominally each nonzero column of the matrix 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 DBOCLS( )
!
!     Use of this option adds 2 to the required length of IOPT(*).
!
!   4
!   -
!     This options allows the user to provide arbitrary (positive)
!     column scaling for the matrix. 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 DBOCLS( )
!
!     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 DBOLS( ). If LP is the processing pointer
!     for IOPT(*),
!
!        IOPT(LP)=5
!        IOPT(LP+1)= Position in IOPT(*) where option array
!                    data for DBOLS( ) begins.
!         .
!        CALL DBOCLS( )
!
!     Use of this option adds 2 to the required length of IOPT(*).
!
!   6
!   -
!     This option is no longer operative.  To pass an option array
!     to the low-level subprogram DBOLSM( ), imbed it within an option
!     array passed to DBOLS() using option 5.
!
!   7
!   -
!     Move the processing pointer (either forward or backward) to the
!     location IOPT(LP+1). The processing pointer moves to locations
!     LP+2 if option number 7 is used with the value -7.  For
!     example to skip over locations 3,...,NCOLS+2,
!
!       IOPT(1)=7
!       IOPT(2)=NCOLS+3
!       (IOPT(I), I=3,...,NCOLS+2 are not defined here.)
!       IOPT(NCOLS+3)=99
!       CALL DBOCLS( )
!
!     CAUTION: Misuse of this option can yield some very hard-to-find
!     bugs. Use it with care. It is intended to be used for passing
!     option arrays to other subprograms.
!
!   8
!   -
!     This option allows the user to suppress the algorithmic feature
!     of DBOCLS( ) that processes the constraint equations C*X = Y and
!     resolves infeasibilities. The steps normally done are to solve
!     C*X - Y = 0 in a least squares sense using the stated bounds on
!     both X and Y. Then the "reachable" vector Y = C*X is computed
!     using the solution X obtained. Finally the stated bounds for Y are
!     enlarged to include C*X. To suppress the feature:
!
!
!       IOPT(LP)=8
!         .
!       CALL DBOCLS( )
!
!     Use of this option adds 1 to the required length of IOPT(*).
!
!   9
!   -
!     This option allows the user to suppress the pretriangularizing
!     step of the least squares matrix that is done within DBOCLS( ).
!     This is primarily a means of enhancing the subprogram efficiency
!     and has little effect on accuracy. To suppress the step, set:
!
!       IOPT(LP)=9
!         .
!       CALL DBOCLS( )
!
!     Use of this option adds 1 to the required length of IOPT(*).
!
!   99
!   --
!     There are no more options to change.
!
!     Only option numbers -99, -9,-8,...,-1, 1,2,...,9, and 99 are
!     permitted. Other values are errors. Options -99,-1,...,-9 mean
!     that the respective options 99,1,...,9 are left at their default
!     values. An example is the option to suppress the preprocessing of
!     contraints:
!
!       IOPT(1)=-8 Option is recognized but not changed
!       IOPT(2)=99
!       CALL DBOCLS( )
!
!***END PROLOGUE  DBOCLS
!     REVISED 880722-1100
!     REVISED YYMMDD-HHMM
!
!    PURPOSE
!    -------
!     THIS IS THE MAIN SUBPROGRAM THAT SOLVES THE LEAST SQUARES
!     PROBLEM CONSISTING OF LINEAR CONSTRAINTS
!
!              C*X = Y
!
!     AND LEAST SQUARES EQUATIONS
!
!              E*X = F
!
!     IN THIS FORMULATION THE VECTORS X AND Y ARE BOTH UNKNOWNS.
!     FURTHER, X AND Y MAY BOTH HAVE USER-SPECIFIED BOUNDS ON EACH
!     COMPONENT.  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(16+NI), IW(2*(NCOLS+MCON))
!
  implicit none
!
  integer mdw
!
  logical, save :: accum
  double precision anorm
  double precision bl(*)
  double precision bu(*)
  logical, save :: checkl
  double precision cnorm
  double precision dasum
  double precision ddot
  double precision dnrm2
  double precision drelpr
  logical filter
  integer i
  integer icase
  integer idope
  integer idum
  integer, save :: igo = 0
  integer iiw
  integer ind(*)
  integer inrows
  integer iopt(*)
  integer ip
  integer irw
  integer iscale
  integer iw(*)
  integer j
  integer jopt(5)
  integer jp
  integer lds
  integer lenx
  integer level
  integer liopt
  integer liw
  integer llb
  integer lliw
  integer llrw
  integer llx
  integer lmdw
  integer lndw
  integer locacc
  integer locdim
  integer lopt
  integer lp
  integer lrw
  integer mcon
  integer mdwl
  integer mnew
  integer mode
  integer modec
  integer mout
  integer mrows
  integer ncols
  integer nerr
  double precision, parameter :: one = 1.0D+00
  logical pretri
  double precision rdum
  double precision rnorm
  double precision rnormc
  double precision rw(*)
  double precision t
  double precision t1
  double precision t2
  double precision w(mdw,*)
  double precision wt
  double precision x(*)
!
  common /dbocom/ idope(5)
!
  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 = 53
          call xerrwv('dbocls(). mdw=(i1) must be positive.', &
                      nerr,level,1,mdw,idum,0,rdum,rdum)
          go to 260
      end if
!
!  SEE THAT NUMBER OF CONSTRAINTS IS NONNEGATIVE.
!
      if ( mcon < 0) then
          nerr = 54
          call xerrwv('dbocls(). mcon=(i1) must be nonnegative.', &
                      nerr,level,1,mcon,idum,0,rdum,rdum)
          go to 260
      end if
!
!  SEE THAT NUMBER OF UNKNOWNS IS POSITIVE.
!
      if ( ncols<=0) then
          nerr = 55
          call xerrwv( &
       'dbocls(). ncols=(i1) the no. of variables must be positive.' &
                      ,nerr,level,1,ncols,idum,0,rdum,rdum)
          go to 260
      end if
!
!  SEE THAT CONSTRAINT INDICATORS ARE ALL WELL-DEFINED.
!
      do j = 1,ncols + mcon
          if ( ind(j) < 1 .or. ind(j) > 4) then
              nerr = 56
              call xerrwv( &
                    'dbocls(). for j=(i1), ind(j)=(i2) must be 1-4.' &
                          ,nerr,level,2,j,ind(j),0,rdum,rdum)
              go to 260
          end if
      end do
!
!  SEE THAT BOUNDS ARE CONSISTENT.
!
      do j = 1,ncols + mcon
          if ( ind(j) == 3) then
              if ( bl(j) > bu(j)) then
                  nerr = 57
                  call xerrwv( &
        'dbocls(). for j=(i1), bound bl(j)=(r1) is  >  bu(j)=(r2).' &
         ,nerr,level,1,j,idum,2,bl(j),bu(j))
                  go to 260
              end if
          end if
      end do
!
!  PROCESS OPTION ARRAY
!
      drelpr = epsilon ( drelpr )
      checkl = .false.
      filter = .true.
      lenx = 2* (ncols+mcon) + 2
      iscale = 1
      igo = 1
      accum = .false.
      pretri = .true.
      lopt = 0
      lp = 0
      lds = 0

   30     continue
      lp = lp + lds
      ip = iopt(lp+1)
      jp = abs(ip)
!
!  TEST FOR NO MORE OPTIONS TO CHANGE.
!
      if ( ip == 99) then
          if ( lopt == 0) lopt = lp+1
!
!  SEND COL. SCALING TO DBOLS().
!
          idope(4)=1
!
!  NOTE THAT DBOLS() WAS CALLED BY DBOCLS()
!
          idope(5)=1
!
!  CHANGE PRETRIANGULARIZATION FACTOR IN DBOLSM().
!
          idope(1) = ncols + mcon + 1
!
!  PASS WEIGHT TO DBOLSM() FOR RANK TEST.
!
          idope(2) = ncols + mcon + 2
          idope(3) = mcon
          go to 50
      else if ( jp == 99) then
          lds = 1
          go to 50
      else if ( jp == 1) then
          if ( ip > 0) then
!
!  SET UP DIRECTION FLAG LOCATION, 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(*).)
!      (MOVE BLOCK OF EQUATIONS INTO W(*,*)  STARTING AT FIRST
!       ROW OF W(*,*) BELOW THE ROWS FOR THE CONSTRAINT MATRIX C.
!       SET IOPT(LOCACC+2)=NO. OF LEAST SQUARES EQUATIONS IN BLOCK.
!
!              LOOP
!              CALL DBOCLS()
!
!                  IF(IOPT(LOCACC) .EQ. 1) THEN
!                      STACK EQUAS. INTO W(*,*), 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
!
              iopt(locacc+1) = mcon + 1
              accum = .true.
              iopt(locacc) = igo
          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.MCON+MAX(MOUT,NCOLS), IF MCON.GT.0 .AND FILTER
!  LMDW.GE.MCON+MOUT, OTHERWISE
!
!  LNDW.GE.NCOLS+MCON+1
!  LLB .GE.NCOLS+MCON
!  LLX .GE.2*(NCOLS+MCON)+2+EXTRA REQD. IN OPTIONS.
!  LLRW.GE.6*NCOLS+5*MCON
!  LLIW.GE.2*(NCOLS+MCON)
!  LIOP.GE. AMOUNT REQD. FOR OPTION 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 = 48
                  call xerrwv('dbocls(). iscale option=(i1) must be 1-3.' &
                    ,nerr,level,1,iscale,idum,0,rdum,rdum)
                  go to 260
              end if
          end if
          lds = 2
          go to 30
!
!  IN THIS OPTION THE USER HAS PROVIDED SCALING.  THE