1      subroutine fptrpe(m,mm,idim,n,nr,sp,p,b,z,a,aa,q,right)
2c  subroutine fptrpe reduces the (m+n-7) x (n-7) cyclic bandmatrix a
3c  to upper triangular form and applies the same givens transformations
4c  to the (m) x (mm) x (idim) matrix z to obtain the (n-7) 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),aa(n,4),q((n-7)*mm*idim)
12     *,
13     * right(mm*idim)
14      integer nr(m)
15c  ..local scalars..
16      real*8 co,pinv,piv,si,one
17      integer i,irot,it,ii,i2,i3,j,jj,l,mid,nmd,m2,m3,
18     * nrold,n4,number,n1,n7,n11,m1
19c  ..local arrays..
20      real*8 h(5),h1(5),h2(4)
21c  ..subroutine references..
22c    fpgivs,fprota
23c  ..
24      one = 1
25      if(p.gt.0.) pinv = one/p
26      n4 = n-4
27      n7 = n-7
28      n11 = n-11
29      mid = mm*idim
30      m2 = m*mm
31      m3 = n7*mm
32      m1 = m-1
33c  we determine the matrix (a) and then we reduce her to
34c  upper triangular form (r) using givens rotations.
35c  we apply the same transformations to the rows of matrix
36c  z to obtain the (mm) x (n-7) matrix g.
37c  we store matrix (r) into a and aa, g into q.
38c  the n7 x n7 upper triangular matrix (r) has the form
39c             | a1 '     |
40c       (r) = |    ' a2  |
41c             |  0 '     |
42c  with (a2) a n7 x 4 matrix and (a1) a n11 x n11 upper
43c  triangular matrix of bandwidth 5.
44c  initialization.
45      nmd = n7*mid
46      do 50 i=1,nmd
47        q(i) = 0.
48  50  continue
49      do 100 i=1,n4
50        a(i,5) = 0.
51        do 100 j=1,4
52          a(i,j) = 0.
53          aa(i,j) = 0.
54 100  continue
55      jper = 0
56      nrold = 0
57      do 760 it=1,m1
58        number = nr(it)
59 120    if(nrold.eq.number) go to 180
60        if(p.le.0.) go to 740
61c  fetch a new row of matrix (b).
62        n1 = nrold+1
63        do 140 j=1,5
64          h(j) = b(n1,j)*pinv
65 140    continue
66c  find the appropriate row of q.
67        do 160 j=1,mid
68          right(j) = 0.
69 160    continue
70        go to 240
71c  fetch a new row of matrix (sp)
72 180    h(5) = 0.
73        do 200 j=1,4
74          h(j) = sp(it,j)
75 200    continue
76c  find the appropriate row of q.
77        j = 0
78        do 220 ii=1,idim
79          l = (ii-1)*m2+(it-1)*mm
80          do 220 jj=1,mm
81            j = j+1
82            l = l+1
83            right(j) = z(l)
84 220    continue
85c  test whether there are non-zero values in the new row of (a)
86c  corresponding to the b-splines n(j,*),j=n7+1,...,n4.
87 240     if(nrold.lt.n11) go to 640
88         if(jper.ne.0) go to 320
89c  initialize the matrix (aa).
90         jk = n11+1
91         do 300 i=1,4
92            ik = jk
93            do 260 j=1,5
94               if(ik.le.0) go to 280
95               aa(ik,i) = a(ik,j)
96               ik = ik-1
97 260        continue
98 280        jk = jk+1
99 300     continue
100         jper = 1
101c  if one of the non-zero elements of the new row corresponds to one of
102c  the b-splines n(j;*),j=n7+1,...,n4,we take account of the periodicity
103c  conditions for setting up this row of (a).
104 320     do 340 i=1,4
105            h1(i) = 0.
106            h2(i) = 0.
107 340     continue
108         h1(5) = 0.
109         j = nrold-n11
110         do 420 i=1,5
111            j = j+1
112            l0 = j
113 360        l1 = l0-4
114            if(l1.le.0) go to 400
115            if(l1.le.n11) go to 380
116            l0 = l1-n11
117            go to 360
118 380        h1(l1) = h(i)
119            go to 420
120 400        h2(l0) = h2(l0) + h(i)
121 420     continue
122c  rotate the new row of (a) into triangle.
123         if(n11.le.0) go to 560
124c  rotations with the rows 1,2,...,n11 of (a).
125         do 540 irot=1,n11
126            piv = h1(1)
127            i2 = min0(n11-irot,4)
128            if(piv.eq.0.) go to 500
129c  calculate the parameters of the givens transformation.
130            call fpgivs(piv,a(irot,1),co,si)
131c  apply that transformation to the columns of matrix q.
132            j = 0
133            do 440 ii=1,idim
134               l = (ii-1)*m3+irot
135               do 440 jj=1,mm
136                 j = j+1
137                 call fprota(co,si,right(j),q(l))
138                 l = l+n7
139 440        continue
140c  apply that transformation to the rows of (a) with respect to aa.
141            do 460 i=1,4
142               call fprota(co,si,h2(i),aa(irot,i))
143 460        continue
144c  apply that transformation to the rows of (a) with respect to a.
145            if(i2.eq.0) go to 560
146            do 480 i=1,i2
147               i1 = i+1
148               call fprota(co,si,h1(i1),a(irot,i1))
149 480        continue
150 500        do 520 i=1,i2
151               h1(i) = h1(i+1)
152 520        continue
153            h1(i2+1) = 0.
154 540     continue
155c  rotations with the rows n11+1,...,n7 of a.
156 560     do 620 irot=1,4
157            ij = n11+irot
158            if(ij.le.0) go to 620
159            piv = h2(irot)
160            if(piv.eq.0.) go to 620
161c  calculate the parameters of the givens transformation.
162            call fpgivs(piv,aa(ij,irot),co,si)
163c  apply that transformation to the columns of matrix q.
164            j = 0
165            do 580 ii=1,idim
166               l = (ii-1)*m3+ij
167               do 580 jj=1,mm
168                 j = j+1
169                 call fprota(co,si,right(j),q(l))
170                 l = l+n7
171 580        continue
172            if(irot.eq.4) go to 620
173c  apply that transformation to the rows of (a) with respect to aa.
174            j1 = irot+1
175            do 600 i=j1,4
176               call fprota(co,si,h2(i),aa(ij,i))
177 600        continue
178 620     continue
179         go to 720
180c  rotation into triangle of the new row of (a), in case the elements
181c  corresponding to the b-splines n(j;*),j=n7+1,...,n4 are all zero.
182 640     irot =nrold
183         do 700 i=1,5
184            irot = irot+1
185            piv = h(i)
186            if(piv.eq.0.) go to 700
187c  calculate the parameters of the givens transformation.
188            call fpgivs(piv,a(irot,1),co,si)
189c  apply that transformation to the columns of matrix g.
190            j = 0
191            do 660 ii=1,idim
192               l = (ii-1)*m3+irot
193               do 660 jj=1,mm
194                 j = j+1
195                 call fprota(co,si,right(j),q(l))
196                 l = l+n7
197 660        continue
198c  apply that transformation to the rows of (a).
199            if(i.eq.5) go to 700
200            i2 = 1
201            i3 = i+1
202            do 680 j=i3,5
203               i2 = i2+1
204               call fprota(co,si,h(j),a(irot,i2))
205 680        continue
206 700     continue
207 720     if(nrold.eq.number) go to 760
208 740     nrold = nrold+1
209         go to 120
210 760  continue
211      return
212      end
213