subroutine fptrpe(m,mm,idim,n,nr,sp,p,b,z,a,aa,q,right) c subroutine fptrpe reduces the (m+n-7) x (n-7) cyclic bandmatrix a c to upper triangular form and applies the same givens transformations c to the (m) x (mm) x (idim) matrix z to obtain the (n-7) 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),aa(n,4),q((n-7)*mm*idim) *, * right(mm*idim) integer nr(m) c ..local scalars.. real*8 co,pinv,piv,si,one integer i,irot,it,ii,i2,i3,j,jj,l,mid,nmd,m2,m3, * nrold,n4,number,n1,n7,n11,m1 c ..local arrays.. real*8 h(5),h1(5),h2(4) c ..subroutine references.. c fpgivs,fprota c .. one = 1 if(p.gt.0.) pinv = one/p n4 = n-4 n7 = n-7 n11 = n-11 mid = mm*idim m2 = m*mm m3 = n7*mm m1 = m-1 c we determine the matrix (a) and then we reduce her to c upper triangular form (r) using givens rotations. c we apply the same transformations to the rows of matrix c z to obtain the (mm) x (n-7) matrix g. c we store matrix (r) into a and aa, g into q. c the n7 x n7 upper triangular matrix (r) has the form c | a1 ' | c (r) = | ' a2 | c | 0 ' | c with (a2) a n7 x 4 matrix and (a1) a n11 x n11 upper c triangular matrix of bandwidth 5. c initialization. nmd = n7*mid do 50 i=1,nmd q(i) = 0. 50 continue do 100 i=1,n4 a(i,5) = 0. do 100 j=1,4 a(i,j) = 0. aa(i,j) = 0. 100 continue jper = 0 nrold = 0 do 760 it=1,m1 number = nr(it) 120 if(nrold.eq.number) go to 180 if(p.le.0.) go to 740 c fetch a new row of matrix (b). n1 = nrold+1 do 140 j=1,5 h(j) = b(n1,j)*pinv 140 continue c find the appropiate row of q. do 160 j=1,mid right(j) = 0. 160 continue go to 240 c fetch a new row of matrix (sp) 180 h(5) = 0. do 200 j=1,4 h(j) = sp(it,j) 200 continue c find the appropiate row of q. j = 0 do 220 ii=1,idim l = (ii-1)*m2+(it-1)*mm do 220 jj=1,mm j = j+1 l = l+1 right(j) = z(l) 220 continue c test whether there are non-zero values in the new row of (a) c corresponding to the b-splines n(j,*),j=n7+1,...,n4. 240 if(nrold.lt.n11) go to 640 if(jper.ne.0) go to 320 c initialize the matrix (aa). jk = n11+1 do 300 i=1,4 ik = jk do 260 j=1,5 if(ik.le.0) go to 280 aa(ik,i) = a(ik,j) ik = ik-1 260 continue 280 jk = jk+1 300 continue jper = 1 c if one of the non-zero elements of the new row corresponds to one of c the b-splines n(j;*),j=n7+1,...,n4,we take account of the periodicity c conditions for setting up this row of (a). 320 do 340 i=1,4 h1(i) = 0. h2(i) = 0. 340 continue h1(5) = 0. j = nrold-n11 do 420 i=1,5 j = j+1 l0 = j 360 l1 = l0-4 if(l1.le.0) go to 400 if(l1.le.n11) go to 380 l0 = l1-n11 go to 360 380 h1(l1) = h(i) go to 420 400 h2(l0) = h2(l0) + h(i) 420 continue c rotate the new row of (a) into triangle. if(n11.le.0) go to 560 c rotations with the rows 1,2,...,n11 of (a). do 540 irot=1,n11 piv = h1(1) i2 = min0(n11-irot,4) if(piv.eq.0.) go to 500 c calculate the parameters of the givens transformation. call fpgivs(piv,a(irot,1),co,si) c apply that transformation to the columns of matrix q. j = 0 do 440 ii=1,idim l = (ii-1)*m3+irot do 440 jj=1,mm j = j+1 call fprota(co,si,right(j),q(l)) l = l+n7 440 continue c apply that transformation to the rows of (a) with respect to aa. do 460 i=1,4 call fprota(co,si,h2(i),aa(irot,i)) 460 continue c apply that transformation to the rows of (a) with respect to a. if(i2.eq.0) go to 560 do 480 i=1,i2 i1 = i+1 call fprota(co,si,h1(i1),a(irot,i1)) 480 continue 500 do 520 i=1,i2 h1(i) = h1(i+1) 520 continue h1(i2+1) = 0. 540 continue c rotations with the rows n11+1,...,n7 of a. 560 do 620 irot=1,4 ij = n11+irot if(ij.le.0) go to 620 piv = h2(irot) if(piv.eq.0.) go to 620 c calculate the parameters of the givens transformation. call fpgivs(piv,aa(ij,irot),co,si) c apply that transformation to the columns of matrix q. j = 0 do 580 ii=1,idim l = (ii-1)*m3+ij do 580 jj=1,mm j = j+1 call fprota(co,si,right(j),q(l)) l = l+n7 580 continue if(irot.eq.4) go to 620 c apply that transformation to the rows of (a) with respect to aa. j1 = irot+1 do 600 i=j1,4 call fprota(co,si,h2(i),aa(ij,i)) 600 continue 620 continue go to 720 c rotation into triangle of the new row of (a), in case the elements c corresponding to the b-splines n(j;*),j=n7+1,...,n4 are all zero. 640 irot =nrold do 700 i=1,5 irot = irot+1 piv = h(i) if(piv.eq.0.) go to 700 c calculate the parameters of the givens transformation. call fpgivs(piv,a(irot,1),co,si) c apply that transformation to the columns of matrix g. j = 0 do 660 ii=1,idim l = (ii-1)*m3+irot do 660 jj=1,mm j = j+1 call fprota(co,si,right(j),q(l)) l = l+n7 660 continue c apply that transformation to the rows of (a). if(i.eq.5) go to 700 i2 = 1 i3 = i+1 do 680 j=i3,5 i2 = i2+1 call fprota(co,si,h(j),a(irot,i2)) 680 continue 700 continue 720 if(nrold.eq.number) go to 760 740 nrold = nrold+1 go to 120 760 continue return end