GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dheqr.f
Go to the documentation of this file.
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
dheqr
subroutine dheqr(A, LDA, N, Q, INFO, IJOB)
Definition: dheqr.f:6
abs
static T abs(T x)
Definition: pr-output.cc:1678