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*> \date November 2011
199*
200*> \ingroup doubleGEcomputational
201*
202*> \par Contributors:
203*  ==================
204*>
205*>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
206*>     Umea University, S-901 87 Umea, Sweden.
207*
208*> \par References:
209*  ================
210*>
211*> \verbatim
212*>
213*>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
214*>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
215*>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
216*>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
217*> \endverbatim
218*>
219*  =====================================================================
220      SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
221     $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
222*
223*  -- LAPACK computational routine (version 3.4.0) --
224*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
225*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226*     November 2011
227*
228*     .. Scalar Arguments ..
229      LOGICAL            WANTQ, WANTZ
230      INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
231*     ..
232*     .. Array Arguments ..
233      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
234     $                   WORK( * ), Z( LDZ, * )
235*     ..
236*
237*  =====================================================================
238*
239*     .. Parameters ..
240      DOUBLE PRECISION   ZERO
241      PARAMETER          ( ZERO = 0.0D+0 )
242*     ..
243*     .. Local Scalars ..
244      LOGICAL            LQUERY
245      INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
246*     ..
247*     .. External Subroutines ..
248      EXTERNAL           DTGEX2, XERBLA
249*     ..
250*     .. Intrinsic Functions ..
251      INTRINSIC          MAX
252*     ..
253*     .. Executable Statements ..
254*
255*     Decode and test input arguments.
256*
257      INFO = 0
258      LQUERY = ( LWORK.EQ.-1 )
259      IF( N.LT.0 ) THEN
260         INFO = -3
261      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
262         INFO = -5
263      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
264         INFO = -7
265      ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
266         INFO = -9
267      ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
268         INFO = -11
269      ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
270         INFO = -12
271      ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
272         INFO = -13
273      END IF
274*
275      IF( INFO.EQ.0 ) THEN
276         IF( N.LE.1 ) THEN
277            LWMIN = 1
278         ELSE
279            LWMIN = 4*N + 16
280         END IF
281         WORK(1) = LWMIN
282*
283         IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
284            INFO = -15
285         END IF
286      END IF
287*
288      IF( INFO.NE.0 ) THEN
289         CALL XERBLA( 'DTGEXC', -INFO )
290         RETURN
291      ELSE IF( LQUERY ) THEN
292         RETURN
293      END IF
294*
295*     Quick return if possible
296*
297      IF( N.LE.1 )
298     $   RETURN
299*
300*     Determine the first row of the specified block and find out
301*     if it is 1-by-1 or 2-by-2.
302*
303      IF( IFST.GT.1 ) THEN
304         IF( A( IFST, IFST-1 ).NE.ZERO )
305     $      IFST = IFST - 1
306      END IF
307      NBF = 1
308      IF( IFST.LT.N ) THEN
309         IF( A( IFST+1, IFST ).NE.ZERO )
310     $      NBF = 2
311      END IF
312*
313*     Determine the first row of the final block
314*     and find out if it is 1-by-1 or 2-by-2.
315*
316      IF( ILST.GT.1 ) THEN
317         IF( A( ILST, ILST-1 ).NE.ZERO )
318     $      ILST = ILST - 1
319      END IF
320      NBL = 1
321      IF( ILST.LT.N ) THEN
322         IF( A( ILST+1, ILST ).NE.ZERO )
323     $      NBL = 2
324      END IF
325      IF( IFST.EQ.ILST )
326     $   RETURN
327*
328      IF( IFST.LT.ILST ) THEN
329*
330*        Update ILST.
331*
332         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
333     $      ILST = ILST - 1
334         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
335     $      ILST = ILST + 1
336*
337         HERE = IFST
338*
339   10    CONTINUE
340*
341*        Swap with next one below.
342*
343         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
344*
345*           Current block either 1-by-1 or 2-by-2.
346*
347            NBNEXT = 1
348            IF( HERE+NBF+1.LE.N ) THEN
349               IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
350     $            NBNEXT = 2
351            END IF
352            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
353     $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
354            IF( INFO.NE.0 ) THEN
355               ILST = HERE
356               RETURN
357            END IF
358            HERE = HERE + NBNEXT
359*
360*           Test if 2-by-2 block breaks into two 1-by-1 blocks.
361*
362            IF( NBF.EQ.2 ) THEN
363               IF( A( HERE+1, HERE ).EQ.ZERO )
364     $            NBF = 3
365            END IF
366*
367         ELSE
368*
369*           Current block consists of two 1-by-1 blocks, each of which
370*           must be swapped individually.
371*
372            NBNEXT = 1
373            IF( HERE+3.LE.N ) THEN
374               IF( A( HERE+3, HERE+2 ).NE.ZERO )
375     $            NBNEXT = 2
376            END IF
377            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
378     $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
379            IF( INFO.NE.0 ) THEN
380               ILST = HERE
381               RETURN
382            END IF
383            IF( NBNEXT.EQ.1 ) THEN
384*
385*              Swap two 1-by-1 blocks.
386*
387               CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
388     $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
389               IF( INFO.NE.0 ) THEN
390                  ILST = HERE
391                  RETURN
392               END IF
393               HERE = HERE + 1
394*
395            ELSE
396*
397*              Recompute NBNEXT in case of 2-by-2 split.
398*
399               IF( A( HERE+2, HERE+1 ).EQ.ZERO )
400     $            NBNEXT = 1
401               IF( NBNEXT.EQ.2 ) THEN
402*
403*                 2-by-2 block did not split.
404*
405                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
406     $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
407     $                         INFO )
408                  IF( INFO.NE.0 ) THEN
409                     ILST = HERE
410                     RETURN
411                  END IF
412                  HERE = HERE + 2
413               ELSE
414*
415*                 2-by-2 block did split.
416*
417                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
418     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
419                  IF( INFO.NE.0 ) THEN
420                     ILST = HERE
421                     RETURN
422                  END IF
423                  HERE = HERE + 1
424                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
425     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
426                  IF( INFO.NE.0 ) THEN
427                     ILST = HERE
428                     RETURN
429                  END IF
430                  HERE = HERE + 1
431               END IF
432*
433            END IF
434         END IF
435         IF( HERE.LT.ILST )
436     $      GO TO 10
437      ELSE
438         HERE = IFST
439*
440   20    CONTINUE
441*
442*        Swap with next one below.
443*
444         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
445*
446*           Current block either 1-by-1 or 2-by-2.
447*
448            NBNEXT = 1
449            IF( HERE.GE.3 ) THEN
450               IF( A( HERE-1, HERE-2 ).NE.ZERO )
451     $            NBNEXT = 2
452            END IF
453            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
454     $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
455     $                   INFO )
456            IF( INFO.NE.0 ) THEN
457               ILST = HERE
458               RETURN
459            END IF
460            HERE = HERE - NBNEXT
461*
462*           Test if 2-by-2 block breaks into two 1-by-1 blocks.
463*
464            IF( NBF.EQ.2 ) THEN
465               IF( A( HERE+1, HERE ).EQ.ZERO )
466     $            NBF = 3
467            END IF
468*
469         ELSE
470*
471*           Current block consists of two 1-by-1 blocks, each of which
472*           must be swapped individually.
473*
474            NBNEXT = 1
475            IF( HERE.GE.3 ) THEN
476               IF( A( HERE-1, HERE-2 ).NE.ZERO )
477     $            NBNEXT = 2
478            END IF
479            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
480     $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
481     $                   INFO )
482            IF( INFO.NE.0 ) THEN
483               ILST = HERE
484               RETURN
485            END IF
486            IF( NBNEXT.EQ.1 ) THEN
487*
488*              Swap two 1-by-1 blocks.
489*
490               CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
491     $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
492               IF( INFO.NE.0 ) THEN
493                  ILST = HERE
494                  RETURN
495               END IF
496               HERE = HERE - 1
497            ELSE
498*
499*             Recompute NBNEXT in case of 2-by-2 split.
500*
501               IF( A( HERE, HERE-1 ).EQ.ZERO )
502     $            NBNEXT = 1
503               IF( NBNEXT.EQ.2 ) THEN
504*
505*                 2-by-2 block did not split.
506*
507                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
508     $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
509                  IF( INFO.NE.0 ) THEN
510                     ILST = HERE
511                     RETURN
512                  END IF
513                  HERE = HERE - 2
514               ELSE
515*
516*                 2-by-2 block did split.
517*
518                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
519     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
520                  IF( INFO.NE.0 ) THEN
521                     ILST = HERE
522                     RETURN
523                  END IF
524                  HERE = HERE - 1
525                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
526     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
527                  IF( INFO.NE.0 ) THEN
528                     ILST = HERE
529                     RETURN
530                  END IF
531                  HERE = HERE - 1
532               END IF
533            END IF
534         END IF
535         IF( HERE.GT.ILST )
536     $      GO TO 20
537      END IF
538      ILST = HERE
539      WORK( 1 ) = LWMIN
540      RETURN
541*
542*     End of DTGEXC
543*
544      END
545