Line data Source code
1 0 : subroutine pda_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 pda_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 pda_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 pda_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 pda_qrsolv and may be of separate interest.
38 : c
39 : c the subroutine statement is
40 : c
41 : c subroutine pda_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 0 : do 20 j = 1, n
90 0 : do 10 i = j, n
91 0 : r(i,j) = r(j,i)
92 0 : 10 continue
93 0 : x(j) = r(j,j)
94 0 : wa(j) = qtb(j)
95 0 : 20 continue
96 : c
97 : c eliminate the diagonal matrix d using a givens rotation.
98 : c
99 0 : 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 0 : l = ipvt(j)
105 0 : if (diag(l) .eq. zero) go to 90
106 0 : do 30 k = j, n
107 0 : sdiag(k) = zero
108 0 : 30 continue
109 0 : 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 0 : qtbpj = zero
116 0 : 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 0 : if (sdiag(k) .eq. zero) go to 70
122 0 : if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
123 0 : cotan = r(k,k)/sdiag(k)
124 0 : sin = p5/dsqrt(p25+p25*cotan**2)
125 0 : cos = sin*cotan
126 0 : go to 50
127 : 40 continue
128 0 : tan = sdiag(k)/r(k,k)
129 0 : cos = p5/dsqrt(p25+p25*tan**2)
130 0 : 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 0 : r(k,k) = cos*r(k,k) + sin*sdiag(k)
137 0 : temp = cos*wa(k) + sin*qtbpj
138 0 : qtbpj = -sin*wa(k) + cos*qtbpj
139 0 : wa(k) = temp
140 : c
141 : c accumulate the tranformation in the row of s.
142 : c
143 0 : kp1 = k + 1
144 0 : if (n .lt. kp1) go to 70
145 0 : do 60 i = kp1, n
146 0 : temp = cos*r(i,k) + sin*sdiag(i)
147 0 : sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
148 0 : r(i,k) = temp
149 0 : 60 continue
150 : 70 continue
151 0 : 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 0 : sdiag(j) = r(j,j)
158 0 : r(j,j) = x(j)
159 0 : 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 0 : nsing = n
165 0 : do 110 j = 1, n
166 0 : if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
167 0 : if (nsing .lt. n) wa(j) = zero
168 0 : 110 continue
169 0 : if (nsing .lt. 1) go to 150
170 0 : do 140 k = 1, nsing
171 0 : j = nsing - k + 1
172 0 : sum = zero
173 0 : jp1 = j + 1
174 0 : if (nsing .lt. jp1) go to 130
175 0 : do 120 i = jp1, nsing
176 0 : sum = sum + r(i,j)*wa(i)
177 0 : 120 continue
178 : 130 continue
179 0 : wa(j) = (wa(j) - sum)/sdiag(j)
180 0 : 140 continue
181 : 150 continue
182 : c
183 : c permute the components of z back to components of x.
184 : c
185 0 : do 160 j = 1, n
186 0 : l = ipvt(j)
187 0 : x(l) = wa(j)
188 0 : 160 continue
189 0 : return
190 : c
191 : c last card of subroutine pda_qrsolv.
192 : c
193 : end
|
