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
|
subroutine qform(m,n,q,ldq,wa)
integer m,n,ldq
double precision q(ldq,m),wa(m)
c **********
c
c subroutine qform
c
c this subroutine proceeds from the computed qr factorization of
c an m by n matrix a to accumulate the m by m orthogonal matrix
c q from its factored form.
c
c the subroutine statement is
c
c subroutine qform(m,n,q,ldq,wa)
c
c where
c
c m is a positive integer input variable set to the number
c of rows of a and the order of q.
c
c n is a positive integer input variable set to the number
c of columns of a.
c
c q is an m by m array. on input the full lower trapezoid in
c the first min(m,n) columns of q contains the factored form.
c on output q has been accumulated into a square matrix.
c
c ldq is a positive integer input variable not less than m
c which specifies the leading dimension of the array q.
c
c wa is a work array of length m.
c
c subprograms called
c
c fortran-supplied ... min0
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jm1,k,l,minmn,np1
double precision one,sum,temp,zero
data one,zero /1.0d0,0.0d0/
c
c zero out upper triangle of q in the first min(m,n) columns.
c
minmn = min0(m,n)
if (minmn .lt. 2) go to 30
do 20 j = 2, minmn
jm1 = j - 1
do 10 i = 1, jm1
q(i,j) = zero
10 continue
20 continue
30 continue
c
c initialize remaining columns to those of the identity matrix.
c
np1 = n + 1
if (m .lt. np1) go to 60
do 50 j = np1, m
do 40 i = 1, m
q(i,j) = zero
40 continue
q(j,j) = one
50 continue
60 continue
c
c accumulate q from its factored form.
c
do 120 l = 1, minmn
k = minmn - l + 1
do 70 i = k, m
wa(i) = q(i,k)
q(i,k) = zero
70 continue
q(k,k) = one
if (wa(k) .eq. zero) go to 110
do 100 j = k, m
sum = zero
do 80 i = k, m
sum = sum + q(i,j)*wa(i)
80 continue
temp = sum/wa(k)
do 90 i = k, m
q(i,j) = q(i,j) - temp*wa(i)
90 continue
100 continue
110 continue
120 continue
return
c
c last card of subroutine qform.
c
end
|
