Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/lapack/cgehd2.f
5204 views
1
SUBROUTINE CGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
2
*
3
* -- LAPACK routine (version 3.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
6
* September 30, 1994
7
*
8
* .. Scalar Arguments ..
9
INTEGER IHI, ILO, INFO, LDA, N
10
* ..
11
* .. Array Arguments ..
12
COMPLEX A( LDA, * ), TAU( * ), WORK( * )
13
* ..
14
*
15
* Purpose
16
* =======
17
*
18
* CGEHD2 reduces a complex general matrix A to upper Hessenberg form H
19
* by a unitary similarity transformation: Q' * A * Q = H .
20
*
21
* Arguments
22
* =========
23
*
24
* N (input) INTEGER
25
* The order of the matrix A. N >= 0.
26
*
27
* ILO (input) INTEGER
28
* IHI (input) INTEGER
29
* It is assumed that A is already upper triangular in rows
30
* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
31
* set by a previous call to CGEBAL; otherwise they should be
32
* set to 1 and N respectively. See Further Details.
33
* 1 <= ILO <= IHI <= max(1,N).
34
*
35
* A (input/output) COMPLEX array, dimension (LDA,N)
36
* On entry, the n by n general matrix to be reduced.
37
* On exit, the upper triangle and the first subdiagonal of A
38
* are overwritten with the upper Hessenberg matrix H, and the
39
* elements below the first subdiagonal, with the array TAU,
40
* represent the unitary matrix Q as a product of elementary
41
* reflectors. See Further Details.
42
*
43
* LDA (input) INTEGER
44
* The leading dimension of the array A. LDA >= max(1,N).
45
*
46
* TAU (output) COMPLEX array, dimension (N-1)
47
* The scalar factors of the elementary reflectors (see Further
48
* Details).
49
*
50
* WORK (workspace) COMPLEX array, dimension (N)
51
*
52
* INFO (output) INTEGER
53
* = 0: successful exit
54
* < 0: if INFO = -i, the i-th argument had an illegal value.
55
*
56
* Further Details
57
* ===============
58
*
59
* The matrix Q is represented as a product of (ihi-ilo) elementary
60
* reflectors
61
*
62
* Q = H(ilo) H(ilo+1) . . . H(ihi-1).
63
*
64
* Each H(i) has the form
65
*
66
* H(i) = I - tau * v * v'
67
*
68
* where tau is a complex scalar, and v is a complex vector with
69
* v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
70
* exit in A(i+2:ihi,i), and tau in TAU(i).
71
*
72
* The contents of A are illustrated by the following example, with
73
* n = 7, ilo = 2 and ihi = 6:
74
*
75
* on entry, on exit,
76
*
77
* ( a a a a a a a ) ( a a h h h h a )
78
* ( a a a a a a ) ( a h h h h a )
79
* ( a a a a a a ) ( h h h h h h )
80
* ( a a a a a a ) ( v2 h h h h h )
81
* ( a a a a a a ) ( v2 v3 h h h h )
82
* ( a a a a a a ) ( v2 v3 v4 h h h )
83
* ( a ) ( a )
84
*
85
* where a denotes an element of the original matrix A, h denotes a
86
* modified element of the upper Hessenberg matrix H, and vi denotes an
87
* element of the vector defining H(i).
88
*
89
* =====================================================================
90
*
91
* .. Parameters ..
92
COMPLEX ONE
93
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
94
* ..
95
* .. Local Scalars ..
96
INTEGER I
97
COMPLEX ALPHA
98
* ..
99
* .. External Subroutines ..
100
EXTERNAL CLARF, CLARFG, XERBLA
101
* ..
102
* .. Intrinsic Functions ..
103
INTRINSIC CONJG, MAX, MIN
104
* ..
105
* .. Executable Statements ..
106
*
107
* Test the input parameters
108
*
109
INFO = 0
110
IF( N.LT.0 ) THEN
111
INFO = -1
112
ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
113
INFO = -2
114
ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
115
INFO = -3
116
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
117
INFO = -5
118
END IF
119
IF( INFO.NE.0 ) THEN
120
CALL XERBLA( 'CGEHD2', -INFO )
121
RETURN
122
END IF
123
*
124
DO 10 I = ILO, IHI - 1
125
*
126
* Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
127
*
128
ALPHA = A( I+1, I )
129
CALL CLARFG( IHI-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAU( I ) )
130
A( I+1, I ) = ONE
131
*
132
* Apply H(i) to A(1:ihi,i+1:ihi) from the right
133
*
134
CALL CLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ),
135
$ A( 1, I+1 ), LDA, WORK )
136
*
137
* Apply H(i)' to A(i+1:ihi,i+1:n) from the left
138
*
139
CALL CLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1,
140
$ CONJG( TAU( I ) ), A( I+1, I+1 ), LDA, WORK )
141
*
142
A( I+1, I ) = ALPHA
143
10 CONTINUE
144
*
145
RETURN
146
*
147
* End of CGEHD2
148
*
149
END
150
151