Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 194 lines (193 sloc) 6.178 kb
f865492 @certik Add current minpack from netlib.org
authored
1 subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
2 integer n,ldr
3 integer ipvt(n)
4 double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
5 c **********
6 c
7 c subroutine qrsolv
8 c
9 c given an m by n matrix a, an n by n diagonal matrix d,
10 c and an m-vector b, the problem is to determine an x which
11 c solves the system
12 c
13 c a*x = b , d*x = 0 ,
14 c
15 c in the least squares sense.
16 c
17 c this subroutine completes the solution of the problem
18 c if it is provided with the necessary information from the
19 c qr factorization, with column pivoting, of a. that is, if
20 c a*p = q*r, where p is a permutation matrix, q has orthogonal
21 c columns, and r is an upper triangular matrix with diagonal
22 c elements of nonincreasing magnitude, then qrsolv expects
23 c the full upper triangle of r, the permutation matrix p,
24 c and the first n components of (q transpose)*b. the system
25 c a*x = b, d*x = 0, is then equivalent to
26 c
27 c t t
28 c r*z = q *b , p *d*p*z = 0 ,
29 c
30 c where x = p*z. if this system does not have full rank,
31 c then a least squares solution is obtained. on output qrsolv
32 c also provides an upper triangular matrix s such that
33 c
34 c t t t
35 c p *(a *a + d*d)*p = s *s .
36 c
37 c s is computed within qrsolv and may be of separate interest.
38 c
39 c the subroutine statement is
40 c
41 c subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
42 c
43 c where
44 c
45 c n is a positive integer input variable set to the order of r.
46 c
47 c r is an n by n array. on input the full upper triangle
48 c must contain the full upper triangle of the matrix r.
49 c on output the full upper triangle is unaltered, and the
50 c strict lower triangle contains the strict upper triangle
51 c (transposed) of the upper triangular matrix s.
52 c
53 c ldr is a positive integer input variable not less than n
54 c which specifies the leading dimension of the array r.
55 c
56 c ipvt is an integer input array of length n which defines the
57 c permutation matrix p such that a*p = q*r. column j of p
58 c is column ipvt(j) of the identity matrix.
59 c
60 c diag is an input array of length n which must contain the
61 c diagonal elements of the matrix d.
62 c
63 c qtb is an input array of length n which must contain the first
64 c n elements of the vector (q transpose)*b.
65 c
66 c x is an output array of length n which contains the least
67 c squares solution of the system a*x = b, d*x = 0.
68 c
69 c sdiag is an output array of length n which contains the
70 c diagonal elements of the upper triangular matrix s.
71 c
72 c wa is a work array of length n.
73 c
74 c subprograms called
75 c
76 c fortran-supplied ... dabs,dsqrt
77 c
78 c argonne national laboratory. minpack project. march 1980.
79 c burton s. garbow, kenneth e. hillstrom, jorge j. more
80 c
81 c **********
82 integer i,j,jp1,k,kp1,l,nsing
83 double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
84 data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
85 c
86 c copy r and (q transpose)*b to preserve input and initialize s.
87 c in particular, save the diagonal elements of r in x.
88 c
89 do 20 j = 1, n
90 do 10 i = j, n
91 r(i,j) = r(j,i)
92 10 continue
93 x(j) = r(j,j)
94 wa(j) = qtb(j)
95 20 continue
96 c
97 c eliminate the diagonal matrix d using a givens rotation.
98 c
99 do 100 j = 1, n
100 c
101 c prepare the row of d to be eliminated, locating the
102 c diagonal element using p from the qr factorization.
103 c
104 l = ipvt(j)
105 if (diag(l) .eq. zero) go to 90
106 do 30 k = j, n
107 sdiag(k) = zero
108 30 continue
109 sdiag(j) = diag(l)
110 c
111 c the transformations to eliminate the row of d
112 c modify only a single element of (q transpose)*b
113 c beyond the first n, which is initially zero.
114 c
115 qtbpj = zero
116 do 80 k = j, n
117 c
118 c determine a givens rotation which eliminates the
119 c appropriate element in the current row of d.
120 c
121 if (sdiag(k) .eq. zero) go to 70
122 if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
123 cotan = r(k,k)/sdiag(k)
124 sin = p5/dsqrt(p25+p25*cotan**2)
125 cos = sin*cotan
126 go to 50
127 40 continue
128 tan = sdiag(k)/r(k,k)
129 cos = p5/dsqrt(p25+p25*tan**2)
130 sin = cos*tan
131 50 continue
132 c
133 c compute the modified diagonal element of r and
134 c the modified element of ((q transpose)*b,0).
135 c
136 r(k,k) = cos*r(k,k) + sin*sdiag(k)
137 temp = cos*wa(k) + sin*qtbpj
138 qtbpj = -sin*wa(k) + cos*qtbpj
139 wa(k) = temp
140 c
141 c accumulate the tranformation in the row of s.
142 c
143 kp1 = k + 1
144 if (n .lt. kp1) go to 70
145 do 60 i = kp1, n
146 temp = cos*r(i,k) + sin*sdiag(i)
147 sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
148 r(i,k) = temp
149 60 continue
150 70 continue
151 80 continue
152 90 continue
153 c
154 c store the diagonal element of s and restore
155 c the corresponding diagonal element of r.
156 c
157 sdiag(j) = r(j,j)
158 r(j,j) = x(j)
159 100 continue
160 c
161 c solve the triangular system for z. if the system is
162 c singular, then obtain a least squares solution.
163 c
164 nsing = n
165 do 110 j = 1, n
166 if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
167 if (nsing .lt. n) wa(j) = zero
168 110 continue
169 if (nsing .lt. 1) go to 150
170 do 140 k = 1, nsing
171 j = nsing - k + 1
172 sum = zero
173 jp1 = j + 1
174 if (nsing .lt. jp1) go to 130
175 do 120 i = jp1, nsing
176 sum = sum + r(i,j)*wa(i)
177 120 continue
178 130 continue
179 wa(j) = (wa(j) - sum)/sdiag(j)
180 140 continue
181 150 continue
182 c
183 c permute the components of z back to components of x.
184 c
185 do 160 j = 1, n
186 l = ipvt(j)
187 x(l) = wa(j)
188 160 continue
189 return
190 c
191 c last card of subroutine qrsolv.
192 c
193 end
Something went wrong with that request. Please try again.