c
c
c***********************************************************************
c
      subroutine VARPRO ( n, kg, m, mmax, npp2, iv, t, y, w,
c
     *                    weight, a, iprint, alf, c,
c
     *                    ierr, itmax, eps1 )
c
c***********************************************************************
c
c  the following comments mostly refer to
c  the batch version of VARPRO.
c
c  given a set of m observations, consisting of
c  values y(1), y(2),...,y(m) of a dependent variable
c  y, where y(i) corresponds to the iv independent
c  variable(s) t(i,1), t(i,2),..., t(i,iv), VARPRO
c  attempts to compute a weighted least squares fit to
c  a function eta (the 'model') which is a linear 
c  combination
c                 n 
c    eta(c,alf;t)  =  sum c * phi (alf;t) + phi  (alf;t)
c                j=1  j     j             n+1
c
c  of nonlinear functions phi(j) (e.g. a sum of
c  exponentials and/or gaussians),  that is, determine
c  the linear parameters c(j) and the vector of
c  nonlinear parameters alf by minimizing
c
c                2     m                           2
c  norm(residual)  =  sum weight * (y - eta(c,alf;t )) .
c                     i=1       i    i             i
c
c  the (n+1)-st term is optional, and is used when it
c  is desired to fix one or more of the c's (rather
c  than let them be determined).  VARPRO requires 
c  first derivatives of the phi's.
c
c                  notes: 
c
c  a)  the above problem is also referred to as
c  'multiple nonlinear regression'.  for use in
c  statistical estimation, VARPRO returns the
c  residuals, the covariance matrix of the linear and
c  nonlinear parameters, and the estimated variance
c  of the observations.  other options are available
c  in the interactive version.
c
c  b)  an eta of the above form is called 'separable'.
c  the case of a nonseparable eta can be handled by
c  setting n = 0 and using phi(n+1).
c
c  c)  VARPRO may also be used to solve linear least
c  squares problems (in that case no iterations are
c  performed).  set kg = 0.
c
c  d)  the main advantage of VARPRO over other least
c  squares programs is that no initial guesses are
c  needed for the linear parameters.  not only does
c  this make it easier to use, but it often leads 
c  to faster convergence.
c
c
c   description of parameters 
c
c  n       number of linear parameters c (must be 
c          .ge. 0). 
c  kg      number of nonlinear parameters alf (must
c          be .ge. 0).
c  m       number of observations, m must be greater
c          than kg + n (i.e., the number of observa-
c          tions must exceed the number of parameters).
c  iv      number of independent variables t.
c  t       real m by iv array of independent variables.
c          t(i,j) contains the value of the i-th
c          observation of the j-th independent
c          variable.
c  y       m-vector of observations, one for each 
c          row of t.
c  w       m-vector of sqrt of weights. 
c  weight  m-vector of nonnegative weights.  should
c          be set to 1's if weights are not desired.
c  inc     kg x (n+1) integer incidence matrix.
c          inc(k,j) = 1 if nonlinear parameter alf(k)
c          appears in the j-th function phi(j),
c          otherwise, inc(k,j) = 0.  if phi(n+1) is
c          not included in the model, the (n+1)-st
c          column of inc should be set to zero.  if
c          phi(n+1) is included, the appropriate
c          elements of the (n+1)-st column should be
c          set to 1's.  inc is not needed when n = 0 or
c          kg = 0.  caution:  the declared row dimension
c          of inc (in ADA) must currently be set to 12.
c          see 'restrictions' below.
c  mmax    the declared row dimension of the matrices
c          a and t.  it must be at least max(m,2*kg+3).
c  npp2    n+p+2, where p is the number of ones in
c          matrix inc..  the declared column
c          dimension of a must be at least npp2.
c          (if n = 0, set npp2 = kg+2.  if kg = 0,
c          set npp2 = n+2.)
c  a       real matrix of size max(m,2*kg+3) by n+p+2.
c          on input it contains the phi(j)'s and
c          their derivatives (see below).  on output
c          the first n+kg rows and columns of a will
c          contain an approximation to the (weighted)
c          covariance matrix at the solution (the 
c          first n rows correspond to the linear
c          parameters, the last kg to the nonlinear
c          ones), column n+kg+1 will contain the
c          weighted residuals (y - eta), a(1,n+kg+2)
c          will contain the (euclidean) norm of the
c          weighted residual, and a(2,n+kg+2) will
c          contain an estimate of the (weighted)
c          variance of the observations,
c          norm(residual)**2/(m-n-kg).
c  iprint  input integer controlling printed output.
c          if iprint is positive, the nonlinear
c          parameters, the weighted norm of the residual, and
c          the marquardt parameter will be output 
c          every iprint-th iteration (and initially,
c          and at the final iteration).  the linear
c          parameters will be printed at the final
c          iteration.  any error messages will also
c          be printed.  (iprint = 1 is recommended
c          at first.)  if iprint = 0,  only the final
c          quantities will be printed, as well as any
c          error messages.  if iprint = -1, no
c          printing will be done.  the user is then
c          responsible for checking the parameter 
c          ierr for errors.
c  alf     kg-vector of estimates of nonlinear
c          parameters (input).  on output it will 
c          contain optimal values of the nonlinear
c          parameters.
c  c       n-vector of linear parameters (output only).
c  ierr    integer error flag (output): 
c          .gt. 0 - successful convergence. ierr is
c          the number of iterations taken.
c          -1  terminated for too many iterations.
c          -2  terminated for ill-conditioning
c              (marquardt parameter too large.) also
c              see ierr = -8 below.
c          -4  input error in parameter n, m, kg, 
c              npp2 or mmax.
c          -5  inc matrix improperly specified, or
c              p disagrees with npp2.
c          -6  weights not nonnegatve.
c          -7  'constant' column was computed more
c              than once.
c          -8  catastrophic failure - a column of the
c              a matrix has become zero.  see
c              'convergence failures' below.
c
c          (if ierr .le. -4, the linear parameters,
c          covariance matrix, etc. are not returned.)
c
c   subroutines required
c
c  seven subroutines, DPA, ORFAC1, ORFAC2, BACSUB,
c  POSTPR, COV, VARERR and ADA are provided.  in addition,
c  the user must provide subroutines CALC, DERIV and WT,
c  which, given alf, will evaluate the functions phi(j) and 
c  their partial derivatives d phi(j)/d alf(k), at the sample
c  points t(i).
c
c  restrictions
c
c        the subroutines DPA and ADA contain the locally
c     dimensioned matrix inc, whose dimensions are
c     currently set for maxima of n+1 = 15, kg = 12.
c     they must be changed for larger problems.
c     data placed in array a is overwritten ('destroyed').
c     data placed in arrays t, y and inc is left intact.
c     the program runs in watfiv, except when n = 0
c     or kg = 0.
c
c        it is assumed that the matrix phi(j, alf; t(i))
c     has full column rank.  this means that the first n
c     columns of the matrix a must be linearly
c     independent.
c
c  convergence failures
c
c        if convergence failures occur, first check
c     for incorrect replies to which terms of the 
c     model depend on which nonlinear parameter.
c     if these are correct,
c     try several starting guesses for alf.  if ADA
c     is coded correctly, and if error returns ierr = -2
c     or -8 persistently occur, this is a sign of ill-
c     conditioning, which may be caused by several things.
c     one is poor scaling of the parameters; another is
c     an unfortunate initial guess for the parameters,
c     still another is a poor choice of the model.
c
c  algorithm
c
c        the residual r is modified to incorporate, for
c     any fixed alf, the optimal linear parameters for
c     that alf.  it is then possible to minimize only on
c     the nonlinear parameters.  after the optimal values
c     of the nonlinear parameters have been determined,
c     the linear parameters can be recovered by linear
c     least squares techniques (see ref. 1).
c
c        the minimization is by a modification of 
c     osborne's (ref. 3) modification of the levenberg-
c     marquardt algorithm.  instead of solving the normal
c     equations with matrix
c
c           t      2
c         (j j + nu  * d),    where  j = d(eta)/d(alf),
c
c     stable orthogonal (householder) reflections are
c     used on a modification of the matrix
c
c                           (   j  )
c                           (------) ,
c                           ( nu*d )
c
c     where d is a diagonal matrix consisting of the
c     lengths of the columns of j.  this marquardt
c     stabilization allows the routine to recover from
c     some rank deficiencies in the jacobian.  osborne's
c     empirical strategy for choosing the marquardt
c     parameter has proven reasonably successful in
c     practice.  (gauss-newton with step control can be
c     obtained by making the change indicated before the
c     instruction labeled 1).  a description can be found
c     in ref. (3), and a flow chart in (2), p. 22.
c
c     for reference, see
c
c     1.  gene h. golub and v. pereyra, 'the differen-
c         tiation of pseudo-inverses and nonlinear least
c         squares problems whose variables separate,'
c         siam j. numer. anal. 10, 413-432 (1973).
c     2.  ------, same title, stanford c.s. report 72-261,
c         feb. 1972.
c     3.  osborne, michael r., 'some aspects of non-
c         linear least squares calculations,' in lootsma,
c         ed., 'numerical methods for non-linear
c         optimization,' academic press, london, 1972.
c     4.  krogh, fred, 'efficient implementation of a
c         variable projection algorithm for nonlinear
c         least squares problems,' comm. acm 17, pp. 167-
c         169 (march, 1974).
c     5.  kaufman, linda, 'a variable projection method
c         for solving separable nonlinear least squares
c         problems', b.i.t. 15,49-57 (1975).
c
c               john bolstad
c               computer science dept., serra house
c               stanford university
c               september, 1976
c
c  ......................................................
c
      real a(mmax,npp2),c(n),alf(kg),t(mmax,iv),weight(m),w(m),y(m),nu
      integer b1,output
      logical skip,ip
      common /blokip/ ip
      data output /6/
c
c  the following three parameters are used in
c  the convergence test:  eps1 is an absolute
c  tolerance for the norm of the projection of
c  the residual onto the range of the jacobian
c  of the variable projection functional.  eps2
c  is a relative and absolute tolerance for the
c  change in the norm of the residual.  itmax is
c  the maximum number of function and derivative
c  evaluations allowed.  caution:  eps1 and
c  eps2 must not be set smaller than 10 times
c  the unit round-off of the machine.
c
      eps2=eps1
      ierr=1
      iter=0
      np1=n+1
      b1=n+2
      nkg2=n+kg+2
      kgp1=kg+1
      kg23=2*kg+3
      skip=.false.
      modit=iprint
      if (iprint.le.0) modit=itmax+2
      nu=0.
c
c  if gauss-newton is desired, remove the
c  next statement.
c
      nu=1.
c
c  check for valid input
c
      if (n.ge.0.and.kg.ge.0.and.n+kg.lt.m.and.nkg2.le.npp2.and.
     2    kg23.le.mmax.and.m.le.mmax.and.iv.gt.0.and..not.(kg.eq.0.and.
     3    n.eq.0)) goto 710
      ierr=-4
c    -----------------------------
      call VARERR (iprint,ierr,1)
c    -----------------------------
      go to 880
c
c           begin outer iteration loop to update alf.
c           calculate the norm of the residual and the
c           derivative of the modified residual the
c           first time, but only the derivative in
c           subsequent iterations.
c
c    -----------------------------------------------------------------
710   call DPA (n,kg,m,mmax,npp2,iv,t,y,w,weight,alf,ierr,iprint,a,c, 
     2          a(1,np1),rnorm)
c   ------------------------------------------------------------------
      gnstep=1.0
      iterin=0
      if (iter.gt.0) go to 720
      if (kg.eq.0) go to 860
      if (ierr.ne.1) go to 880
c
      if (iprint.le.0) go to 720
      write (output,900) iter
      if (ip) write (10,900) iter
      write (output,940) (a(k,b1),k=1,kg)
      if (ip) write (10,940) (a(k,b1),k=1,kg)
      write (output,920) iterin,rnorm
      if (ip) write (10,920) iterin,rnorm
      write (output,890) nu
      if (ip) write (10,890) nu
c
c       begin two-stage orthogonal factorization
c    --------------------------------------------------------
720   call ORFAC1 (kgp1,mmax,m,n,iprint,a(1,b1),prjres,ierr)
c    --------------------------------------------------------
c
      if (ierr.lt.0) go to 880
      ierr=2
      if (nu.eq.0.) go to 740 
c
c       begin inner iteration loop for generating 
c       new alf and testing it for acceptance.
c
c    ------------------------------------
730   call ORFAC2 (kgp1,mmax,nu,a(1,b1))
c    ------------------------------------
c
c       solve a kg x kg upper triangular system
c        for delta-alf.  the transformed residual 
c       (in col. nkg2 of a) is overwritten by
c       the result delta-alf. 
c
c    -----------------------------------------
740   call BACSUB (mmax,kg,a(1,b1),a(1,nkg2))
c    -----------------------------------------
c
      do 750 k=1,kg 
750      a(k,b1)=alf(k)+a(k,nkg2)
c
c       new alf(k) = alf(k) + delta-alf(k)
c
c            step to the new point new alf, and compute
c            the new norm of the residual.  new alf
c            is stored in column b1 of a.
c
c    -------------------------------------------------------------------
760   call DPA (n,kg,m,mmax,npp2,iv,t,y,w,weight,a(1,b1),ierr,iprint,a,
     2          c,a(1,np1),rnew)
c   --------------------------------------------------------------------
      if (ierr.ne.2) go to 880
      iter=iter+1
      iterin=iterin+1
      skip=mod(iter,modit).ne.0
      if (skip) go to 770
      write (output,900) iter 
      if (ip) write (10,900) iter
      write (output,940) (a(k,b1),k=1,kg)
      if (ip) write (10,940) (a(k,b1),k=1,kg)
      write (output,920) iterin,rnew
      if (ip) write (10,920) iterin,rnew
c
770   if (iter.lt.itmax) go to 780
      ierr=-1
c    -----------------------------
      call VARERR (iprint,ierr,1)
c    -----------------------------
      go to 870
780   if (rnew.lt.rnorm) go to 830
c
c             retract the step just taken
c
      if (nu.ne.0) go to 800
c             gauss-newton option only
      gnstep=0.5*gnstep
      if (gnstep.lt.eps1) go to 870
      do 790 k=1,kg 
790      a(k,b1)=alf(k)+gnstep*a(k,nkg2)
      go to 760
c        enlarge the marquardt parameter
800   nu=1.5*nu
      if (.not.skip) write (output,910) nu
      if ((.not.skip).and.ip) write (10,910) nu
      if (nu.le.100) go to 810
      ierr=-2
c    -----------------------------
      call VARERR (iprint,ierr,1)
c    -----------------------------
      go to 870
c        retrieve upper triangular form 
c        and residual of first stage.
810   do 820 k=1,kg 
         ksub=np1+k 
         do 820 j=k,kgp1
            jsub=np1+j
            isub=kgp1+j
820         a(k,jsub)=a(isub,ksub)
      go to 730
c
c       end of inner iteration loop
c
c     accept the step just taken
c
830   delr=rnorm-rnew
      rnorm=rnew
      alfnrm=0.
      dalfnm=0.
      do 840 k=1,kg 
         alf(k)=a(k,b1)
         alfnrm=alfnrm+alf(k)**2
840      dalfnm=dalfnm+a(k,nkg2)**2
      alfnrm=sqrt(alfnrm)
      dalfnm=gnstep*sqrt(dalfnm)
      acum=dalfnm/alfnrm
c
c        if iterin is greater than 1, a step was
c        retracted during this outer iteration.
c
      if (iterin.eq.1) nu=0.5*nu
      if (skip) go to 850
      write (output,890) nu
      if (ip) write (10,890) nu
      write (output,930) acum 
      if (ip) write (10,930) acum
850   ierr=3
      if (prjres.gt.eps1.or.delr.gt.(rnorm+1.d0)*eps2) go to 710
c        end of outer iteration loop
c
c        calculate final quantities -- linear parameters,
c        residuals, covariance matrix, etc.
c
860   ierr=iter
c                 ---------------------------------------------------
870   if (kg.gt.0) call DPA (n,kg,m,mmax,npp2,iv,t,y,w,weight,alf,4,
     2                       iprint,a,c,a(1,np1),rnorm)
c                 ---------------------------------------------------
c
c    -----------------------------------------------------------
      call POSTPR (n,kg,m,nkg2,rnorm,iprint,alf,w,a,a(1,np1),c)
c    -----------------------------------------------------------
c
880   return
c
890   format (1h ,29x,'nu =',e16.7)
900   format (/1x/3x,'iteration',i4,4x,'nonlinear ',10hparameters)
910   format (1h ,4x,'step retracted, nu =',1pe15.7)
920   format (' '/i5,29h  weighted norm of residual =,1pe16.7)
930   format (' ',4x,'norm(delta-alf) / norm(alf) =',1pe12.3)
940   format (' '/' ', 1p4e18.7 )
c
      end 
