1 SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
2 $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
3 $ LWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBVSL, JOBVSR
12 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
16 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
17 $ VSR( LDVSR, * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * This routine is deprecated and has been replaced by routine DGGES.
24 *
25 * DGEGS computes the eigenvalues, real Schur form, and, optionally,
26 * left and or/right Schur vectors of a real matrix pair (A,B).
27 * Given two square matrices A and B, the generalized real Schur
28 * factorization has the form
29 *
30 * A = Q*S*Z**T, B = Q*T*Z**T
31 *
32 * where Q and Z are orthogonal matrices, T is upper triangular, and S
33 * is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
34 * blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
35 * of eigenvalues of (A,B). The columns of Q are the left Schur vectors
36 * and the columns of Z are the right Schur vectors.
37 *
38 * If only the eigenvalues of (A,B) are needed, the driver routine
39 * DGEGV should be used instead. See DGEGV for a description of the
40 * eigenvalues of the generalized nonsymmetric eigenvalue problem
41 * (GNEP).
42 *
43 * Arguments
44 * =========
45 *
46 * JOBVSL (input) CHARACTER*1
47 * = 'N': do not compute the left Schur vectors;
48 * = 'V': compute the left Schur vectors (returned in VSL).
49 *
50 * JOBVSR (input) CHARACTER*1
51 * = 'N': do not compute the right Schur vectors;
52 * = 'V': compute the right Schur vectors (returned in VSR).
53 *
54 * N (input) INTEGER
55 * The order of the matrices A, B, VSL, and VSR. N >= 0.
56 *
57 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
58 * On entry, the matrix A.
59 * On exit, the upper quasi-triangular matrix S from the
60 * generalized real Schur factorization.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of A. LDA >= max(1,N).
64 *
65 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
66 * On entry, the matrix B.
67 * On exit, the upper triangular matrix T from the generalized
68 * real Schur factorization.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of B. LDB >= max(1,N).
72 *
73 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
74 * The real parts of each scalar alpha defining an eigenvalue
75 * of GNEP.
76 *
77 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
78 * The imaginary parts of each scalar alpha defining an
79 * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
80 * eigenvalue is real; if positive, then the j-th and (j+1)-st
81 * eigenvalues are a complex conjugate pair, with
82 * ALPHAI(j+1) = -ALPHAI(j).
83 *
84 * BETA (output) DOUBLE PRECISION array, dimension (N)
85 * The scalars beta that define the eigenvalues of GNEP.
86 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
87 * beta = BETA(j) represent the j-th eigenvalue of the matrix
88 * pair (A,B), in one of the forms lambda = alpha/beta or
89 * mu = beta/alpha. Since either lambda or mu may overflow,
90 * they should not, in general, be computed.
91 *
92 * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
93 * If JOBVSL = 'V', the matrix of left Schur vectors Q.
94 * Not referenced if JOBVSL = 'N'.
95 *
96 * LDVSL (input) INTEGER
97 * The leading dimension of the matrix VSL. LDVSL >=1, and
98 * if JOBVSL = 'V', LDVSL >= N.
99 *
100 * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
101 * If JOBVSR = 'V', the matrix of right Schur vectors Z.
102 * Not referenced if JOBVSR = 'N'.
103 *
104 * LDVSR (input) INTEGER
105 * The leading dimension of the matrix VSR. LDVSR >= 1, and
106 * if JOBVSR = 'V', LDVSR >= N.
107 *
108 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
109 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
110 *
111 * LWORK (input) INTEGER
112 * The dimension of the array WORK. LWORK >= max(1,4*N).
113 * For good performance, LWORK must generally be larger.
114 * To compute the optimal value of LWORK, call ILAENV to get
115 * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
116 * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
117 * The optimal LWORK is 2*N + N*(NB+1).
118 *
119 * If LWORK = -1, then a workspace query is assumed; the routine
120 * only calculates the optimal size of the WORK array, returns
121 * this value as the first entry of the WORK array, and no error
122 * message related to LWORK is issued by XERBLA.
123 *
124 * INFO (output) INTEGER
125 * = 0: successful exit
126 * < 0: if INFO = -i, the i-th argument had an illegal value.
127 * = 1,...,N:
128 * The QZ iteration failed. (A,B) are not in Schur
129 * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
130 * be correct for j=INFO+1,...,N.
131 * > N: errors that usually indicate LAPACK problems:
132 * =N+1: error return from DGGBAL
133 * =N+2: error return from DGEQRF
134 * =N+3: error return from DORMQR
135 * =N+4: error return from DORGQR
136 * =N+5: error return from DGGHRD
137 * =N+6: error return from DHGEQZ (other than failed
138 * iteration)
139 * =N+7: error return from DGGBAK (computing VSL)
140 * =N+8: error return from DGGBAK (computing VSR)
141 * =N+9: error return from DLASCL (various places)
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146 DOUBLE PRECISION ZERO, ONE
147 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
148 * ..
149 * .. Local Scalars ..
150 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
151 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
152 $ IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
153 $ LWKOPT, NB, NB1, NB2, NB3
154 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
155 $ SAFMIN, SMLNUM
156 * ..
157 * .. External Subroutines ..
158 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
159 $ DLASCL, DLASET, DORGQR, DORMQR, XERBLA
160 * ..
161 * .. External Functions ..
162 LOGICAL LSAME
163 INTEGER ILAENV
164 DOUBLE PRECISION DLAMCH, DLANGE
165 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
166 * ..
167 * .. Intrinsic Functions ..
168 INTRINSIC INT, MAX
169 * ..
170 * .. Executable Statements ..
171 *
172 * Decode the input arguments
173 *
174 IF( LSAME( JOBVSL, 'N' ) ) THEN
175 IJOBVL = 1
176 ILVSL = .FALSE.
177 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
178 IJOBVL = 2
179 ILVSL = .TRUE.
180 ELSE
181 IJOBVL = -1
182 ILVSL = .FALSE.
183 END IF
184 *
185 IF( LSAME( JOBVSR, 'N' ) ) THEN
186 IJOBVR = 1
187 ILVSR = .FALSE.
188 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
189 IJOBVR = 2
190 ILVSR = .TRUE.
191 ELSE
192 IJOBVR = -1
193 ILVSR = .FALSE.
194 END IF
195 *
196 * Test the input arguments
197 *
198 LWKMIN = MAX( 4*N, 1 )
199 LWKOPT = LWKMIN
200 WORK( 1 ) = LWKOPT
201 LQUERY = ( LWORK.EQ.-1 )
202 INFO = 0
203 IF( IJOBVL.LE.0 ) THEN
204 INFO = -1
205 ELSE IF( IJOBVR.LE.0 ) THEN
206 INFO = -2
207 ELSE IF( N.LT.0 ) THEN
208 INFO = -3
209 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
210 INFO = -5
211 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
212 INFO = -7
213 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
214 INFO = -12
215 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
216 INFO = -14
217 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
218 INFO = -16
219 END IF
220 *
221 IF( INFO.EQ.0 ) THEN
222 NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
223 NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
224 NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
225 NB = MAX( NB1, NB2, NB3 )
226 LOPT = 2*N + N*( NB+1 )
227 WORK( 1 ) = LOPT
228 END IF
229 *
230 IF( INFO.NE.0 ) THEN
231 CALL XERBLA( 'DGEGS ', -INFO )
232 RETURN
233 ELSE IF( LQUERY ) THEN
234 RETURN
235 END IF
236 *
237 * Quick return if possible
238 *
239 IF( N.EQ.0 )
240 $ RETURN
241 *
242 * Get machine constants
243 *
244 EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
245 SAFMIN = DLAMCH( 'S' )
246 SMLNUM = N*SAFMIN / EPS
247 BIGNUM = ONE / SMLNUM
248 *
249 * Scale A if max element outside range [SMLNUM,BIGNUM]
250 *
251 ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
252 ILASCL = .FALSE.
253 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
254 ANRMTO = SMLNUM
255 ILASCL = .TRUE.
256 ELSE IF( ANRM.GT.BIGNUM ) THEN
257 ANRMTO = BIGNUM
258 ILASCL = .TRUE.
259 END IF
260 *
261 IF( ILASCL ) THEN
262 CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
263 IF( IINFO.NE.0 ) THEN
264 INFO = N + 9
265 RETURN
266 END IF
267 END IF
268 *
269 * Scale B if max element outside range [SMLNUM,BIGNUM]
270 *
271 BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
272 ILBSCL = .FALSE.
273 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
274 BNRMTO = SMLNUM
275 ILBSCL = .TRUE.
276 ELSE IF( BNRM.GT.BIGNUM ) THEN
277 BNRMTO = BIGNUM
278 ILBSCL = .TRUE.
279 END IF
280 *
281 IF( ILBSCL ) THEN
282 CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
283 IF( IINFO.NE.0 ) THEN
284 INFO = N + 9
285 RETURN
286 END IF
287 END IF
288 *
289 * Permute the matrix to make it more nearly triangular
290 * Workspace layout: (2*N words -- "work..." not actually used)
291 * left_permutation, right_permutation, work...
292 *
293 ILEFT = 1
294 IRIGHT = N + 1
295 IWORK = IRIGHT + N
296 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
297 $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
298 IF( IINFO.NE.0 ) THEN
299 INFO = N + 1
300 GO TO 10
301 END IF
302 *
303 * Reduce B to triangular form, and initialize VSL and/or VSR
304 * Workspace layout: ("work..." must have at least N words)
305 * left_permutation, right_permutation, tau, work...
306 *
307 IROWS = IHI + 1 - ILO
308 ICOLS = N + 1 - ILO
309 ITAU = IWORK
310 IWORK = ITAU + IROWS
311 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
312 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
313 IF( IINFO.GE.0 )
314 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
315 IF( IINFO.NE.0 ) THEN
316 INFO = N + 2
317 GO TO 10
318 END IF
319 *
320 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
321 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
322 $ LWORK+1-IWORK, IINFO )
323 IF( IINFO.GE.0 )
324 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
325 IF( IINFO.NE.0 ) THEN
326 INFO = N + 3
327 GO TO 10
328 END IF
329 *
330 IF( ILVSL ) THEN
331 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
332 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
333 $ VSL( ILO+1, ILO ), LDVSL )
334 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
335 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
336 $ IINFO )
337 IF( IINFO.GE.0 )
338 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
339 IF( IINFO.NE.0 ) THEN
340 INFO = N + 4
341 GO TO 10
342 END IF
343 END IF
344 *
345 IF( ILVSR )
346 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
347 *
348 * Reduce to generalized Hessenberg form
349 *
350 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
351 $ LDVSL, VSR, LDVSR, IINFO )
352 IF( IINFO.NE.0 ) THEN
353 INFO = N + 5
354 GO TO 10
355 END IF
356 *
357 * Perform QZ algorithm, computing Schur vectors if desired
358 * Workspace layout: ("work..." must have at least 1 word)
359 * left_permutation, right_permutation, work...
360 *
361 IWORK = ITAU
362 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
363 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
364 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
365 IF( IINFO.GE.0 )
366 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
367 IF( IINFO.NE.0 ) THEN
368 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
369 INFO = IINFO
370 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
371 INFO = IINFO - N
372 ELSE
373 INFO = N + 6
374 END IF
375 GO TO 10
376 END IF
377 *
378 * Apply permutation to VSL and VSR
379 *
380 IF( ILVSL ) THEN
381 CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
382 $ WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
383 IF( IINFO.NE.0 ) THEN
384 INFO = N + 7
385 GO TO 10
386 END IF
387 END IF
388 IF( ILVSR ) THEN
389 CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
390 $ WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
391 IF( IINFO.NE.0 ) THEN
392 INFO = N + 8
393 GO TO 10
394 END IF
395 END IF
396 *
397 * Undo scaling
398 *
399 IF( ILASCL ) THEN
400 CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
401 IF( IINFO.NE.0 ) THEN
402 INFO = N + 9
403 RETURN
404 END IF
405 CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
406 $ IINFO )
407 IF( IINFO.NE.0 ) THEN
408 INFO = N + 9
409 RETURN
410 END IF
411 CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
412 $ IINFO )
413 IF( IINFO.NE.0 ) THEN
414 INFO = N + 9
415 RETURN
416 END IF
417 END IF
418 *
419 IF( ILBSCL ) THEN
420 CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
421 IF( IINFO.NE.0 ) THEN
422 INFO = N + 9
423 RETURN
424 END IF
425 CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
426 IF( IINFO.NE.0 ) THEN
427 INFO = N + 9
428 RETURN
429 END IF
430 END IF
431 *
432 10 CONTINUE
433 WORK( 1 ) = LWKOPT
434 *
435 RETURN
436 *
437 * End of DGEGS
438 *
439 END
2 $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
3 $ LWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBVSL, JOBVSR
12 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
16 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
17 $ VSR( LDVSR, * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * This routine is deprecated and has been replaced by routine DGGES.
24 *
25 * DGEGS computes the eigenvalues, real Schur form, and, optionally,
26 * left and or/right Schur vectors of a real matrix pair (A,B).
27 * Given two square matrices A and B, the generalized real Schur
28 * factorization has the form
29 *
30 * A = Q*S*Z**T, B = Q*T*Z**T
31 *
32 * where Q and Z are orthogonal matrices, T is upper triangular, and S
33 * is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
34 * blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
35 * of eigenvalues of (A,B). The columns of Q are the left Schur vectors
36 * and the columns of Z are the right Schur vectors.
37 *
38 * If only the eigenvalues of (A,B) are needed, the driver routine
39 * DGEGV should be used instead. See DGEGV for a description of the
40 * eigenvalues of the generalized nonsymmetric eigenvalue problem
41 * (GNEP).
42 *
43 * Arguments
44 * =========
45 *
46 * JOBVSL (input) CHARACTER*1
47 * = 'N': do not compute the left Schur vectors;
48 * = 'V': compute the left Schur vectors (returned in VSL).
49 *
50 * JOBVSR (input) CHARACTER*1
51 * = 'N': do not compute the right Schur vectors;
52 * = 'V': compute the right Schur vectors (returned in VSR).
53 *
54 * N (input) INTEGER
55 * The order of the matrices A, B, VSL, and VSR. N >= 0.
56 *
57 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
58 * On entry, the matrix A.
59 * On exit, the upper quasi-triangular matrix S from the
60 * generalized real Schur factorization.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of A. LDA >= max(1,N).
64 *
65 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
66 * On entry, the matrix B.
67 * On exit, the upper triangular matrix T from the generalized
68 * real Schur factorization.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of B. LDB >= max(1,N).
72 *
73 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
74 * The real parts of each scalar alpha defining an eigenvalue
75 * of GNEP.
76 *
77 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
78 * The imaginary parts of each scalar alpha defining an
79 * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
80 * eigenvalue is real; if positive, then the j-th and (j+1)-st
81 * eigenvalues are a complex conjugate pair, with
82 * ALPHAI(j+1) = -ALPHAI(j).
83 *
84 * BETA (output) DOUBLE PRECISION array, dimension (N)
85 * The scalars beta that define the eigenvalues of GNEP.
86 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
87 * beta = BETA(j) represent the j-th eigenvalue of the matrix
88 * pair (A,B), in one of the forms lambda = alpha/beta or
89 * mu = beta/alpha. Since either lambda or mu may overflow,
90 * they should not, in general, be computed.
91 *
92 * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
93 * If JOBVSL = 'V', the matrix of left Schur vectors Q.
94 * Not referenced if JOBVSL = 'N'.
95 *
96 * LDVSL (input) INTEGER
97 * The leading dimension of the matrix VSL. LDVSL >=1, and
98 * if JOBVSL = 'V', LDVSL >= N.
99 *
100 * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
101 * If JOBVSR = 'V', the matrix of right Schur vectors Z.
102 * Not referenced if JOBVSR = 'N'.
103 *
104 * LDVSR (input) INTEGER
105 * The leading dimension of the matrix VSR. LDVSR >= 1, and
106 * if JOBVSR = 'V', LDVSR >= N.
107 *
108 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
109 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
110 *
111 * LWORK (input) INTEGER
112 * The dimension of the array WORK. LWORK >= max(1,4*N).
113 * For good performance, LWORK must generally be larger.
114 * To compute the optimal value of LWORK, call ILAENV to get
115 * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
116 * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
117 * The optimal LWORK is 2*N + N*(NB+1).
118 *
119 * If LWORK = -1, then a workspace query is assumed; the routine
120 * only calculates the optimal size of the WORK array, returns
121 * this value as the first entry of the WORK array, and no error
122 * message related to LWORK is issued by XERBLA.
123 *
124 * INFO (output) INTEGER
125 * = 0: successful exit
126 * < 0: if INFO = -i, the i-th argument had an illegal value.
127 * = 1,...,N:
128 * The QZ iteration failed. (A,B) are not in Schur
129 * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
130 * be correct for j=INFO+1,...,N.
131 * > N: errors that usually indicate LAPACK problems:
132 * =N+1: error return from DGGBAL
133 * =N+2: error return from DGEQRF
134 * =N+3: error return from DORMQR
135 * =N+4: error return from DORGQR
136 * =N+5: error return from DGGHRD
137 * =N+6: error return from DHGEQZ (other than failed
138 * iteration)
139 * =N+7: error return from DGGBAK (computing VSL)
140 * =N+8: error return from DGGBAK (computing VSR)
141 * =N+9: error return from DLASCL (various places)
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146 DOUBLE PRECISION ZERO, ONE
147 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
148 * ..
149 * .. Local Scalars ..
150 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
151 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
152 $ IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
153 $ LWKOPT, NB, NB1, NB2, NB3
154 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
155 $ SAFMIN, SMLNUM
156 * ..
157 * .. External Subroutines ..
158 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
159 $ DLASCL, DLASET, DORGQR, DORMQR, XERBLA
160 * ..
161 * .. External Functions ..
162 LOGICAL LSAME
163 INTEGER ILAENV
164 DOUBLE PRECISION DLAMCH, DLANGE
165 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
166 * ..
167 * .. Intrinsic Functions ..
168 INTRINSIC INT, MAX
169 * ..
170 * .. Executable Statements ..
171 *
172 * Decode the input arguments
173 *
174 IF( LSAME( JOBVSL, 'N' ) ) THEN
175 IJOBVL = 1
176 ILVSL = .FALSE.
177 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
178 IJOBVL = 2
179 ILVSL = .TRUE.
180 ELSE
181 IJOBVL = -1
182 ILVSL = .FALSE.
183 END IF
184 *
185 IF( LSAME( JOBVSR, 'N' ) ) THEN
186 IJOBVR = 1
187 ILVSR = .FALSE.
188 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
189 IJOBVR = 2
190 ILVSR = .TRUE.
191 ELSE
192 IJOBVR = -1
193 ILVSR = .FALSE.
194 END IF
195 *
196 * Test the input arguments
197 *
198 LWKMIN = MAX( 4*N, 1 )
199 LWKOPT = LWKMIN
200 WORK( 1 ) = LWKOPT
201 LQUERY = ( LWORK.EQ.-1 )
202 INFO = 0
203 IF( IJOBVL.LE.0 ) THEN
204 INFO = -1
205 ELSE IF( IJOBVR.LE.0 ) THEN
206 INFO = -2
207 ELSE IF( N.LT.0 ) THEN
208 INFO = -3
209 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
210 INFO = -5
211 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
212 INFO = -7
213 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
214 INFO = -12
215 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
216 INFO = -14
217 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
218 INFO = -16
219 END IF
220 *
221 IF( INFO.EQ.0 ) THEN
222 NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
223 NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
224 NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
225 NB = MAX( NB1, NB2, NB3 )
226 LOPT = 2*N + N*( NB+1 )
227 WORK( 1 ) = LOPT
228 END IF
229 *
230 IF( INFO.NE.0 ) THEN
231 CALL XERBLA( 'DGEGS ', -INFO )
232 RETURN
233 ELSE IF( LQUERY ) THEN
234 RETURN
235 END IF
236 *
237 * Quick return if possible
238 *
239 IF( N.EQ.0 )
240 $ RETURN
241 *
242 * Get machine constants
243 *
244 EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
245 SAFMIN = DLAMCH( 'S' )
246 SMLNUM = N*SAFMIN / EPS
247 BIGNUM = ONE / SMLNUM
248 *
249 * Scale A if max element outside range [SMLNUM,BIGNUM]
250 *
251 ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
252 ILASCL = .FALSE.
253 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
254 ANRMTO = SMLNUM
255 ILASCL = .TRUE.
256 ELSE IF( ANRM.GT.BIGNUM ) THEN
257 ANRMTO = BIGNUM
258 ILASCL = .TRUE.
259 END IF
260 *
261 IF( ILASCL ) THEN
262 CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
263 IF( IINFO.NE.0 ) THEN
264 INFO = N + 9
265 RETURN
266 END IF
267 END IF
268 *
269 * Scale B if max element outside range [SMLNUM,BIGNUM]
270 *
271 BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
272 ILBSCL = .FALSE.
273 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
274 BNRMTO = SMLNUM
275 ILBSCL = .TRUE.
276 ELSE IF( BNRM.GT.BIGNUM ) THEN
277 BNRMTO = BIGNUM
278 ILBSCL = .TRUE.
279 END IF
280 *
281 IF( ILBSCL ) THEN
282 CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
283 IF( IINFO.NE.0 ) THEN
284 INFO = N + 9
285 RETURN
286 END IF
287 END IF
288 *
289 * Permute the matrix to make it more nearly triangular
290 * Workspace layout: (2*N words -- "work..." not actually used)
291 * left_permutation, right_permutation, work...
292 *
293 ILEFT = 1
294 IRIGHT = N + 1
295 IWORK = IRIGHT + N
296 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
297 $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
298 IF( IINFO.NE.0 ) THEN
299 INFO = N + 1
300 GO TO 10
301 END IF
302 *
303 * Reduce B to triangular form, and initialize VSL and/or VSR
304 * Workspace layout: ("work..." must have at least N words)
305 * left_permutation, right_permutation, tau, work...
306 *
307 IROWS = IHI + 1 - ILO
308 ICOLS = N + 1 - ILO
309 ITAU = IWORK
310 IWORK = ITAU + IROWS
311 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
312 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
313 IF( IINFO.GE.0 )
314 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
315 IF( IINFO.NE.0 ) THEN
316 INFO = N + 2
317 GO TO 10
318 END IF
319 *
320 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
321 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
322 $ LWORK+1-IWORK, IINFO )
323 IF( IINFO.GE.0 )
324 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
325 IF( IINFO.NE.0 ) THEN
326 INFO = N + 3
327 GO TO 10
328 END IF
329 *
330 IF( ILVSL ) THEN
331 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
332 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
333 $ VSL( ILO+1, ILO ), LDVSL )
334 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
335 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
336 $ IINFO )
337 IF( IINFO.GE.0 )
338 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
339 IF( IINFO.NE.0 ) THEN
340 INFO = N + 4
341 GO TO 10
342 END IF
343 END IF
344 *
345 IF( ILVSR )
346 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
347 *
348 * Reduce to generalized Hessenberg form
349 *
350 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
351 $ LDVSL, VSR, LDVSR, IINFO )
352 IF( IINFO.NE.0 ) THEN
353 INFO = N + 5
354 GO TO 10
355 END IF
356 *
357 * Perform QZ algorithm, computing Schur vectors if desired
358 * Workspace layout: ("work..." must have at least 1 word)
359 * left_permutation, right_permutation, work...
360 *
361 IWORK = ITAU
362 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
363 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
364 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
365 IF( IINFO.GE.0 )
366 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
367 IF( IINFO.NE.0 ) THEN
368 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
369 INFO = IINFO
370 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
371 INFO = IINFO - N
372 ELSE
373 INFO = N + 6
374 END IF
375 GO TO 10
376 END IF
377 *
378 * Apply permutation to VSL and VSR
379 *
380 IF( ILVSL ) THEN
381 CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
382 $ WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
383 IF( IINFO.NE.0 ) THEN
384 INFO = N + 7
385 GO TO 10
386 END IF
387 END IF
388 IF( ILVSR ) THEN
389 CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
390 $ WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
391 IF( IINFO.NE.0 ) THEN
392 INFO = N + 8
393 GO TO 10
394 END IF
395 END IF
396 *
397 * Undo scaling
398 *
399 IF( ILASCL ) THEN
400 CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
401 IF( IINFO.NE.0 ) THEN
402 INFO = N + 9
403 RETURN
404 END IF
405 CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
406 $ IINFO )
407 IF( IINFO.NE.0 ) THEN
408 INFO = N + 9
409 RETURN
410 END IF
411 CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
412 $ IINFO )
413 IF( IINFO.NE.0 ) THEN
414 INFO = N + 9
415 RETURN
416 END IF
417 END IF
418 *
419 IF( ILBSCL ) THEN
420 CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
421 IF( IINFO.NE.0 ) THEN
422 INFO = N + 9
423 RETURN
424 END IF
425 CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
426 IF( IINFO.NE.0 ) THEN
427 INFO = N + 9
428 RETURN
429 END IF
430 END IF
431 *
432 10 CONTINUE
433 WORK( 1 ) = LWKOPT
434 *
435 RETURN
436 *
437 * End of DGEGS
438 *
439 END