|
2329
|
1 SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 INTEGER IHI, ILO, INFO, LDA, N |
|
|
9 * .. |
|
|
10 * .. Array Arguments .. |
|
|
11 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) |
|
|
12 * .. |
|
|
13 * |
|
|
14 * Purpose |
|
|
15 * ======= |
|
|
16 * |
|
|
17 * DGEHD2 reduces a real general matrix A to upper Hessenberg form H by |
|
|
18 * an orthogonal similarity transformation: Q' * A * Q = H . |
|
|
19 * |
|
|
20 * Arguments |
|
|
21 * ========= |
|
|
22 * |
|
|
23 * N (input) INTEGER |
|
|
24 * The order of the matrix A. N >= 0. |
|
|
25 * |
|
|
26 * ILO (input) INTEGER |
|
|
27 * IHI (input) INTEGER |
|
|
28 * It is assumed that A is already upper triangular in rows |
|
|
29 * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally |
|
|
30 * set by a previous call to DGEBAL; otherwise they should be |
|
|
31 * set to 1 and N respectively. See Further Details. |
|
|
32 * 1 <= ILO <= IHI <= max(1,N). |
|
|
33 * |
|
|
34 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) |
|
|
35 * On entry, the n by n general matrix to be reduced. |
|
|
36 * On exit, the upper triangle and the first subdiagonal of A |
|
|
37 * are overwritten with the upper Hessenberg matrix H, and the |
|
|
38 * elements below the first subdiagonal, with the array TAU, |
|
|
39 * represent the orthogonal matrix Q as a product of elementary |
|
|
40 * reflectors. See Further Details. |
|
|
41 * |
|
|
42 * LDA (input) INTEGER |
|
|
43 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
44 * |
|
|
45 * TAU (output) DOUBLE PRECISION array, dimension (N-1) |
|
|
46 * The scalar factors of the elementary reflectors (see Further |
|
|
47 * Details). |
|
|
48 * |
|
|
49 * WORK (workspace) DOUBLE PRECISION array, dimension (N) |
|
|
50 * |
|
|
51 * INFO (output) INTEGER |
|
|
52 * = 0: successful exit. |
|
|
53 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
54 * |
|
|
55 * Further Details |
|
|
56 * =============== |
|
|
57 * |
|
|
58 * The matrix Q is represented as a product of (ihi-ilo) elementary |
|
|
59 * reflectors |
|
|
60 * |
|
|
61 * Q = H(ilo) H(ilo+1) . . . H(ihi-1). |
|
|
62 * |
|
|
63 * Each H(i) has the form |
|
|
64 * |
|
|
65 * H(i) = I - tau * v * v' |
|
|
66 * |
|
|
67 * where tau is a real scalar, and v is a real vector with |
|
|
68 * v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on |
|
|
69 * exit in A(i+2:ihi,i), and tau in TAU(i). |
|
|
70 * |
|
|
71 * The contents of A are illustrated by the following example, with |
|
|
72 * n = 7, ilo = 2 and ihi = 6: |
|
|
73 * |
|
|
74 * on entry, on exit, |
|
|
75 * |
|
|
76 * ( a a a a a a a ) ( a a h h h h a ) |
|
|
77 * ( a a a a a a ) ( a h h h h a ) |
|
|
78 * ( a a a a a a ) ( h h h h h h ) |
|
|
79 * ( a a a a a a ) ( v2 h h h h h ) |
|
|
80 * ( a a a a a a ) ( v2 v3 h h h h ) |
|
|
81 * ( a a a a a a ) ( v2 v3 v4 h h h ) |
|
|
82 * ( a ) ( a ) |
|
|
83 * |
|
|
84 * where a denotes an element of the original matrix A, h denotes a |
|
|
85 * modified element of the upper Hessenberg matrix H, and vi denotes an |
|
|
86 * element of the vector defining H(i). |
|
|
87 * |
|
|
88 * ===================================================================== |
|
|
89 * |
|
|
90 * .. Parameters .. |
|
|
91 DOUBLE PRECISION ONE |
|
|
92 PARAMETER ( ONE = 1.0D+0 ) |
|
|
93 * .. |
|
|
94 * .. Local Scalars .. |
|
|
95 INTEGER I |
|
|
96 DOUBLE PRECISION AII |
|
|
97 * .. |
|
|
98 * .. External Subroutines .. |
|
|
99 EXTERNAL DLARF, DLARFG, XERBLA |
|
|
100 * .. |
|
|
101 * .. Intrinsic Functions .. |
|
|
102 INTRINSIC MAX, MIN |
|
|
103 * .. |
|
|
104 * .. Executable Statements .. |
|
|
105 * |
|
|
106 * Test the input parameters |
|
|
107 * |
|
|
108 INFO = 0 |
|
|
109 IF( N.LT.0 ) THEN |
|
|
110 INFO = -1 |
|
|
111 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN |
|
|
112 INFO = -2 |
|
|
113 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN |
|
|
114 INFO = -3 |
|
|
115 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
116 INFO = -5 |
|
|
117 END IF |
|
|
118 IF( INFO.NE.0 ) THEN |
|
|
119 CALL XERBLA( 'DGEHD2', -INFO ) |
|
|
120 RETURN |
|
|
121 END IF |
|
|
122 * |
|
|
123 DO 10 I = ILO, IHI - 1 |
|
|
124 * |
|
|
125 * Compute elementary reflector H(i) to annihilate A(i+2:ihi,i) |
|
|
126 * |
|
|
127 CALL DLARFG( IHI-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, |
|
|
128 $ TAU( I ) ) |
|
|
129 AII = A( I+1, 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 DLARF( '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 DLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1, TAU( I ), |
|
|
140 $ A( I+1, I+1 ), LDA, WORK ) |
|
|
141 * |
|
|
142 A( I+1, I ) = AII |
|
|
143 10 CONTINUE |
|
|
144 * |
|
|
145 RETURN |
|
|
146 * |
|
|
147 * End of DGEHD2 |
|
|
148 * |
|
|
149 END |
