|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB) |
|
|
6 C |
|
|
7 C***BEGIN PROLOGUE DHEQR |
|
|
8 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
9 C***REVISION DATE 900926 (YYMMDD) |
|
|
10 C |
|
|
11 C----------------------------------------------------------------------- |
|
|
12 C***DESCRIPTION |
|
|
13 C |
|
|
14 C This routine performs a QR decomposition of an upper |
|
|
15 C Hessenberg matrix A. There are two options available: |
|
|
16 C |
|
|
17 C (1) performing a fresh decomposition |
|
|
18 C (2) updating the QR factors by adding a row and A |
|
|
19 C column to the matrix A. |
|
|
20 C |
|
|
21 C DHEQR decomposes an upper Hessenberg matrix by using Givens |
|
|
22 C rotations. |
|
|
23 C |
|
|
24 C On entry |
|
|
25 C |
|
|
26 C A DOUBLE PRECISION(LDA, N) |
|
|
27 C The matrix to be decomposed. |
|
|
28 C |
|
|
29 C LDA INTEGER |
|
|
30 C The leading dimension of the array A. |
|
|
31 C |
|
|
32 C N INTEGER |
|
|
33 C A is an (N+1) by N Hessenberg matrix. |
|
|
34 C |
|
|
35 C IJOB INTEGER |
|
|
36 C = 1 Means that a fresh decomposition of the |
|
|
37 C matrix A is desired. |
|
|
38 C .GE. 2 Means that the current decomposition of A |
|
|
39 C will be updated by the addition of a row |
|
|
40 C and a column. |
|
|
41 C On return |
|
|
42 C |
|
|
43 C A The upper triangular matrix R. |
|
|
44 C The factorization can be written Q*A = R, where |
|
|
45 C Q is a product of Givens rotations and R is upper |
|
|
46 C triangular. |
|
|
47 C |
|
|
48 C Q DOUBLE PRECISION(2*N) |
|
|
49 C The factors C and S of each Givens rotation used |
|
|
50 C in decomposing A. |
|
|
51 C |
|
|
52 C INFO INTEGER |
|
|
53 C = 0 normal value. |
|
|
54 C = K If A(K,K) .EQ. 0.0. This is not an error |
|
|
55 C condition for this subroutine, but it does |
|
|
56 C indicate that DHELS will divide by zero |
|
|
57 C if called. |
|
|
58 C |
|
|
59 C Modification of LINPACK. |
|
|
60 C Peter Brown, Lawrence Livermore Natl. Lab. |
|
|
61 C |
|
|
62 C----------------------------------------------------------------------- |
|
|
63 C***ROUTINES CALLED (NONE) |
|
|
64 C |
|
|
65 C***END PROLOGUE DHEQR |
|
|
66 C |
|
|
67 INTEGER LDA, N, INFO, IJOB |
|
|
68 DOUBLE PRECISION A(LDA,*), Q(*) |
|
|
69 INTEGER I, IQ, J, K, KM1, KP1, NM1 |
|
|
70 DOUBLE PRECISION C, S, T, T1, T2 |
|
|
71 C |
|
|
72 IF (IJOB .GT. 1) GO TO 70 |
|
|
73 C----------------------------------------------------------------------- |
|
|
74 C A new factorization is desired. |
|
|
75 C----------------------------------------------------------------------- |
|
|
76 C |
|
|
77 C QR decomposition without pivoting. |
|
|
78 C |
|
|
79 INFO = 0 |
|
|
80 DO 60 K = 1, N |
|
|
81 KM1 = K - 1 |
|
|
82 KP1 = K + 1 |
|
|
83 C |
|
|
84 C Compute Kth column of R. |
|
|
85 C First, multiply the Kth column of A by the previous |
|
|
86 C K-1 Givens rotations. |
|
|
87 C |
|
|
88 IF (KM1 .LT. 1) GO TO 20 |
|
|
89 DO 10 J = 1, KM1 |
|
|
90 I = 2*(J-1) + 1 |
|
|
91 T1 = A(J,K) |
|
|
92 T2 = A(J+1,K) |
|
|
93 C = Q(I) |
|
|
94 S = Q(I+1) |
|
|
95 A(J,K) = C*T1 - S*T2 |
|
|
96 A(J+1,K) = S*T1 + C*T2 |
|
|
97 10 CONTINUE |
|
|
98 C |
|
|
99 C Compute Givens components C and S. |
|
|
100 C |
|
|
101 20 CONTINUE |
|
|
102 IQ = 2*KM1 + 1 |
|
|
103 T1 = A(K,K) |
|
|
104 T2 = A(KP1,K) |
|
|
105 IF (T2 .NE. 0.0D0) GO TO 30 |
|
|
106 C = 1.0D0 |
|
|
107 S = 0.0D0 |
|
|
108 GO TO 50 |
|
|
109 30 CONTINUE |
|
|
110 IF (ABS(T2) .LT. ABS(T1)) GO TO 40 |
|
|
111 T = T1/T2 |
|
|
112 S = -1.0D0/SQRT(1.0D0+T*T) |
|
|
113 C = -S*T |
|
|
114 GO TO 50 |
|
|
115 40 CONTINUE |
|
|
116 T = T2/T1 |
|
|
117 C = 1.0D0/SQRT(1.0D0+T*T) |
|
|
118 S = -C*T |
|
|
119 50 CONTINUE |
|
|
120 Q(IQ) = C |
|
|
121 Q(IQ+1) = S |
|
|
122 A(K,K) = C*T1 - S*T2 |
|
|
123 IF (A(K,K) .EQ. 0.0D0) INFO = K |
|
|
124 60 CONTINUE |
|
|
125 RETURN |
|
|
126 C----------------------------------------------------------------------- |
|
|
127 C The old factorization of A will be updated. A row and a column |
|
|
128 C has been added to the matrix A. |
|
|
129 C N by N-1 is now the old size of the matrix. |
|
|
130 C----------------------------------------------------------------------- |
|
|
131 70 CONTINUE |
|
|
132 NM1 = N - 1 |
|
|
133 C----------------------------------------------------------------------- |
|
|
134 C Multiply the new column by the N previous Givens rotations. |
|
|
135 C----------------------------------------------------------------------- |
|
|
136 DO 100 K = 1,NM1 |
|
|
137 I = 2*(K-1) + 1 |
|
|
138 T1 = A(K,N) |
|
|
139 T2 = A(K+1,N) |
|
|
140 C = Q(I) |
|
|
141 S = Q(I+1) |
|
|
142 A(K,N) = C*T1 - S*T2 |
|
|
143 A(K+1,N) = S*T1 + C*T2 |
|
|
144 100 CONTINUE |
|
|
145 C----------------------------------------------------------------------- |
|
|
146 C Complete update of decomposition by forming last Givens rotation, |
|
|
147 C and multiplying it times the column vector (A(N,N),A(NP1,N)). |
|
|
148 C----------------------------------------------------------------------- |
|
|
149 INFO = 0 |
|
|
150 T1 = A(N,N) |
|
|
151 T2 = A(N+1,N) |
|
|
152 IF (T2 .NE. 0.0D0) GO TO 110 |
|
|
153 C = 1.0D0 |
|
|
154 S = 0.0D0 |
|
|
155 GO TO 130 |
|
|
156 110 CONTINUE |
|
|
157 IF (ABS(T2) .LT. ABS(T1)) GO TO 120 |
|
|
158 T = T1/T2 |
|
|
159 S = -1.0D0/SQRT(1.0D0+T*T) |
|
|
160 C = -S*T |
|
|
161 GO TO 130 |
|
|
162 120 CONTINUE |
|
|
163 T = T2/T1 |
|
|
164 C = 1.0D0/SQRT(1.0D0+T*T) |
|
|
165 S = -C*T |
|
|
166 130 CONTINUE |
|
|
167 IQ = 2*N - 1 |
|
|
168 Q(IQ) = C |
|
|
169 Q(IQ+1) = S |
|
|
170 A(N,N) = C*T1 - S*T2 |
|
|
171 IF (A(N,N) .EQ. 0.0D0) INFO = N |
|
|
172 RETURN |
|
|
173 C |
|
|
174 C------END OF SUBROUTINE DHEQR------------------------------------------ |
|
|
175 END |
