GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
subroutine dheqr(A, LDA, N, Q, INFO, IJOB)
Definition: dheqr.f:5
T abs(T x)
Definition: pr-output.cc:3062
octave_value sqrt(void) const
Definition: ov.h:1200