program sbvq
c
c
c  Strict bounds on Linear Functionals of Bounded Variables Subject
c     to a Quadratic Misfit Criterion on Linear Data Constraints
c
c  P.B. Stark      Version 1/88
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  a specified weighted sum-of-squares tolerance.
c
c  i.e.      min and max  p .x     j=1,...,nv
c                          j
c                               2       2
c      such that   || A.x - d ||  <= chi     (+- small tolerance chitol)
c                               2
c
c         and      bl(i) <= x(i) <= bu(i) , i=1,...,lvls
c
c  for a user provided list of vectors p .
c                                       j
c
c  x, bl, bu, and p  are  lvls-vectors,  A  is an  nd  by  lvls  matrix,
c                  j
c  d  is an  nd-vector.
c
c
c  Structured so that the user may easily replace many IO routines with
c  subroutines to generate the mapping matrix, a priori bounds, and
c  objective vectors.
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, DOT, FITMIS, FORWRD, GETBND,
c* GETMAP, GETVEC, PREP1, PREP2, SETDAT, and SPACE.
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                                  mean datum    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     relative weighting of      1.d-4 seems to work consistently.
c       penalty functional       1.e-3 works in CRAY single precision.
c
c     expected chi-squared       set for desired confidence level for 
c                                  number of data used
c
c     tolerance for fitting      set to about .001 chi-squared
c       chi
c
c     initial step length for    make this big enough to be sure chi**2
c       sweeping objective         is violated in 2-3 steps
c       values
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=50000, nimax=500, maxit=500)
c  nmax: maximum allocated array storage.
c  maxit: maximum number of iterations when hunting for bounds.
      common /c/ c(nmax)
      common /srch/ scale,chi20,chitol,del0,scz
c  scale: relative weight of constraints versus penalty functional.
c  chi20: expected chi**2 to achieve.
c  chitol: tolerance for fitting chi20.
c  del0: initial step size for the objective constraint.
c
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
c
c  nd: number of data.
c  lvls: number of elements in the model.
c  dbar: data means.
c  dsv: data standard deviations.
c  nv: number penalty vectors for which to find extremal bounds.
c  zne: objective value the model is constrained to achieve.
c
      common /icons/ nz0,nz,nbnd,map,npv,naa,nbb,nbl,nbu,
     +  nw,nzz,nact,ntot
c
c  nz0: starting location of best fitting model in /c/.
c  nz:  starting location of the current solution in /c/.
c  nd:  starting location of normalized data and target depth.
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/
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
      print '(/)'
      print *,'Real array space required ',ntot,' reserved ',nmax
      print *,'Integer space required ',lvls+1, ' reserved ',nimax
      if(ntot.gt.nmax.or.lvls+1.gt.nimax) go to 100
      print '(/)'
c  Find best-fitting model for reference.
c
c  Get a priori upper/lower bounds on model, and data mapping matrix.
      call getbnd(lvls,c(nbl),c(nbu))
      call getmap(c(map),nd,lvls)
c  Open penalty vector file.
      call getvec
c  Form all scaled vectors and matrices for bvls.
      call prep1(nd,lvls,c(map),c(naa),c(nbb))
c  Solve.
      key=0
      call bvls(key,nd,lvls,c(naa),c(nbb),c(nbl),c(nbu),c(nz0),
     +  c(nw),c(nact),c(nzz),istate,its)
      chi2=fitmis(nd,lvls,c(map),c(nz0))
      print *,' minimum misfit ', chi2, ' bvls misfit ', c(nw)
      if (chi2.gt.chi20+chitol) then
         print *, ' *** infeasible problem *** '
         go to 100
      endif
C  Forward model, print best-fitting model if requested.
      if (mode.ge.3) call forwrd (c(nz0), lvls, c(map), nd)
      if (mode.eq.4) print '(/a/(i4,6x,g10.4))',
     +  '  i        m(i)', (l,c(nz0+l-1),l=1,lvls)
      key=1
c  Loop for minimization and maximization.
      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  Find the objective value achieved by the best-fitting model.
      call dot(lvls,c(npv),c(nz0),bnd0)
      print '(/)'
      print *,'vector # ',i,' best-fitting objective ',bnd0
c  Find max and min possible penalty values for a model satisfying
c  the a priori bounds.
      call dot(lvls,c(npv),c(nbl),bndlo)
      call dot(lvls,c(npv),c(nbu),bndhi)
c  Set new bvls parameters, etc.
      call prep2(nd,lvls,c(map),c(npv),c(naa),c(nbb))
c  Initialize the search direction.
      del=del0*(-1.0)**(j+1)
c  Set the target objective value; scale and load into bb.
      zne=bnd0+del
      c(nbb)=zne*scz
C  Outer loop with constant step size.
      do 300 it1=1,maxit
      call bvls(key,nd+1,lvls,c(naa),c(nbb),c(nbl),c(nbu),c(nz),
     + c(nw),c(nact),c(nzz),istate,its)
c  Find the misfit.
      chi2=fitmis(nd,lvls,c(map),c(nz))
c  Compute the extremal bound.
      call dot(lvls,c(nz),c(npv),bnd)
      if (mode.gt.1)
     +  print *,'target value ',zne,'vector # ',i,' achieved ',bnd,
     +   ' chi**2 ',chi2,' bvls iterations ', its, 
     +   ' bvls misfit ', c(nw)
      if ( abs(chi2-chi20) .le. chitol .or.
     +  (chi2 .lt. chi20 .and. (zne .le. bndlo .or.
     +     zne .ge. bndhi ))) then
c  Within the tolerance--save and print.
        go to 400
c  If chi2 < chi20 take another step.
      else if (chi2.lt.chi20) then
        zne=zne+del
        c(nbb)=zne*scz
c  If chi2 > chi20, [ zne, zne-del ] has a root: bisect.
      else 
        del=del/2.0
        zne=zne-del
        c(nbb)=zne*scz
c  Inner loop for bisection.
        do 200 it=1,maxit
          call bvls(key,nd+1,lvls,c(naa),c(nbb),c(nbl),c(nbu),
     +      c(nz),c(nw),c(nact),c(nzz),istate,its)
          chi2=fitmis(nd,lvls,c(map),c(nz))
c  Compute and print the extremal bound.
          call dot(lvls,c(nz),c(npv),bnd)
          if (mode.gt.1)
     +      print *,'target ',zne,'vector # ',i,' achieved ',bnd,
     +       ' chi**2 ',chi2, ' bvls iterations ', its, 
     +       ' bvls misfit ', c(nw)
          del=del/2.0
          if ( abs(chi2-chi20) .le. chitol  .or.
     +  (chi2 .lt. chi20 .and. (zne .le. bndlo .or.
     +     zne .gt. bndhi ))) then
c  Within the tolerance--save and print.
            go to 400
          else if (chi2.gt.chi20) then
c  [zne,zne-2*del] contains a root.
            zne=zne-del
            c(nbb)=zne*scz
          else
c  [zne,zne+2*del] contains a root.
            zne=zne+del
            c(nbb)=zne*scz
          end if
c  Iterate again.
 200    continue
        print '(/a)','Maximum number of steps exceeded.'
        go to 450
      end if
 300  continue
      print '(/a)','Maximum number of steps exceeded.'
      go to 450
c  Save and print the answer for this target.
 400  print '(/a/)','Bound found:'
      print *,'target ',zne,'vector # ',i,' achieved ',bnd,
     + ' chi**2 ',chi2, ' bvls iterations ', its, 
     + ' bvls misfit ', c(nw)
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.
 450  continue
 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)
      stop
      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**********************************************************************
      double precision function fitmis(m,n,a,x)
cs    real function fitmis(m,n,a,x)
c**********************************************************************
c                       2
c  Computes ||a.x-dbar||     .
c                       dsv
c
      implicit double precision (a-h, o-z)
      dimension a(m,n), x(n)
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
      fitmis=0.0
      do 20 i=1,nd
      p=0.0
      do 10 j=1,lvls
10    p=p+a(i,j)*x(j)
20    fitmis=fitmis+((p-dbar(i))/dsv(i))**2
      return
      end
c*********************************************************************
      subroutine forwrd(zp,m,a,nr)
c*********************************************************************
c  Calculates forward predictions from ZP using matrix A, and
c  displays the predictions between TXLO and TXHI at each p in PL.
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      dimension zp(m), a(nr,m)
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
      print '(/)'
      print *,'  i         sigma         prediction       mean  '
      do 10 i=1,nd
      tx=0.0
      do 5 j=1,m
      tx=tx+a(i,j)*zp(j)
5     continue
10    print '(i4,2x,3(g15.7,3x))',i,dsv(i),tx,dbar(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/ nd,lvls,dbar(100),dsv(100),nv,zne
      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/ nd,lvls,dbar(100),dsv(100),nv,zne
      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 *,'Enter file name for penalty vectors.'
      read '(a64)', name
      open (unit=id, file=name)
      return
      end
c*********************************************************************
      subroutine prep1( m, n, a, aa, bb )
c*********************************************************************
c  Makes normalized matrices and vectors for bvls for
c  best-fitting model problem.
C  Normalizes the data and rows of the mapping by their relative errors.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
      common /srch/ scale,chi20,chitol,del0,scz
      dimension a(m,n),aa(m,n),bb(m)
c
c  Find the normalization for the data constraints.
      rmin=dsv(1)
      do 10 i=2,m
 10   rmin=min(rmin,dsv(i))
c  Enter the mapping a, properly scaled, into rows 1 to m of aa.
      do 50 i=1, m
c  Build bb from the scaled data.
      bb(i)=rmin*dbar(i)/dsv(i)
      do 50 j=1, n
50    aa(i,  j)= a(i,j)*rmin/dsv(i)
      return
      end
c*********************************************************************
      subroutine prep2(m, n, a, pv, aa, bb )
c*********************************************************************
c  Makes normalized matrix aa and vector bb for bvls from data mapping 
c  matrix a, penalty vector pv, data vectors dbar and dsv, and the
c  scale factor scale.
c  bb(1) is left undefined; scz is returned (in common) to allow 
c  bb(1) to be constructed later in the program.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
      common /srch/ scale,chi20,chitol,del0,scz
      dimension a(m,n),pv(n),aa(m+1,n),bb(m+1)
c
c  Find the normalization for the data constraints.
      rmin=dsv(1)
      do 10 i=2,m
 10   rmin=min(rmin,dsv(i))
c  Multiply by the smallest error so the rows do not get too big.
      sc= scale*rmin
c  Find the norm of the penalty vector.
      scz=0.0
      do 20 i=1,n
 20   scz=scz+pv(i)**2
      scz=1.e0/sqrt(scz)
c  Enter the scaled objective functional in row 1 of aa.
      do 40 i=1,n
40    aa(1,i)=scz*pv(i)
c
c  Rows 2 to m+1 of aa get a, scaled.
      do 50 i=1, m
c  Elements 2 to m+1 of bb are the scaled data; element 1 will get the
c    scaled target depth later.
      bb(i+1)=sc*dbar(i)/dsv(i)
      do 50 j=1, n
50    aa(i+1,  j)= a(i,j)*sc/dsv(i)
      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()  mean, standard deviation of measured data values.
c  nv is the number of penalty vectors to be bounded.
c
c* Calls no other routines.
      implicit double precision (a-h,o-z)
      common /datset/ nd,lvls,dbar(100),dsv(100),nv,zne
      common/print/mode,id
      common /srch/ scale,chi20,chitol,del0,scz
      character*64 name
      print '(/''Input data file name.''/
     + ''Maximum of 100 data in the form:''/
     + ''   mean    standard deviation '')'
      read '(a64)', name
      print *,'Enter # data, # of model elements (<500), ',
     +  '# of penalty vectors to bound, penalty weight, ',
     +  'chi**2, chi tolerance, and step length'
      read *,nd,lvls,nv,scale,chi20,chitol,del0
      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,
     + ''mean'',8x,''sigma'')'
      do 200 i=1,nd
  200 print '(i4,2x,2g12.6)',i,dbar(i),dsv(i)
      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/ nd,lvls,dbar(100),dsv(100),nv,zne
      common /icons/ nz0,nz,nbnd,map,npv,naa,nbb,nbl,nbu,nw,
     +  nzz,nact,ntot
      nz0=1
      nz=nz0+lvls                      
      nbnd=nz+lvls                     
      map=nbnd+2*nv                    
      npv=map+lvls*nd
      naa=npv+lvls                     
      nbb=naa+lvls*(nd+1)           
      nbl=nbb+nd+1                  
      nbu=nbl+lvls
      nw=nbu+lvls
      nzz=nw+lvls                      
      nact=nzz+nd+1                 
      ntot=nzz+(nd+1)*(nd+3)     
      return
      end