program sbl1
c
c
c  Strict bounds on Linear Functionals of Bounded Variables Subject
c     to a One-Norm Misfit Criterion on Linear Data Constraints
c
c  P.B. Stark      Version 12/29/92
c
c  Finds bounds on linear functionals of a model required to satisfy
c  a priori upper and lower bounds and a set of linear equations within
c  an l-1 norm error tolerance:
c
c
c  min pv.z   such that             bl<=z<=bu and 
c  max                    sum_i |A(i).z-dbar(i)|/dsv(i) <= chi*(1+tol) 
c
c  for a list of vectors pv.
c
c  Quadratic programming with heavy weights on the data constraints
c  is used to simulate linear programming with inequality constraints.
c  The use of quadratic programming thus necessitates the introduction
c  of copious slack variables, which subroutine  prep  handles.
c
c* Uses BVLS (Bounded-Variable Least-Squares) quadratic programming
c* routine by R.L. Parker and P.B. Stark, 1987.
c*
c*
c* Calls  BVLS, CHECK, DOT, FORWRD, GETBND, GETMAP, GETVEC, 
c*        PREP, PREPV, SETDAT, SLACKS, and SPACE.
c
c*********************************************************************
c******************************WARNING********************************
c
c  The program contains two ``tuning constants:''   scale  and  tol.
c  They are set in the data statement to 1.d-2 .
c  These values have worked in a wide variety of test cases.
c  ``scale'' is the relative weight of the (bottom) row of the matrix  
c  aa sent to bvls.  It specifies the tradeoff between minimizing the 
c  difference between pv.z and its target value, and fitting the data.
c  tol specifies the relative tolerance with which the target value
c  of the misfit,  chi,  should be achieved.
c
c*********************************************************************
c
c  Input required:
c
c     print mode integer (1-4)   gives increasing amounts of output.
c
c     data file name             file should contain pairs:
c                                  central value   standard deviation
c
c     number of data
c
c     number of levels in the
c       model
c
c     number of objective
c       vectors to find
c       bounds for
c
c     file name for a priori     file should have pairs:
c       upper and lower bounds     bl(i), bu(i)
c
c     file name for              file should contain the matrix,
c       data mapping matrix        row by row (natural order)
c
c     file name for objective    file should contain the vectors,
c       vectors                    one after another
c
c
      implicit double precision (a-h,o-z)
      parameter (nmax=15000, nimax=500 )
c  nmax:   maximum allocated array storage.
c  nimax:  maximum integer array storage.
      common /c/ c(nmax)
      common /s/ scale, pmin, pmax
c  scale:  relative weight of constraints versus penalty functional.
c          bounds on pv.z
c  pmin:   a priori lower bound on pv.z
c  pmax:   a priori upper bound on pv.z

      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
c  dbar:   data central values.
c  dsv:    data standard deviations.
c  chi:    maximum acceptable l1-norm of misfit.  Should be chosen
c          to provide a 1-alpha confidence region for the noise-free
c          data values.  Percentage points for the 1-norm of iid
c          Gaussian variables are tabulated in Parker and McNutt, JGR, 
c          v85, p4429, 1980.
c  nd:     number of data.
c  lvls:   number of elements in the model.
c  nv:     number penalty vectors for which to find extremal bounds.
c
      common /icons/ nz,nbnd,map,npv,naa,nbb,nbl,nbu,
     +  nw,nzz,nact,ntot
c
c  nz:    starting location of the current solution in /c/.
c  nbnd:  starting location of the tabulated bounds.  First nv elements
c         of /c/ from nbnd are lower bounds; next nv elements are upper.
c  map:   starting location of the data mapping matrix.
c  npv:   starting location of the penalty vector.
c  naa:   starting location of matrix aa for bvls.
c  nbb:   starting location of vector bb for bvls.
c  nbl:   starting location of lower bounds on the model.
c  nbu:   starting location of upper bounds on the model.
c  nw:    starting location for bvls working array w.
c  nzz:   starting location for bvls working array zz.
c  nact:  starting location for bvls working array act.
c  ntot:  total required storage for the above arrays.
      common /print/mode,id
c  mode:  determines the amount of printed output.
c  id:    disk fortran unit number.
      dimension istate (nimax)
c  istate:integer working array for bvls.
      data id/7/, scale /1.d-2/, tol/1.d-2/
c  tol:   the fitting tolerance for the data (a fraction of chi)
c
c------------------------first executable statement--------------------
c
 100  print '(/a/a)','Enter print mode integer (1-4): ',
     +  'larger values yield more printing. Enter -1 to quit.'
      read *, mode
      if ( mode .lt. 1 ) stop
      call setdat
c  Partition the main common block.
      call space
      iused=lvls+2*nd+4
      print '(/)'
      print *,'Real array space required ',ntot,' reserved ',nmax
      print *,'Integer space required ',iused,' reserved ',nimax
      if(ntot.gt.nmax.or.iused.gt.nimax) go to 100
      print '(/)'
c  Get a priori upper/lower bounds on model.
      call getbnd(lvls,c(nbl),c(nbu))
c  Set the bounds on the slack variables using the data tolerances.
      call slacks(lvls+2*nd+3, c(nbl),c(nbu))
c  Get the data mapping matrix.
      call getmap(c(map),nd,lvls)
c
c  Test for problem feasibility.
c  Set up the matrices and bounds without penalty vector.
      call prep(nd,lvls,nd+1,lvls+2*nd+3,c(map),c(naa),c(nbb))
c Set the bounds on the slacks for the objective vector.
      c(nbl+lvls+2*nd+1)=0.0
      c(nbl+lvls+2*nd+2)=0.0
      c(nbu+lvls+2*nd+1)=0.0
      c(nbu+lvls+2*nd+2)=0.0
      key=0
      call bvls(key,nd+1,lvls+2*nd+3,c(naa),c(nbb),c(nbl),c(nbu),
     + c(nz),c(nw),c(nact),c(nzz),istate,its)
      call check(c(nz),lvls,c(map),nd,viol)
      if (mode.ge.2) print *,' l-1 misfit 1st trial ', viol
      if (mode.ge.3) call forwrd(c(nz),lvls,c(map),nd)
      if (viol .ge. chi*(1.0+tol) ) then
        print '(/)'
        print *,'**********************************************'
        print *,' DATA AND A PRIORI CONSTRAINTS INCONSISTENT.'
        print *,' ONE-NORM MISFIT ', viol, ' ALLOWED ',chi
        print *,' Check data, mapping matrix, and bounds.'
        print *,'**********************************************'
        go to 100
      endif
c
c  Open the penalty vector file.
      call getvec
c  Set up the bvls matrices  aa  and  bb.
      call prep(nd,lvls,nd+3,lvls+2*nd+3,c(map),c(naa),c(nbb))
c  Future calls to bvls are warm starts.
      key=1
c
c------------Loop for minimization and maximization.----------------
c
      do 600 j=0,1
c  Loop for target values of objective functional.
      do 500 i=1,nv
c  Get the next penalty vector.
      read (id,*),(c(npv+ii),ii=0,lvls-1)
c  Load the (scaled) vector into  aa and the target value into  bb.
      call prepv(lvls,nd+3,lvls+2*nd+3,c(npv),c(nbl),c(nbu),c(naa),
     +  c(nbb),j)
c  Update the upper bounds on the slacks for the objective value.
      c(nbu+lvls+2*nd+1)=4.0*(pmax-pmin)
      c(nbu+lvls+2*nd+2)=4.0*(pmax-pmin)
c  Loop to achieve target value of chi---changes to  scale .
      do 300 k=1,20
        call scalobj(nd+3,lvls+2*nd+3,c(naa))
        call bvls(key,nd+3,lvls+2*nd+3,c(naa),c(nbb),c(nbl),c(nbu),
     +   c(nz),c(nw),c(nact),c(nzz),istate,its)
c  Find the l-1 misfit for this model.
        call check(c(nz),lvls,c(map),nd,viol)
c  If the misfit is less than  chi , need to test whether the bound
c  was really found.  The optimal value must be achieved on the 
c  boundary of the constraint set.  Thus if the misfit constraint is
c  not active, the a priori bounds must be, and the value of the
c  objective functional should be the a priori bound.  If not, adjust 
c  the value of  scale  and return to the bvls call.
c  If the misfit is greater than  chi , the value of  scale  is too
c  large  (since the problem is already known to be feasible).  Adjust
c  the value of   scale  and try again.
c
c  First calculate the achieved value of the objective functional.
        call dot(lvls,c(nz),c(npv),bnd)
        if ( (j .eq. 0 .and. bnd .eq. pmin)
     +    .or. (j .eq. 1 .and. bnd .eq. pmax) 
     +    .or. ( (viol .gt. chi*(1.0-tol)) .and. 
     +         (viol .lt. chi*(1.0+tol)) ) ) 
     +    then
            go to 400
          else if (viol .lt. chi*(1.0-tol) ) then
            scale=2.0*scale
          else 
            scale=scale/2.0
        endif
        if (mode .gt. 2) print *,' l-1 misfit ', viol,
     +    ' Changing  scale  to ', scale
 300  continue
      print *,'**********************************************'
      print *,' Target misfit not achieved. '
      print *,' Target: ', chi, ' Actual: ', viol
      print *,' Scale may need to be changed. '
      print *,'**********************************************'
      go to 100
 400  continue
c  Print as indicated and store the extremal bound.
      if (mode.gt.1) then
        print '(/)'
        print *,'vector # ',i,' bound ',bnd,
     +  ' one-norm misfit ', viol,
     +  ' a priori bounds on objective:', pmin, pmax,
     +   ' bvls iterations ', its
      endif
c  Save the bound.
      c(nbnd+j*nv+i-1)=bnd
c  Forward model, print model for appropriate mode.
      if (mode.ge.3) call forwrd(c(nz),lvls,c(map),nd)
      if (mode.ge.4) print '(/a/(i4,6x,g10.4))',
     +  '  i        m(i)', (l,c(nz+l-1),l=1,lvls)
c  Next target.
 500  continue
c  Start the maximization with the first penalty vector.
      rewind (unit=id)
 600  continue
c  Close the penalty vector file.
      close (unit=id)
c  Print the results.
      print '(/a/a/(i4,8x,g10.4))','Extremal Bounds',
     +  ' vector #     min ', (j,c(nbnd+j-1),j=1,nv)
      print '(/a/(i4,8x,g10.4))',' vector #     max ',
     +  (j,c(nbnd+nv+j-1),j=1,nv)
      go to 100
      end
c*********************************************************************
      subroutine check(z, m, a, nr, viol )
c*********************************************************************
c  Calculates one-norm violation of the data constraints.
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      dimension z(m), a(nr,m)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      viol=0.0
      do 10 i=1,nd
      tx=0.0
      do 5 j=1,lvls
      tx=tx+a(i,j)*z(j)
5     continue
      viol=viol+abs(tx-dbar(i))/dsv(i)
10    continue
      return
      end
c**********************************************************************
      subroutine dot(n,a,b,d)
c**********************************************************************
c  Computes double precision dot product of n-vectors a and b.
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      dimension a(n),b(n)
      d=0.0
      do 10 i=1,n
 10   d=d+a(i)*b(i)
      return
      end
c*********************************************************************
      subroutine forwrd(z,m,a,nr)
c*********************************************************************
c  Calculates forward predictions from Z using matrix A, and
c  displays the predictions, dbar and dsv.
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      dimension z(m), a(nr,m)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      print '(/)'
      print *,'  i    prediction         datum            std. dev.  '
      do 10 i=1,nd
      tx=0.0
      do 5 j=1,m
      tx=tx+a(i,j)*z(j)
5     continue
10    print '(i4,2x,3(g15.7,3x))',i,tx,dbar(i),dsv(i)
      return
      end
c**********************************************************************
      subroutine getbnd(m,bl,bu)
c***********************************************************************
c  Fetches upper and lower bounds.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      common /print/mode,id
      dimension bl(m),bu(m)
      character*64 name
      print *,'Enter file name for upper and lower bounds on the',
     + 'model elements. Lower bound for each element should',
     + ' precede upper bound:  bl(i), bu(i).'
      read '(a64)', name
      open (unit=id, file=name)
      read (id,*) (bl(j),bu(j),j=1,m)
      close (unit=id)
      return
      end
c**********************************************************************
      subroutine getmap (a,m,n)
c**********************************************************************
c  Fetches the data mapping matrix.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      common/print/mode,id
      dimension a(m,n)
      character*64 name
      print '(/a/a)','Enter file name for data mapping matrix.',
     + 'File should be in row-by-row format.'
      read '(a64)', name
      open (unit=id, file=name)
      read (id,*) ((a(i,j),j=1,n),i=1,m)
      close (unit=id)
      return
      end
c**********************************************************************
      subroutine getvec
c***********************************************************************
c  Opens the file containing the penalty vectors.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /print/mode,id
      character*64 name
      print *,''
      print *,'Enter file name for penalty vectors.'
      read '(a64)', name
      open (unit=id, file=name)
      return
      end
c*********************************************************************
      subroutine prep(m, n, mm, nn, a, aa, bb)
c*********************************************************************
c  This subroutine builds the greater part of the  bvls  arrays aa, bb.
c
C  The matrix  aa   and the ``data'' vector  bb  look like 
c
c  [ A  I  -I  0  0  0]                         [ dbar ]
c  [ 0  1   1  1  0  0]                         [ chi  ]
c  [pv  0   0  0  1 -1]                         [ p    ]
c  [ 0  0   0  0  1  1]                         [ 0    ]
c
c  Each of the first nd rows of A and dbar is multiplied by the 
c  smallest standard deviation and divided by its individual standard 
c  deviation.  The identity matrices are multiplied by the smallest 
c  standard deviation.
c  This routine builds the first nd+1 rows of  aa  and  bb , and the 
c  last 0 in  bb; the last two rows of  aa  and the next-to-last 
c  element of  bb  are built by prepv.
c  The last row of aa is modified by scalobj.
c
c* Calls no other routines.
c
      implicit double precision (a-h,o-z)
      dimension a(m,n), aa(mm,nn), bb(mm)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
c
c  Initialize aa .
      do 10 i=1, mm
      do 10 j=1,nn
 10   aa(i,j)=0.0
c  Find the scaling.
      rmin=dsv(1)
      do 20 i=2,nd
 20   rmin=min(rmin,dsv(i))
      do 40 i=1,nd
         bb(i)=dbar(i)*rmin/dsv(i)
         aa(i,i+lvls)=rmin
         aa(i,i+lvls+nd)=-rmin
         aa(nd+1,i+lvls)=1.0
         aa(nd+1,i+lvls+nd)=1.0
         do 30 j=1,lvls
 30         aa(i,j)=a(i,j)*rmin/dsv(i)
 40   continue
      aa(nd+1,lvls+2*nd+1)=1.0
      bb(nd+1)=chi
      bb(nd+3)=0.0
      return
      end
c*********************************************************************
      subroutine prepv(n, m, nn, pv, bl, bu, aa, bb, key)
c*********************************************************************
c
c  Constructs the bottom two rows of  aa  using  scale  and the penalty
c  vector  pv.
c  Loads the scaled target value into the bottom of bb.
c  Finds a ``reasonable'' target value by computing a priori bounds
c  on pv.z given the upper and lower bounds on z.
c
c* Calls no other routines.
      implicit double precision (a-h, o-z)
      common /s/ scale, pmin, pmax
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      dimension pv(n), bl(nn), bu(nn), aa(m,nn), bb(m)
c
c  Find a priori limits on pv.z given that bl<=z<=bu.
c
      pmax=0.0
      pmin=0.0
      pnorm=0.0
      do 10 i=1,lvls
      pnorm=pnorm+pv(i)**2
      if (pv(i) .ge. 0.0 ) then
        pmax=pmax+bu(i)*pv(i)
        pmin=pmin+bl(i)*pv(i)
      else
        pmax=pmax+bl(i)*pv(i)
        pmin=pmin+bu(i)*pv(i)
      endif
 10   continue
c  Scale and enter the elements of pv into aa and the target into bb.
      sc=1.0/sqrt( pnorm )
      if (key .eq. 0) then
        bb(nd+2)=(pmin-3.0*(pmax-pmin))*sc
      else
        bb(nd+2)=(pmax+3.0*(pmax-pmin))*sc
      endif
      do 20 i=1,lvls
 20      aa (nd+2, i)=pv(i)*sc
c  Set up slacks for objective value.
      aa (nd+2, lvls+2*nd+2)=1.d0
      aa (nd+2, lvls+2*nd+3)=-1.d0
      return
      end
c*********************************************************************
      subroutine scalobj(m,n,aa)
c*********************************************************************
c  Sets up the elements of  aa  to minimize the deviation from the
c  objective value.
c
c* Calls no other routines.
c
      implicit double precision (a-h, o-z)
      common /s/ scale, pmin, pmax
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      dimension aa(m,n)
      aa (nd+3, lvls+2*nd+2)=scale
      aa (nd+3, lvls+2*nd+3)=scale
      return
      end
c*********************************************************************
      subroutine setdat
c*********************************************************************
c  Reads data, upper and lower bounds on the model, 
c       and defines the variables in /datset/.
c
c  nd  is the number data.
c  lvls  is the number of elements in the model
c  xl() is an array of values at which data have been measured.
c  dbar(),dsv()  data values and their standard deviations.
c  nv is the number of penalty vectors to be bounded.
c  chi is the acceptable weighted one-norm misfit.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      common/print/mode,id
      common /s/ scale, pmin, pmax
      character*64 name
      print '(/''Input data file name.''/
     + ''Maximum of 50 data in the form:''/
     + ''   datum        standard deviation '')'
      read '(a64)', name
      print *,'Enter # data, # of model elements (<500), ',
     +  ' # of penalty vectors to bound,', 
     +  ' and allowable one-norm misfit.'
      read *,nd,lvls,nv,chi
      open(unit=id,file=name)
      do 100 i=1,nd
 100  read (id,*) dbar(i),dsv(i)
      close (unit=id)
c  Print the data.
      print '(/,''Data to be inverted:'',/,2x,''i'',5x,
     + ''datum '',4x,''standard dev.'')'
      do 200 i=1,nd
  200 print '(i4,7x,2g12.6)',i,dbar(i),dsv(i)
      print *,'Maximum misfit ', chi
      return
      end
c*********************************************************************
      subroutine slacks(m,bl,bu)
c*********************************************************************
c  Sets up the upper and lower bounds on all but the last 2 slack 
c  variables to impose the data constraints.
C  The last two slacks are set up in prepv.
c* Calls no other routines.
      implicit double precision (a-h, o-z)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      dimension bl(m), bu(m)
c  The lower bounds on all the slacks are zero, the upper bounds are
c  chi.
      do 10 i=1, nd
         bl(lvls+i)=0.0
         bl(lvls+nd+i)=0.0
         bu(lvls+i)=chi
 10      bu(lvls+nd+i)=chi
      bl(lvls+2*nd+1)=0.0
      bu(lvls+2*nd+1)=chi
      return
      end
c*********************************************************************
      subroutine space
c*********************************************************************
c  Allocates space for storage of data mapping matrix, bounds, 
c  solutions, and objective functional.
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ dbar(50),dsv(50),chi,nd,lvls,nv
      common /icons/ nz,nbnd,map,npv,naa,nbb,nbl,nbu,nw,
     +  nzz,nact,ntot
      nz=1
      nbnd=nz+lvls+2*nd+3 
      map=nbnd+2*nv                    
      npv=map+lvls*nd
      naa=npv+lvls                     
      nbb=naa+(lvls+2*nd+3)*(nd+3)           
      nbl=nbb+nd+3                  
      nbu=nbl+lvls+2*nd+3
      nw=nbu+lvls+2*nd+3
      nzz=nw+lvls+2*nd+3
      nact=nzz+nd+3                 
      ntot=nzz+(nd+3)*(nd+5)     
      return
      end