c
c
c***********************************************************************
c
      subroutine DPA ( n, kg, m, mmax, npp2, iv, t, y, w,
c
     *                 weight, alf, isel, iprint, a,
c
     *                 u, r, rnorm )
c
c------------------------------------------------------------------------
c
c     compute the norm of the residual (if isel = 1 or 2),
c     or the (m-n) x kg derivative of the modified
c     residual (m-n) vector q2*y (if isel = 1 or 3).
c     here q * phi = s, i.e., 
c
c    n     ( q1 ) (     .   .        )   ( s  . r1 .  f1  ) 
c          (----) ( phi . y . d(phi) ) = (--- . -- . ---- ) 
c    m-n   ( q2 ) (     .   .        )   ( 0  . r2 .  f2  ) 
c
c            m       n    1      p         n     1     p
c
c     where q is m x m orthogonal, and s is n x n upper
c     triangular.  the norm of the residual = norm(r2),
c     and the desired derivative according to ref.(5), is
c
c                                         -1
c              d(q2 * y) = -q2 * d(phi)* s  * q1* y.
c
c  ......................................................
c
      PARAMETER ( maxn=100, maxkg=50 )
c
      save
c
      common /blkinc/ inc
c
      real a(mmax,npp2),alf(kg),t(mmax,iv),w(m),weight(m),y(m),r(m),u(n)
c
      integer firstc,firstr,p,inc(maxkg,maxn+1)
c
      logical phinp1
c
c-------------------------------------------------------------------------------
c
      np1 = n + 1
c    ----------------------------------------------------------------
      call ADA (np1,n,kg,m,mmax,npp2,iv,a,t,alf,min0(isel,3),phinp1,
     2          nconp1)
c    ----------------------------------------------------------------
c
      go to (1100,1090,1080,1080), isel 
c                                      isel = 3 or 4
1080  firstc=np2
      lastc=npp1
      firstr=(4-isel)*n+1
      go to 1200
c                                      isel = 2
1090  firstc=nconp1 
      lastc=np1
      if (ncon.eq.0) go to 1150
      if (a(1,ncon).eq.save) go to 1150 
      isel=-7
c    --------------------------------
      call VARERR (iprint,isel,ncon)
c    --------------------------------
      go to 1320
c                             isel = 1 (first time only)
1100  np2=n+2
      nkg2=n+2+kg
      npp1=npp2-1
      firstc=1
      lastc=npp1
      firstr=np1
      do 1110 i=1,m 
         if (weight(i).ge.0.) go to 1110
c                                      error in weights
         isel=-6
         ii=i
c       ------------------------------
         call VARERR (iprint,isel,ii)
c       ------------------------------
         go to 1320 
1110     w(i)=sqrt(weight(i)) 
c
      ncon=n
      nconp1=np1
      if (kg.eq.0) go to 1160 
c
      if (n.ne.0) go to 1120
      phinp1=.true. 
      go to 1180
c                      determine number of constant fcns.
1120  p=0 
      do 1130 j=1,np1
         if (p.eq.0) nconp1=j 
         do 1130 k=1,kg
            inckj=inc(k,j)
            if (inckj.eq.1) p=p+1
1130  continue
c
      ncon=nconp1-1 
c                determine if phi(n+1) is in the model.
      phinp1=.false.
      do 1140 k=1,kg
         if (inc(k,np1).eq.1) phinp1=.true.
1140  continue
c                                     isel = 1 or 2
1150  if (phinp1) go to 1180
1160  do 1170 i=1,m 
1170     r(i)=y(i)
      go to 1200
1180  do 1190 i=1,m 
1190     r(i)=y(i)-r(i)
c                             weight appropriate columns
1200  do 1210 i=1,m 
         acum=w(i)
         do 1210 j=firstc,lastc
1210        a(i,j)=a(i,j)*acum
c
c        compute orthogonal factorizations by householder
c        reflections.  if isel = 1 or 2, reduce phi
c        (stored in the first n columns of the matrix a)
c        to upper triangular form, (q*phi = s), and
c        transform y (stored in column n+1), getting
c        q*y = r.  if isel = 1, also transform j = d phi
c        (stored in columns n+2 through n+p+1 of the
c        matrix a), getting q*j = f.  if isel = 3 or 4,
c        phi has already been reduced, transform only j.
c        s, r, and f overwrite phi, y, and j, respec-
c        tively and a factored form of q is saved in u
c        and the lower triangle of phi. 
c
      if (n.eq.0) go to 1260
      do 1250 k=1,n 
         kp1=k+1
         if (isel.ge.3.or.(isel.eq.2.and.k.lt.nconp1)) go to 1230
         alpha=0.
         do 1220 i=k,m
1220        alpha=alpha+a(i,k)**2
         alpha=sign(sqrt(alpha),a(k,k)) 
         u(k)=a(k,k)+alpha
         a(k,k)=-alpha
         firstc=kp1 
         if (alpha.ne.0.0) go to 1230
         isel=-8
         kk=k
c      --------------------------------
         call VARERR (iprint,isel,kk)
c      --------------------------------
         go to 1320 
c                              apply reflections to columns 
c                              firstc to lastc.
1230     beta=-a(k,k)*u(k)
         do 1250 j=firstc,lastc
            acum=u(k)*a(k,j)
            do 1240 i=kp1,m
1240           acum=acum+a(i,k)*a(i,j)
            acum=acum/beta
            a(k,j)=a(k,j)-u(k)*acum
            do 1250 i=kp1,m
1250           a(i,j)=a(i,j)-a(i,k)*acum
c
1260  if (isel.ge.3) go to 1280
      rnorm=0.
      do 1270 i=np1,m
1270     rnorm=rnorm+r(i)**2
      rnorm=sqrt(rnorm)
      if (isel.eq.2) go to 1320
      if (ncon.gt.0) save=a(1,ncon)
c
c        f2 is now contained in rows n+1 to m and columns
c        n+2 to n+p+1 of the matrix a.  now solve the
c        n x n upper triangular system s*c = r1 for the
c        linear parameters c.  c overwrites r1.
c
c                 --------------------------
1280  if (n.gt.0)  call BACSUB (mmax,n,a,r)
c                 --------------------------
c
c        major part of kaufman's simplification occurs
c        here.  compute the derivative of eta with
c        respect to the nonlinear parameters
c
c    t  d eta      t   n       d phi(j)   d phi(n+1)
c   q *-------- = q *(sum c(j) -------- + ----------) =  f2*c
c      d alf(k)        j=1      d alf(k)   d alf(k)
c
c        and store the result in columns n+2 to n+kg+1.
c        if isel not = 4, the first n rows are omitted.
c        this is -d(q2)*y.  if isel not = 4 the residual
c        r2 = q2*y (in col. n+1) is copied to column
c        n+kg+2.  otherwise all of column n+1 is copied.
c
      do 1310 i=firstr,m
         if (n.eq.ncon) go to 1310
         l=np1
         do 1300 k=1,kg
            acum=0. 
            do 1290 j=nconp1,n
               if (inc(k,j).eq.0) go to 1290
               l=l+1
               acum=acum+a(i,l)*r(j)
1290        continue
            ksub=np1+k
            if (inc(k,np1).eq.0) go to 1300
            l=l+1
            acum=acum+a(i,l)
1300        a(i,ksub)=acum
1310     a(i,nkg2)=r(i)
c
1320  return
c
      end 
