LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
stgexc.f
Go to the documentation of this file.
1 *> \brief \b STGEXC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STGEXC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgexc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgexc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgexc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STGEXC( 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 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30 * $ WORK( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> STGEXC 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 SGGES), 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 REAL 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 REAL 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 REAL array, dimension (LDZ,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 REAL 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 REAL 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 realGEcomputational
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 stgexc( 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  REAL a( lda, * ), b( ldb, * ), q( ldq, * ),
234  $ work( * ), z( ldz, * )
235 * ..
236 *
237 * =====================================================================
238 *
239 * .. Parameters ..
240  REAL zero
241  parameter( zero = 0.0e+0 )
242 * ..
243 * .. Local Scalars ..
244  LOGICAL lquery
245  INTEGER here, lwmin, nbf, nbl, nbnext
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL stgex2, 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( 'STGEXC', -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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 STGEXC
543 *
544  END