www.pudn.com > TSAI30B3.rar > QRSOLV.F
subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
integer n,ldr
integer ipvt(n)
double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
c **********
c
c subroutine qrsolv
c
c given an m by n matrix a, an n by n diagonal matrix d,
c and an m-vector b, the problem is to determine an x which
c solves the system
c
c a*x = b , d*x = 0 ,
c
c in the least squares sense.
c
c this subroutine completes the solution of the problem
c if it is provided with the necessary information from the
c qr factorization, with column pivoting, of a. that is, if
c a*p = q*r, where p is a permutation matrix, q has orthogonal
c columns, and r is an upper triangular matrix with diagonal
c elements of nonincreasing magnitude, then qrsolv expects
c the full upper triangle of r, the permutation matrix p,
c and the first n components of (q transpose)*b. the system
c a*x = b, d*x = 0, is then equivalent to
c
c t t
c r*z = q *b , p *d*p*z = 0 ,
c
c where x = p*z. if this system does not have full rank,
c then a least squares solution is obtained. on output qrsolv
c also provides an upper triangular matrix s such that
c
c t t t
c p *(a *a + d*d)*p = s *s .
c
c s is computed within qrsolv and may be of separate interest.
c
c the subroutine statement is
c
c subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
c
c where
c
c n is a positive integer input variable set to the order of r.
c
c r is an n by n array. on input the full upper triangle
c must contain the full upper triangle of the matrix r.
c on output the full upper triangle is unaltered, and the
c strict lower triangle contains the strict upper triangle
c (transposed) of the upper triangular matrix s.
c
c ldr is a positive integer input variable not less than n
c which specifies the leading dimension of the array r.
c
c ipvt is an integer input array of length n which defines the
c permutation matrix p such that a*p = q*r. column j of p
c is column ipvt(j) of the identity matrix.
c
c diag is an input array of length n which must contain the
c diagonal elements of the matrix d.
c
c qtb is an input array of length n which must contain the first
c n elements of the vector (q transpose)*b.
c
c x is an output array of length n which contains the least
c squares solution of the system a*x = b, d*x = 0.
c
c sdiag is an output array of length n which contains the
c diagonal elements of the upper triangular matrix s.
c
c wa is a work array of length n.
c
c subprograms called
c
c fortran-supplied ... dabs,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jp1,k,kp1,l,nsing
double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
c
c copy r and (q transpose)*b to preserve input and initialize s.
c in particular, save the diagonal elements of r in x.
c
do 20 j = 1, n
do 10 i = j, n
r(i,j) = r(j,i)
10 continue
x(j) = r(j,j)
wa(j) = qtb(j)
20 continue
c
c eliminate the diagonal matrix d using a givens rotation.
c
do 100 j = 1, n
c
c prepare the row of d to be eliminated, locating the
c diagonal element using p from the qr factorization.
c
l = ipvt(j)
if (diag(l) .eq. zero) go to 90
do 30 k = j, n
sdiag(k) = zero
30 continue
sdiag(j) = diag(l)
c
c the transformations to eliminate the row of d
c modify only a single element of (q transpose)*b
c beyond the first n, which is initially zero.
c
qtbpj = zero
do 80 k = j, n
c
c determine a givens rotation which eliminates the
c appropriate element in the current row of d.
c
if (sdiag(k) .eq. zero) go to 70
if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
cotan = r(k,k)/sdiag(k)
sin = p5/dsqrt(p25+p25*cotan**2)
cos = sin*cotan
go to 50
40 continue
tan = sdiag(k)/r(k,k)
cos = p5/dsqrt(p25+p25*tan**2)
sin = cos*tan
50 continue
c
c compute the modified diagonal element of r and
c the modified element of ((q transpose)*b,0).
c
r(k,k) = cos*r(k,k) + sin*sdiag(k)
temp = cos*wa(k) + sin*qtbpj
qtbpj = -sin*wa(k) + cos*qtbpj
wa(k) = temp
c
c accumulate the tranformation in the row of s.
c
kp1 = k + 1
if (n .lt. kp1) go to 70
do 60 i = kp1, n
temp = cos*r(i,k) + sin*sdiag(i)
sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
r(i,k) = temp
60 continue
70 continue
80 continue
90 continue
c
c store the diagonal element of s and restore
c the corresponding diagonal element of r.
c
sdiag(j) = r(j,j)
r(j,j) = x(j)
100 continue
c
c solve the triangular system for z. if the system is
c singular, then obtain a least squares solution.
c
nsing = n
do 110 j = 1, n
if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
if (nsing .lt. n) wa(j) = zero
110 continue
if (nsing .lt. 1) go to 150
do 140 k = 1, nsing
j = nsing - k + 1
sum = zero
jp1 = j + 1
if (nsing .lt. jp1) go to 130
do 120 i = jp1, nsing
sum = sum + r(i,j)*wa(i)
120 continue
130 continue
wa(j) = (wa(j) - sum)/sdiag(j)
140 continue
150 continue
c
c permute the components of z back to components of x.
c
do 160 j = 1, n
l = ipvt(j)
x(l) = wa(j)
160 continue
return
c
c last card of subroutine qrsolv.
c
end