Line data Source code
1 0 : subroutine pda_qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
2 : integer m,n,lda,lipvt
3 : integer ipvt(lipvt)
4 : logical pivot
5 : double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
6 : c **********
7 : c
8 : c subroutine pda_qrfac
9 : c
10 : c this subroutine uses householder transformations with column
11 : c pivoting (optional) to compute a qr factorization of the
12 : c m by n matrix a. that is, pda_qrfac determines an orthogonal
13 : c matrix q, a permutation matrix p, and an upper trapezoidal
14 : c matrix r with diagonal elements of nonincreasing magnitude,
15 : c such that a*p = q*r. the householder transformation for
16 : c column k, k = 1,2,...,min(m,n), is of the form
17 : c
18 : c t
19 : c i - (1/u(k))*u*u
20 : c
21 : c where u has zeros in the first k-1 positions. the form of
22 : c this transformation and the method of pivoting first
23 : c appeared in the corresponding linpack subroutine.
24 : c
25 : c the subroutine statement is
26 : c
27 : c subroutine pda_qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
28 : c
29 : c where
30 : c
31 : c m is a positive integer input variable set to the number
32 : c of rows of a.
33 : c
34 : c n is a positive integer input variable set to the number
35 : c of columns of a.
36 : c
37 : c a is an m by n array. on input a contains the matrix for
38 : c which the qr factorization is to be computed. on output
39 : c the strict upper trapezoidal part of a contains the strict
40 : c upper trapezoidal part of r, and the lower trapezoidal
41 : c part of a contains a factored form of q (the non-trivial
42 : c elements of the u vectors described above).
43 : c
44 : c lda is a positive integer input variable not less than m
45 : c which specifies the leading dimension of the array a.
46 : c
47 : c pivot is a logical input variable. if pivot is set true,
48 : c then column pivoting is enforced. if pivot is set false,
49 : c then no column pivoting is done.
50 : c
51 : c ipvt is an integer output array of length lipvt. ipvt
52 : c defines the permutation matrix p such that a*p = q*r.
53 : c column j of p is column ipvt(j) of the identity matrix.
54 : c if pivot is false, ipvt is not referenced.
55 : c
56 : c lipvt is a positive integer input variable. if pivot is false,
57 : c then lipvt may be as small as 1. if pivot is true, then
58 : c lipvt must be at least n.
59 : c
60 : c rdiag is an output array of length n which contains the
61 : c diagonal elements of r.
62 : c
63 : c acnorm is an output array of length n which contains the
64 : c norms of the corresponding columns of the input matrix a.
65 : c if this information is not needed, then acnorm can coincide
66 : c with rdiag.
67 : c
68 : c wa is a work array of length n. if pivot is false, then wa
69 : c can coincide with rdiag.
70 : c
71 : c subprograms called
72 : c
73 : c minpack-supplied ... pda_dpmpar,pda_enorm
74 : c
75 : c fortran-supplied ... dmax1,dsqrt,min0
76 : c
77 : c argonne national laboratory. minpack project. march 1980.
78 : c burton s. garbow, kenneth e. hillstrom, jorge j. more
79 : c
80 : c **********
81 : integer i,j,jp1,k,kmax,minmn
82 : double precision ajnorm,epsmch,one,p05,sum,temp,zero
83 : double precision pda_dpmpar,pda_enorm
84 : data one,p05,zero /1.0d0,5.0d-2,0.0d0/
85 : c
86 : c epsmch is the machine precision.
87 : c
88 0 : epsmch = pda_dpmpar(1)
89 : c
90 : c compute the initial column norms and initialize several arrays.
91 : c
92 0 : do 10 j = 1, n
93 0 : acnorm(j) = pda_enorm(m,a(1,j))
94 0 : rdiag(j) = acnorm(j)
95 0 : wa(j) = rdiag(j)
96 0 : if (pivot) ipvt(j) = j
97 0 : 10 continue
98 : c
99 : c reduce a to r with householder transformations.
100 : c
101 0 : minmn = min0(m,n)
102 0 : do 110 j = 1, minmn
103 0 : if (.not.pivot) go to 40
104 : c
105 : c bring the column of largest norm into the pivot position.
106 : c
107 0 : kmax = j
108 0 : do 20 k = j, n
109 0 : if (rdiag(k) .gt. rdiag(kmax)) kmax = k
110 0 : 20 continue
111 0 : if (kmax .eq. j) go to 40
112 0 : do 30 i = 1, m
113 0 : temp = a(i,j)
114 0 : a(i,j) = a(i,kmax)
115 0 : a(i,kmax) = temp
116 0 : 30 continue
117 0 : rdiag(kmax) = rdiag(j)
118 0 : wa(kmax) = wa(j)
119 0 : k = ipvt(j)
120 0 : ipvt(j) = ipvt(kmax)
121 0 : ipvt(kmax) = k
122 : 40 continue
123 : c
124 : c compute the householder transformation to reduce the
125 : c j-th column of a to a multiple of the j-th unit vector.
126 : c
127 0 : ajnorm = pda_enorm(m-j+1,a(j,j))
128 0 : if (ajnorm .eq. zero) go to 100
129 0 : if (a(j,j) .lt. zero) ajnorm = -ajnorm
130 0 : do 50 i = j, m
131 0 : a(i,j) = a(i,j)/ajnorm
132 0 : 50 continue
133 0 : a(j,j) = a(j,j) + one
134 : c
135 : c apply the transformation to the remaining columns
136 : c and update the norms.
137 : c
138 0 : jp1 = j + 1
139 0 : if (n .lt. jp1) go to 100
140 0 : do 90 k = jp1, n
141 0 : sum = zero
142 0 : do 60 i = j, m
143 0 : sum = sum + a(i,j)*a(i,k)
144 0 : 60 continue
145 0 : temp = sum/a(j,j)
146 0 : do 70 i = j, m
147 0 : a(i,k) = a(i,k) - temp*a(i,j)
148 0 : 70 continue
149 0 : if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
150 0 : temp = a(j,k)/rdiag(k)
151 0 : rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
152 0 : if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
153 0 : rdiag(k) = pda_enorm(m-j,a(jp1,k))
154 0 : wa(k) = rdiag(k)
155 : 80 continue
156 0 : 90 continue
157 : 100 continue
158 0 : rdiag(j) = -ajnorm
159 0 : 110 continue
160 0 : return
161 : c
162 : c last card of subroutine pda_qrfac.
163 : c
164 : end
|
