1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
|
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
|
