dqed.f90

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

!  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 = 49
                  call xerrwv('dbocls(). 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 260
              end if
              call dcopy(ncols,x(ncols+iopt(lp+2)),1,rw,1)
              lenx = lenx + ncols
              do j = 1,ncols
                  if ( rw(j)<=0.0D+00 ) then
                      nerr = 50
                      call xerrwv('dbocls(). each provided col. scale ' // &
                        'factor must be positive. comp. (i1)   now = (r1).', &
                        nerr,level,1,j,idum,1,rw(j),rdum)
                      go to 260
                  end if
              end do
          end if
          lds = 2
          go to 30
!
!  IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO DBOLS().
!
      else if ( jp == 5) then
          if ( ip > 0) then
              lopt = iopt(lp+2)
          end if
          lds = 2
          go to 30
!
!  IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO DBOLSM().
!  (NO LONGER USED.) OPTION NOW MUST BE PASSED IMBEDDED IN
!  OPTION ARRAY FOR DBOLS().
!
      else if ( jp == 6) then
          lds = 2
          go to 30
!
!  THIS OPTION USES THE NEXT LOC OF IOPT(*) AS A
!  POINTER VALUE TO SKIP TO NEXT.
!
      else if ( jp == 7) then
          if ( ip > 0) then
              lp = iopt(lp+2)-1
              lds = 0
          else
              lds = 2
          end if
          go to 30
!
!  THIS OPTION AVOIDS THE CONSTRAINT RESOLVING PHASE FOR
!  THE LINEAR CONSTRAINTS C*X=Y.
!
      else if ( jp == 8) then
          filter = .not. (ip > 0)
          lds = 1
          go to 30
!
!  THIS OPTION SUPPRESSES PRETRIANGULARIZATION OF THE LEAST
!  SQUARES EQATIONS.
!
      else if ( jp == 9) then
          pretri = .not. (ip > 0)
          lds = 1
          go to 30
!
!  NO VALID OPTION NUMBER WAS NOTED. THIS IS AN ERROR CONDITION.
!
      else
          nerr = 51
          call xerrwv('dbocls(). the option number=(i1) is not defined.' &
            ,nerr,level,1,jp,idum,0,rdum,rdum)
          go to 260
      end if

   50     continue

      if ( checkl) then
!
!  CHECK LENGTHS OF ARRAYS
!
!  THIS FEATURE ALLOWS THE USER TO MAKE SURE THAT THE
!  ARRAYS ARE LONG ENOUGH FOR THE INTENDED PROBLEM SIZE AND USE.
!
       if ( filter .and. .not.accum) then
         mdwl=mcon+max(mrows,ncols)
       else if ( accum) then
         mdwl=mcon+ncols+1
       else
         mdwl=mcon+ncols
       end if

          if ( lmdw < mdwl) then
              nerr = 41
              call xerrwv('dbocls(). the row dimension of w(,)=(i1) ' // &
                'must be >= the number of effective rows=(i2).', &
                nerr,level,2,lmdw,mdwl,0,rdum,rdum)
              go to 260
          end if
          if ( lndw < ncols+mcon+1) then
              nerr = 42
              call xerrwv('dbocls(). the column dimension of w(,)=(i1) ' // &
                'must be >= ncols+mcon+1=(i2).',nerr,level,2,lndw, &
                ncols+mcon+1,0,rdum,rdum)
              go to 260
          end if
          if ( llb < ncols+mcon) then
              nerr = 43
              call xerrwv('dbocls(). the dimensions of the arrays bl(),' // &
                'bu(), and ind()=(i1) must be >= ncols+mcon=(i2).', &
                nerr,level,2,llb,ncols+mcon,0,rdum,rdum)
              go to 260
          end if

          if ( llx < lenx) then
              nerr = 44
              call xerrwv( 'dbocls(). the dimension of x()=(i1) must be ' // &
                '>= the reqd. length=(i2).',nerr,level,2,llx,lenx, &
                0,rdum,rdum)
              go to 260
          end if

          if ( llrw < 6*ncols+5*mcon) then
              nerr = 45
              call xerrwv('dbocls(). the dimension of rw()=(i1) must be ' // &
                '>= 6*ncols+5*mcon=(i2).',nerr,level,2, &
                llrw,6*ncols+5*mcon,0,rdum,rdum)
              go to 260
          end if

          if ( lliw < 2*ncols+2*mcon) then
              nerr = 46
              call xerrwv('dbocls() the dimension of iw()=(i1) must be ' // &
                '>= 2*ncols+2*mcon=(i2).',nerr,level,2,lliw, &
                2*ncols+2*mcon,0,rdum,rdum)
              go to 260
          end if

          if ( liopt < lp+17) then
              nerr = 47
              call xerrwv('dbocls(). the dimension of iopt()=(i1) must be ' // &
                '>= the reqd. len.=(i2).',nerr,level,2,liopt,lp+17, 0,rdum,rdum)
              go to 260
          end if
      end if
  end if
!
!  OPTIONALLY GO BACK TO THE USER FOR ACCUMULATION OF LEAST SQUARES
!  EQUATIONS AND DIRECTIONS FOR PROCESSING THESE EQUATIONS.
!
!  ACCUMULATE LEAST SQUARES EQUATIONS
!
  if ( accum) then
      mrows = iopt(locacc+1) - 1 - mcon
      inrows = iopt(locacc+2)
      mnew = mrows + inrows
      if ( mnew < 0 .or. mnew+mcon > mdw) then
          nerr = 52
          call xerrwv('dbocls(). no. of rows=(i1) must be >= 0 ' // &
            '.and. <=mdw-mcon=(i2)',nerr,level,2,mnew,mdw-mcon,0,rdum,rdum)
          go to 260
      end if
  end if
!
!  USE THE SOFTWARE OF DBOLS( ) FOR THE TRIANGULARIZATION OF THE
!  LEAST SQUARES MATRIX.  THIS MAY INVOLVE A SYSTALTIC INTERCHANGE
!  OF PROCESSING POINTERS BETWEEN THE CALLING AND CALLED (DBOLS())
!  PROGRAM UNITS.
!
  jopt(01) = 1
  jopt(02) = 2
  jopt(04) = mrows
  jopt(05) = 99
  irw = ncols + 1
  iiw = 1
  if ( accum .or. pretri) then
!
!  NOTE THAT DBOLS() WAS CALLED BY DBOCLS()
!
          idope(5)=0
      call dbols(w(mcon+1,1),mdw,mout,ncols,bl,bu,ind,jopt,x,rnorm, &
             mode,rw(irw),iw(iiw))
  else
      mout = mrows
  end if

  if ( accum) then
    accum = iopt(locacc)  ==  1
    iopt(locacc+1) = jopt(03) + mcon
    mrows = min(ncols+1,mnew)
  end if

  if ( accum) return
!
!  SOLVE CONSTRAINED AND BOUNDED LEAST SQUARES PROBLEM
!
!  MOVE RIGHT HAND SIDE OF LEAST SQUARES EQUATIONS.
!
  call dcopy(mout,w(mcon+1,ncols+1),1,w(mcon+1,ncols+mcon+1),1)
  if ( mcon > 0 .and. filter) then
!
!  PROJECT THE LINEAR CONSTRAINTS INTO A REACHABLE SET.
!
      do i = 1,mcon
        call dcopy(ncols,w(i,1),mdw,w(mcon+1,ncols+i),1)
      end do
!
!  PLACE (-)IDENTITY MATRIX AFTER CONSTRAINT DATA.
!
      do j = ncols + 1,ncols + mcon + 1
        w(1:mcon,j) = 0.0D+00
      end do
      do j = ncols + 1,ncols + mcon + 1
        w(j-ncols,j) = -1.0D+00
      end do
!
!  OBTAIN A 'FEASIBLE POINT' FOR THE LINEAR CONSTRAINTS.
!
      jopt(01) = 99
      irw = ncols + 1
      iiw = 1
!
!  NOTE THAT DBOLS() WAS CALLED BY DBOCLS()
!
          idope(5)=0
      call dbols(w,mdw,mcon,ncols+mcon,bl,bu,ind,jopt,x,rnormc, &
        modec,rw(irw),iw(iiw))
!
!  ENLARGE THE BOUNDS SET, IF REQUIRED, TO INCLUDE POINTS THAT
!  CAN BE REACHED.
!
      do j = ncols + 1,ncols + mcon
          icase = ind(j)
          if ( icase < 4) then
              t = ddot ( ncols, w(mcon+1,j), 1, x, 1 )
          end if
          go to (80,90,100,110),icase
          go to 120
   80         bl(j) = min(t,bl(j))
          go to 120
   90         bu(j) = max(t,bu(j))
          go to 120
  100         bl(j) = min(t,bl(j))
          bu(j) = max(t,bu(j))
          go to 120
  110         continue
  120         continue
      end do
!
!  MOVE CONSTRAINT DATA BACK TO THE ORIGINAL AREA.
!
      do j = ncols + 1,ncols + mcon
          call dcopy(ncols,w(mcon+1,j),1,w(j-ncols,1),mdw)
      end do

  end if

  if ( mcon > 0) then
      do j = ncols + 1,ncols + mcon
        w(mcon+1:mcon+mout,j) = 0.0D+00
      end do
!
!  PUT IN (-)IDENTITY MATRIX (POSSIBLY) ONCE AGAIN.
!
      do j = ncols + 1,ncols + mcon + 1
        w(1:mcon,j) = 0.0D+00
      end do

      do i = 1, mcon
        w(i,ncols+i) = - 1.0D+00
      end do

  end if
!
!  COMPUTE NOMINAL COLUMN SCALING FOR THE UNWEIGHTED MATRIX.
!
  cnorm = 0.0D+00
  anorm = 0.0D+00
  do j = 1,ncols
      t1 = dasum(mcon,w(1,j),1)
      t2 = dasum(mout,w(mcon+1,1),1)
      t = t1 + t2
      if ( t == 0.0D+00 ) t = one
      cnorm = max(cnorm,t1)
      anorm = max(anorm,t2)
      x(ncols+mcon+j) = one/t
  end do

  go to (180,190,210,220),iscale
  go to 230
  180 continue
  go to 230
!
!  SCALE COLS. (BEFORE WEIGHTING) TO HAVE LENGTH ONE.
!
  190 continue

  do j = 1,ncols
    t = dnrm2(mcon+mout,w(1,j),1)
    if ( t == 0.0D+00 ) t = one
    x(ncols+mcon+j) = one/t
  end do

  go to 230
!
!  SUPPRESS SCALING (USE UNIT MATRIX).
!
  210 continue

  x(ncols+mcon+1) = one
  call dcopy(ncols,x(ncols+mcon+1),0,x(ncols+mcon+1),1)
  go to 230
!
!  THE USER HAS PROVIDED SCALING.
!
  220 call dcopy(ncols,rw,1,x(ncols+mcon+1),1)
  230 continue

  do j = ncols + 1,ncols + mcon
    x(ncols+mcon+j) = one
  end do
!
!  WEIGHT THE LEAST SQUARES EQUATIONS.
!
  wt = sqrt(drelpr)
  if ( anorm > 0.0D+00 ) wt = wt/anorm
  if ( cnorm > 0.0D+00 ) wt = wt*cnorm

  do i = 1,mout
      call dscal(ncols,wt,w(i+mcon,1),mdw)
  end do
  call dscal(mout,wt,w(mcon+1,mcon+ncols+1),1)
  lrw = 1
  liw = 1
!
!  SET THE NEW TRIANGULARIZATION FACTOR.
!
  x(ncols+mcon+idope(1))= 0.0D+00
!
!  SET THE WEIGHT TO USE IN COMPONENTS .GT. MCON,
!  WHEN MAKING LINEAR INDEPENDENCE TEST.
!
  x(ncols+mcon+idope(2))= one/wt
  idope(5)=1
  call dbols(w,mdw,mout+mcon,ncols+mcon,bl,bu,ind,iopt(lopt),x, &
    rnorm,mode,rw(lrw),iw(liw))

  260 if ( mode>=0)mode=-nerr
  igo = 0
  return
end
subroutine dbols ( w, mdw, mrows, ncols, bl, bu, ind, iopt, x, rnorm, &
  mode, rw, iw )
!
!***********************************************************************
!
!! DBOLS solves the linear system E*X = F in the least squares sense.
!
!
!  Discussion:
!
!    This routine solves the problem
!
!      E*X = F 
!
!    in the least squares sense with bounds on selected X values.
!
!  Author:
!
!    Richard Hanson, Sandia National Laboratory
!
!***DESCRIPTION
!
!     The user must have dimension statements of the form:
!
!       DIMENSION W(MDW,NCOLS+1), BL(NCOLS), BU(NCOLS),
!      * X(NCOLS+NX), RW(5*NCOLS)
!       integer IND(NCOLS), IOPT(1+NI), IW(2*NCOLS)
!
!     (here NX=number of extra locations required for option 4; NX=0
!     for no options; NX=NCOLS if this option is in use. Here NI=number
!     of extra locations required for options 1-6; NI=0 for no
!     options.)
!
!   INPUT
!   -----
!
!    --------------------
!    W(MDW,*),MROWS,NCOLS
!    --------------------
!     The array W(*,*) contains the matrix [E:F] on entry. The matrix
!     [E:F] has MROWS rows and NCOLS+1 columns. This data is placed in
!     the array W(*,*) with E occupying the first NCOLS columns and the
!     right side vector F in column NCOLS+1. The row dimension, MDW, of
!     the array W(*,*) must satisfy the inequality MDW >= MROWS.
!     Other values of MDW are errrors. The values of MROWS and NCOLS
!     must be positive. Other values are errors. There is an exception
!     to this 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.  MROWS contains the number of rows in the matrix after
!     triangularizing several blocks of equations. This is an OUTPUT
!     parameter ONLY when option 1 is used. 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.
!
!    1.    For IND(J)=1, require X(J) >= BL(J).
!          (the value of BU(J) is not used.)
!    2.    For IND(J)=2, require X(J) <= 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).
!    4.    For IND(J)=4, no bounds on X(J) 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.
!
!    -------
!    IOPT(*)
!    -------
!     This is the array where the user can specify nonstandard options
!     for DBOLSM( ). 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.
!           2         Check lengths of all arrays used in the
!                     subprogram.
!           3         Standard scaling of the data matrix, E.
!           4         User provides column scaling for matrix E.
!           5         Provide option array to the low-level
!                     subprogram DBOLSM( ).
!           6         Move the IOPT(*) processing pointer.
!           7         User has called DBOLS() directly.
!          99         No more options to change.
!
!    ----
!    X(*)
!    ----