LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtgexc.f
Go to the documentation of this file.
1*> \brief \b DTGEXC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DTGEXC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgexc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgexc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgexc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
22* LDZ, IFST, ILST, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* LOGICAL WANTQ, WANTZ
26* INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30* $ WORK( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DTGEXC reorders the generalized real Schur decomposition of a real
40*> matrix pair (A,B) using an orthogonal equivalence transformation
41*>
42*> (A, B) = Q * (A, B) * Z**T,
43*>
44*> so that the diagonal block of (A, B) with row index IFST is moved
45*> to row ILST.
46*>
47*> (A, B) must be in generalized real Schur canonical form (as returned
48*> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
49*> diagonal blocks. B is upper triangular.
50*>
51*> Optionally, the matrices Q and Z of generalized Schur vectors are
52*> updated.
53*>
54*> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
55*> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
56*>
57*> \endverbatim
58*
59* Arguments:
60* ==========
61*
62*> \param[in] WANTQ
63*> \verbatim
64*> WANTQ is LOGICAL
65*> .TRUE. : update the left transformation matrix Q;
66*> .FALSE.: do not update Q.
67*> \endverbatim
68*>
69*> \param[in] WANTZ
70*> \verbatim
71*> WANTZ is LOGICAL
72*> .TRUE. : update the right transformation matrix Z;
73*> .FALSE.: do not update Z.
74*> \endverbatim
75*>
76*> \param[in] N
77*> \verbatim
78*> N is INTEGER
79*> The order of the matrices A and B. N >= 0.
80*> \endverbatim
81*>
82*> \param[in,out] A
83*> \verbatim
84*> A is DOUBLE PRECISION array, dimension (LDA,N)
85*> On entry, the matrix A in generalized real Schur canonical
86*> form.
87*> On exit, the updated matrix A, again in generalized
88*> real Schur canonical form.
89*> \endverbatim
90*>
91*> \param[in] LDA
92*> \verbatim
93*> LDA is INTEGER
94*> The leading dimension of the array A. LDA >= max(1,N).
95*> \endverbatim
96*>
97*> \param[in,out] B
98*> \verbatim
99*> B is DOUBLE PRECISION array, dimension (LDB,N)
100*> On entry, the matrix B in generalized real Schur canonical
101*> form (A,B).
102*> On exit, the updated matrix B, again in generalized
103*> real Schur canonical form (A,B).
104*> \endverbatim
105*>
106*> \param[in] LDB
107*> \verbatim
108*> LDB is INTEGER
109*> The leading dimension of the array B. LDB >= max(1,N).
110*> \endverbatim
111*>
112*> \param[in,out] Q
113*> \verbatim
114*> Q is DOUBLE PRECISION array, dimension (LDQ,N)
115*> On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
116*> On exit, the updated matrix Q.
117*> If WANTQ = .FALSE., Q is not referenced.
118*> \endverbatim
119*>
120*> \param[in] LDQ
121*> \verbatim
122*> LDQ is INTEGER
123*> The leading dimension of the array Q. LDQ >= 1.
124*> If WANTQ = .TRUE., LDQ >= N.
125*> \endverbatim
126*>
127*> \param[in,out] Z
128*> \verbatim
129*> Z is DOUBLE PRECISION array, dimension (LDZ,N)
130*> On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
131*> On exit, the updated matrix Z.
132*> If WANTZ = .FALSE., Z is not referenced.
133*> \endverbatim
134*>
135*> \param[in] LDZ
136*> \verbatim
137*> LDZ is INTEGER
138*> The leading dimension of the array Z. LDZ >= 1.
139*> If WANTZ = .TRUE., LDZ >= N.
140*> \endverbatim
141*>
142*> \param[in,out] IFST
143*> \verbatim
144*> IFST is INTEGER
145*> \endverbatim
146*>
147*> \param[in,out] ILST
148*> \verbatim
149*> ILST is INTEGER
150*> Specify the reordering of the diagonal blocks of (A, B).
151*> The block with row index IFST is moved to row ILST, by a
152*> sequence of swapping between adjacent blocks.
153*> On exit, if IFST pointed on entry to the second row of
154*> a 2-by-2 block, it is changed to point to the first row;
155*> ILST always points to the first row of the block in its
156*> final position (which may differ from its input value by
157*> +1 or -1). 1 <= IFST, ILST <= N.
158*> \endverbatim
159*>
160*> \param[out] WORK
161*> \verbatim
162*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
163*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164*> \endverbatim
165*>
166*> \param[in] LWORK
167*> \verbatim
168*> LWORK is INTEGER
169*> The dimension of the array WORK.
170*> LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
171*>
172*> If LWORK = -1, then a workspace query is assumed; the routine
173*> only calculates the optimal size of the WORK array, returns
174*> this value as the first entry of the WORK array, and no error
175*> message related to LWORK is issued by XERBLA.
176*> \endverbatim
177*>
178*> \param[out] INFO
179*> \verbatim
180*> INFO is INTEGER
181*> =0: successful exit.
182*> <0: if INFO = -i, the i-th argument had an illegal value.
183*> =1: The transformed matrix pair (A, B) would be too far
184*> from generalized Schur form; the problem is ill-
185*> conditioned. (A, B) may have been partially reordered,
186*> and ILST points to the first row of the current
187*> position of the block being moved.
188*> \endverbatim
189*
190* Authors:
191* ========
192*
193*> \author Univ. of Tennessee
194*> \author Univ. of California Berkeley
195*> \author Univ. of Colorado Denver
196*> \author NAG Ltd.
197*
198*> \ingroup tgexc
199*
200*> \par Contributors:
201* ==================
202*>
203*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
204*> Umea University, S-901 87 Umea, Sweden.
205*
206*> \par References:
207* ================
208*>
209*> \verbatim
210*>
211*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
212*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
213*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
214*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
215*> \endverbatim
216*>
217* =====================================================================
218 SUBROUTINE dtgexc( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
219 $ LDZ, IFST, ILST, WORK, LWORK, INFO )
220*
221* -- LAPACK computational routine --
222* -- LAPACK is a software package provided by Univ. of Tennessee, --
223* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224*
225* .. Scalar Arguments ..
226 LOGICAL WANTQ, WANTZ
227 INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
228* ..
229* .. Array Arguments ..
230 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
231 $ work( * ), z( ldz, * )
232* ..
233*
234* =====================================================================
235*
236* .. Parameters ..
237 DOUBLE PRECISION ZERO
238 parameter( zero = 0.0d+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY
242 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
243* ..
244* .. External Subroutines ..
245 EXTERNAL dtgex2, xerbla
246* ..
247* .. Intrinsic Functions ..
248 INTRINSIC max
249* ..
250* .. Executable Statements ..
251*
252* Decode and test input arguments.
253*
254 info = 0
255 lquery = ( lwork.EQ.-1 )
256 IF( n.LT.0 ) THEN
257 info = -3
258 ELSE IF( lda.LT.max( 1, n ) ) THEN
259 info = -5
260 ELSE IF( ldb.LT.max( 1, n ) ) THEN
261 info = -7
262 ELSE IF( ldq.LT.1 .OR. wantq .AND. ( ldq.LT.max( 1, n ) ) ) THEN
263 info = -9
264 ELSE IF( ldz.LT.1 .OR. wantz .AND. ( ldz.LT.max( 1, n ) ) ) THEN
265 info = -11
266 ELSE IF( ifst.LT.1 .OR. ifst.GT.n ) THEN
267 info = -12
268 ELSE IF( ilst.LT.1 .OR. ilst.GT.n ) THEN
269 info = -13
270 END IF
271*
272 IF( info.EQ.0 ) THEN
273 IF( n.LE.1 ) THEN
274 lwmin = 1
275 ELSE
276 lwmin = 4*n + 16
277 END IF
278 work(1) = lwmin
279*
280 IF (lwork.LT.lwmin .AND. .NOT.lquery) THEN
281 info = -15
282 END IF
283 END IF
284*
285 IF( info.NE.0 ) THEN
286 CALL xerbla( 'DTGEXC', -info )
287 RETURN
288 ELSE IF( lquery ) THEN
289 RETURN
290 END IF
291*
292* Quick return if possible
293*
294 IF( n.LE.1 )
295 $ RETURN
296*
297* Determine the first row of the specified block and find out
298* if it is 1-by-1 or 2-by-2.
299*
300 IF( ifst.GT.1 ) THEN
301 IF( a( ifst, ifst-1 ).NE.zero )
302 $ ifst = ifst - 1
303 END IF
304 nbf = 1
305 IF( ifst.LT.n ) THEN
306 IF( a( ifst+1, ifst ).NE.zero )
307 $ nbf = 2
308 END IF
309*
310* Determine the first row of the final block
311* and find out if it is 1-by-1 or 2-by-2.
312*
313 IF( ilst.GT.1 ) THEN
314 IF( a( ilst, ilst-1 ).NE.zero )
315 $ ilst = ilst - 1
316 END IF
317 nbl = 1
318 IF( ilst.LT.n ) THEN
319 IF( a( ilst+1, ilst ).NE.zero )
320 $ nbl = 2
321 END IF
322 IF( ifst.EQ.ilst )
323 $ RETURN
324*
325 IF( ifst.LT.ilst ) THEN
326*
327* Update ILST.
328*
329 IF( nbf.EQ.2 .AND. nbl.EQ.1 )
330 $ ilst = ilst - 1
331 IF( nbf.EQ.1 .AND. nbl.EQ.2 )
332 $ ilst = ilst + 1
333*
334 here = ifst
335*
336 10 CONTINUE
337*
338* Swap with next one below.
339*
340 IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
341*
342* Current block either 1-by-1 or 2-by-2.
343*
344 nbnext = 1
345 IF( here+nbf+1.LE.n ) THEN
346 IF( a( here+nbf+1, here+nbf ).NE.zero )
347 $ nbnext = 2
348 END IF
349 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
350 $ ldz, here, nbf, nbnext, work, lwork, info )
351 IF( info.NE.0 ) THEN
352 ilst = here
353 RETURN
354 END IF
355 here = here + nbnext
356*
357* Test if 2-by-2 block breaks into two 1-by-1 blocks.
358*
359 IF( nbf.EQ.2 ) THEN
360 IF( a( here+1, here ).EQ.zero )
361 $ nbf = 3
362 END IF
363*
364 ELSE
365*
366* Current block consists of two 1-by-1 blocks, each of which
367* must be swapped individually.
368*
369 nbnext = 1
370 IF( here+3.LE.n ) THEN
371 IF( a( here+3, here+2 ).NE.zero )
372 $ nbnext = 2
373 END IF
374 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
375 $ ldz, here+1, 1, nbnext, work, lwork, info )
376 IF( info.NE.0 ) THEN
377 ilst = here
378 RETURN
379 END IF
380 IF( nbnext.EQ.1 ) THEN
381*
382* Swap two 1-by-1 blocks.
383*
384 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
385 $ ldz, here, 1, 1, work, lwork, info )
386 IF( info.NE.0 ) THEN
387 ilst = here
388 RETURN
389 END IF
390 here = here + 1
391*
392 ELSE
393*
394* Recompute NBNEXT in case of 2-by-2 split.
395*
396 IF( a( here+2, here+1 ).EQ.zero )
397 $ nbnext = 1
398 IF( nbnext.EQ.2 ) THEN
399*
400* 2-by-2 block did not split.
401*
402 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
403 $ z, ldz, here, 1, nbnext, work, lwork,
404 $ info )
405 IF( info.NE.0 ) THEN
406 ilst = here
407 RETURN
408 END IF
409 here = here + 2
410 ELSE
411*
412* 2-by-2 block did split.
413*
414 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
415 $ z, ldz, here, 1, 1, work, lwork, info )
416 IF( info.NE.0 ) THEN
417 ilst = here
418 RETURN
419 END IF
420 here = here + 1
421 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
422 $ z, ldz, here, 1, 1, work, lwork, info )
423 IF( info.NE.0 ) THEN
424 ilst = here
425 RETURN
426 END IF
427 here = here + 1
428 END IF
429*
430 END IF
431 END IF
432 IF( here.LT.ilst )
433 $ GO TO 10
434 ELSE
435 here = ifst
436*
437 20 CONTINUE
438*
439* Swap with next one below.
440*
441 IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
442*
443* Current block either 1-by-1 or 2-by-2.
444*
445 nbnext = 1
446 IF( here.GE.3 ) THEN
447 IF( a( here-1, here-2 ).NE.zero )
448 $ nbnext = 2
449 END IF
450 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
451 $ ldz, here-nbnext, nbnext, nbf, work, lwork,
452 $ info )
453 IF( info.NE.0 ) THEN
454 ilst = here
455 RETURN
456 END IF
457 here = here - nbnext
458*
459* Test if 2-by-2 block breaks into two 1-by-1 blocks.
460*
461 IF( nbf.EQ.2 ) THEN
462 IF( a( here+1, here ).EQ.zero )
463 $ nbf = 3
464 END IF
465*
466 ELSE
467*
468* Current block consists of two 1-by-1 blocks, each of which
469* must be swapped individually.
470*
471 nbnext = 1
472 IF( here.GE.3 ) THEN
473 IF( a( here-1, here-2 ).NE.zero )
474 $ nbnext = 2
475 END IF
476 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
477 $ ldz, here-nbnext, nbnext, 1, work, lwork,
478 $ info )
479 IF( info.NE.0 ) THEN
480 ilst = here
481 RETURN
482 END IF
483 IF( nbnext.EQ.1 ) THEN
484*
485* Swap two 1-by-1 blocks.
486*
487 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
488 $ ldz, here, nbnext, 1, work, lwork, info )
489 IF( info.NE.0 ) THEN
490 ilst = here
491 RETURN
492 END IF
493 here = here - 1
494 ELSE
495*
496* Recompute NBNEXT in case of 2-by-2 split.
497*
498 IF( a( here, here-1 ).EQ.zero )
499 $ nbnext = 1
500 IF( nbnext.EQ.2 ) THEN
501*
502* 2-by-2 block did not split.
503*
504 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
505 $ z, ldz, here-1, 2, 1, work, lwork, info )
506 IF( info.NE.0 ) THEN
507 ilst = here
508 RETURN
509 END IF
510 here = here - 2
511 ELSE
512*
513* 2-by-2 block did split.
514*
515 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
516 $ z, ldz, here, 1, 1, work, lwork, info )
517 IF( info.NE.0 ) THEN
518 ilst = here
519 RETURN
520 END IF
521 here = here - 1
522 CALL dtgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
523 $ z, ldz, here, 1, 1, work, lwork, info )
524 IF( info.NE.0 ) THEN
525 ilst = here
526 RETURN
527 END IF
528 here = here - 1
529 END IF
530 END IF
531 END IF
532 IF( here.GT.ilst )
533 $ GO TO 20
534 END IF
535 ilst = here
536 work( 1 ) = lwmin
537 RETURN
538*
539* End of DTGEXC
540*
541 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dtgex2
subroutine dtgex2(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, n1, n2, work, lwork, info)
DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition dtgex2.f:221
dtgexc
subroutine dtgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
DTGEXC
Definition dtgexc.f:220