subroutine fptrnp(m,mm,idim,n,nr,sp,p,b,z,a,q,right)
c  subroutine fptrnp reduces the (m+n-7) x (n-4) matrix a to upper
c  triangular form and applies the same givens transformations to
c  the (m) x (mm) x (idim) matrix z to obtain the (n-4) x (mm) x
c  (idim) matrix q
c  ..
c  ..scalar arguments..
      real*8 p
      integer m,mm,idim,n
c  ..array arguments..
      real*8 sp(m,4),b(n,5),z(m*mm*idim),a(n,5),q((n-4)*mm*idim),
     * right(mm*idim)
      integer nr(m)
c  ..local scalars..
      real*8 cos,pinv,piv,sin,one
      integer i,iband,irot,it,ii,i2,i3,j,jj,l,mid,nmd,m2,m3,
     * nrold,n4,number,n1
c  ..local arrays..
      real*8 h(7)
c  ..subroutine references..
c    fpgivs,fprota
c  ..
      one = 1
      if(p.gt.0.) pinv = one/p
      n4 = n-4
      mid = mm*idim
      m2 = m*mm
      m3 = n4*mm
c  reduce the matrix (a) to upper triangular form (r) using givens
c  rotations. apply the same transformations to the rows of matrix z
c  to obtain the mm x (n-4) matrix g.
c  store matrix (r) into (a) and g into q.
c  initialization.
      nmd = n4*mid
      do 50 i=1,nmd
        q(i) = 0.
  50  continue
      do 100 i=1,n4
        do 100 j=1,5
          a(i,j) = 0.
 100  continue
      nrold = 0
c  iband denotes the bandwidth of the matrices (a) and (r).
      iband = 4
      do 750 it=1,m
        number = nr(it)
 150    if(nrold.eq.number) go to 300
        if(p.le.0.) go to 700
        iband = 5
c  fetch a new row of matrix (b).
        n1 = nrold+1
        do 200 j=1,5
          h(j) = b(n1,j)*pinv
 200    continue
c  find the appropriate column of q.
        do 250 j=1,mid
          right(j) = 0.
 250    continue
        irot = nrold
        go to 450
c  fetch a new row of matrix (sp).
 300    h(iband) = 0.
        do 350 j=1,4
          h(j) = sp(it,j)
 350    continue
c  find the appropriate column of q.
        j = 0
        do 400 ii=1,idim
          l = (ii-1)*m2+(it-1)*mm
          do 400 jj=1,mm
            j = j+1
            l = l+1
            right(j) = z(l)
 400    continue
        irot = number
c  rotate the new row of matrix (a) into triangle.
 450    do 600 i=1,iband
          irot = irot+1
          piv = h(i)
          if(piv.eq.0.) go to 600
c  calculate the parameters of the givens transformation.
          call fpgivs(piv,a(irot,1),cos,sin)
c  apply that transformation to the rows of matrix q.
          j = 0
          do 500 ii=1,idim
            l = (ii-1)*m3+irot
            do 500 jj=1,mm
              j = j+1
              call fprota(cos,sin,right(j),q(l))
              l = l+n4
 500      continue
c  apply that transformation to the columns of (a).
          if(i.eq.iband) go to 650
          i2 = 1
          i3 = i+1
          do 550 j=i3,iband
            i2 = i2+1
            call fprota(cos,sin,h(j),a(irot,i2))
 550      continue
 600    continue
 650    if(nrold.eq.number) go to 750
 700    nrold = nrold+1
        go to 150
 750  continue
      return
      end