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
|
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
|
