1      subroutine fptrnp(m,mm,idim,n,nr,sp,p,b,z,a,q,right)
2c  subroutine fptrnp reduces the (m+n-7) x (n-4) matrix a to upper
3c  triangular form and applies the same givens transformations to
4c  the (m) x (mm) x (idim) matrix z to obtain the (n-4) x (mm) x
5c  (idim) matrix q
6c  ..
7c  ..scalar arguments..
8      real*8 p
9      integer m,mm,idim,n
10c  ..array arguments..
11      real*8 sp(m,4),b(n,5),z(m*mm*idim),a(n,5),q((n-4)*mm*idim),
12     * right(mm*idim)
13      integer nr(m)
14c  ..local scalars..
15      real*8 cos,pinv,piv,sin,one
16      integer i,iband,irot,it,ii,i2,i3,j,jj,l,mid,nmd,m2,m3,
17     * nrold,n4,number,n1
18c  ..local arrays..
19      real*8 h(7)
20c  ..subroutine references..
21c    fpgivs,fprota
22c  ..
23      one = 1
24      if(p.gt.0.) pinv = one/p
25      n4 = n-4
26      mid = mm*idim
27      m2 = m*mm
28      m3 = n4*mm
29c  reduce the matrix (a) to upper triangular form (r) using givens
30c  rotations. apply the same transformations to the rows of matrix z
31c  to obtain the mm x (n-4) matrix g.
32c  store matrix (r) into (a) and g into q.
33c  initialization.
34      nmd = n4*mid
35      do 50 i=1,nmd
36        q(i) = 0.
37  50  continue
38      do 100 i=1,n4
39        do 100 j=1,5
40          a(i,j) = 0.
41 100  continue
42      nrold = 0
43c  iband denotes the bandwidth of the matrices (a) and (r).
44      iband = 4
45      do 750 it=1,m
46        number = nr(it)
47 150    if(nrold.eq.number) go to 300
48        if(p.le.0.) go to 700
49        iband = 5
50c  fetch a new row of matrix (b).
51        n1 = nrold+1
52        do 200 j=1,5
53          h(j) = b(n1,j)*pinv
54 200    continue
55c  find the appropriate column of q.
56        do 250 j=1,mid
57          right(j) = 0.
58 250    continue
59        irot = nrold
60        go to 450
61c  fetch a new row of matrix (sp).
62 300    h(iband) = 0.
63        do 350 j=1,4
64          h(j) = sp(it,j)
65 350    continue
66c  find the appropriate column of q.
67        j = 0
68        do 400 ii=1,idim
69          l = (ii-1)*m2+(it-1)*mm
70          do 400 jj=1,mm
71            j = j+1
72            l = l+1
73            right(j) = z(l)
74 400    continue
75        irot = number
76c  rotate the new row of matrix (a) into triangle.
77 450    do 600 i=1,iband
78          irot = irot+1
79          piv = h(i)
80          if(piv.eq.0.) go to 600
81c  calculate the parameters of the givens transformation.
82          call fpgivs(piv,a(irot,1),cos,sin)
83c  apply that transformation to the rows of matrix q.
84          j = 0
85          do 500 ii=1,idim
86            l = (ii-1)*m3+irot
87            do 500 jj=1,mm
88              j = j+1
89              call fprota(cos,sin,right(j),q(l))
90              l = l+n4
91 500      continue
92c  apply that transformation to the columns of (a).
93          if(i.eq.iband) go to 650
94          i2 = 1
95          i3 = i+1
96          do 550 j=i3,iband
97            i2 = i2+1
98            call fprota(cos,sin,h(j),a(irot,i2))
99 550      continue
100 600    continue
101 650    if(nrold.eq.number) go to 750
102 700    nrold = nrold+1
103        go to 150
104 750  continue
105      return
106      end
107