V3FIT
qrfac.f
1  SUBROUTINE qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
2  USE stel_kinds
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D u m m y A r g u m e n t s
6 C-----------------------------------------------
7  INTEGER, INTENT(in) :: m, n, lda, lipvt
8  LOGICAL, INTENT(in) :: pivot
9  INTEGER, DIMENSION(lipvt), INTENT(out) :: ipvt
10  REAL(rprec), DIMENSION(lda,n), INTENT(inout) :: a
11  REAL(rprec), DIMENSION(n), INTENT(out) :: rdiag, acnorm, wa
12 C-----------------------------------------------
13 C L o c a l P a r a m e t e r s
14 C-----------------------------------------------
15  REAL(rprec), PARAMETER :: one = 1, p05 = 0.05_dp,
16  1 zero = 0
17 C-----------------------------------------------
18 C L o c a l V a r i a b l e s
19 C-----------------------------------------------
20  INTEGER :: j, jp1, k, kmax, minmn
21  REAL(rprec) :: ajnorm, epsmch, sum0, temp, enorm, dpmpar
22  REAL(rprec) , ALLOCATABLE :: temp1u(:)
23 C-----------------------------------------------
24 C E x t e r n a l F u n c t i o n s
25 C-----------------------------------------------
26  EXTERNAL dpmpar,enorm
27 C-----------------------------------------------
28 c
29 c SUBROUTINE qrfac
30 c
31 c this SUBROUTINE uses householder transformations with column
32 c pivoting (optional) to compute a qr factorization of the
33 c m by n matrix a. that is, qrfac determines an orthogonal
34 c matrix q, a permutation matrix p, and an upper trapezoidal
35 c matrix r with diagonal elements of nonincreasing magnitude,
36 c such that a*p = q*r. the householder transformation for
37 c column k, k = 1,2,...,MIN(m,n), is of the form
38 c
39 c t
40 c i - (1/u(k))*u*u
41 c
42 c WHERE u has zeros in the first k-1 positions. the form of
43 c this transformation and the method of pivoting first
44 c appeared in the corresponding linpack SUBROUTINE.
45 c
46 c the SUBROUTINE statement is
47 c
48 c SUBROUTINE qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
49 c
50 c WHERE
51 c
52 c m is a positive INTEGER input variable set to the number
53 c of rows of a.
54 c
55 c n is a positive INTEGER input variable set to the number
56 c of columns of a.
57 c
58 c a is an m by n array. on input a CONTAINS the matrix for
59 c which the qr factorization is to be computed. on output
60 c the strict upper trapezoidal part of a CONTAINS the strict
61 c upper trapezoidal part of r, and the lower trapezoidal
62 c part of a CONTAINS a factored form of q (the non-trivial
63 c elements of the u vectors described above).
64 c
65 c lda is a positive INTEGER input variable not less than m
66 c which specifies the leading DIMENSION of the array a.
67 c
68 c pivot is a LOGICAL input variable. IF pivot is set true,
69 c THEN column pivoting is enforced. IF pivot is set false,
70 c THEN no column pivoting is done.
71 c
72 c ipvt is an INTEGER output array of length lipvt. ipvt
73 c defines the permutation matrix p such that a*p = q*r.
74 c column j of p is column ipvt(j) of the identity matrix.
75 c IF pivot is false, ipvt is not referenced.
76 c
77 c lipvt is a positive integer input variable. if pivot is false,
78 c then lipvt may be as small as 1. if pivot is true, then
79 c lipvt must be at least n.
80 c
81 c rdiag is an output array of length n which CONTAINS the
82 c diagonal elements of r.
83 c
84 c acnorm is an output array of length n which CONTAINS the
85 c norms of the corresponding columns of the input matrix a.
86 c IF this information is not needed, THEN acnorm can coincide
87 c with rdiag.
88 c
89 c wa is a work array of length n. IF pivot is false, THEN wa
90 c can coincide with rdiag.
91 c
92 c subprograms called
93 c
94 c MINpack-supplied ... dpmpar,enorm
95 c
96 c fortran-supplied ... MAX,sqrt,min
97 c
98 c argonne national laboratory. MINpack project. march 1980.
99 c burton s. garbow, kenneth e. hillstrom, jorge j. more
100 c
101 c **********
102 c
103 c epsmch is the machine precision.
104 c
105  epsmch = dpmpar(1)
106 c
107 c compute the initial column norms and initialize several arrays.
108 c
109  DO j = 1, n
110  acnorm(j) = enorm(m,a(1,j))
111  rdiag(j) = acnorm(j)
112  wa(j) = rdiag(j)
113  IF (pivot) ipvt(j) = j
114  END DO
115 c
116 c reduce a to r with householder transformations.
117 c
118  minmn = min(m,n)
119  ALLOCATE (temp1u(m))
120 
121  DO j = 1, minmn
122  IF (pivot) THEN
123 c
124 c bring the column of largest norm into the pivot position.
125 c
126  kmax = j
127  DO k = j,n
128  IF (rdiag(k) .gt. rdiag(kmax)) kmax = k
129  END DO
130  IF (kmax .ne. j) THEN
131  temp1u(:m) = a(:m,j)
132  a(:m,j) = a(:m,kmax)
133  a(:m,kmax) = temp1u(:m)
134  rdiag(kmax) = rdiag(j)
135  wa(kmax) = wa(j)
136  k = ipvt(j)
137  ipvt(j) = ipvt(kmax)
138  ipvt(kmax) = k
139  END IF
140  END IF
141 c
142 c compute the householder transformation to reduce the
143 c j-th column of a to a multiple of the j-th unit vector.
144 c
145  ajnorm = enorm(m-j+1,a(j,j))
146  IF (ajnorm .ne. zero) THEN
147  IF (a(j,j) .lt. zero) ajnorm = -ajnorm
148  a(j:m,j) = a(j:m,j)/ajnorm
149  a(j,j) = a(j,j) + one
150 c
151 c apply the transformation to the remaining columns
152 c and update the norms.
153 c
154  jp1 = j + 1
155  IF (n .ge. jp1) THEN
156  DO k = jp1, n
157  sum0 = sum(a(j:m,j)*a(j:m,k))
158  temp = sum0/a(j,j)
159  a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
160  IF (pivot .and. rdiag(k).ne.zero) THEN
161  temp = a(j,k)/rdiag(k)
162  rdiag(k) = rdiag(k)*sqrt(max(zero,one-temp**2))
163  IF (p05*(rdiag(k)/wa(k))**2 .le. epsmch) THEN
164  rdiag(k) = enorm(m - j,a(jp1,k))
165  wa(k) = rdiag(k)
166  END IF
167  END IF
168  END DO
169  END IF
170  END IF
171  rdiag(j) = -ajnorm
172  END DO
173 
174  DEALLOCATE (temp1u)
175 
176  END SUBROUTINE qrfac
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17