LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 (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 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 *> \ingroup realGEcomputational
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 stgexc( 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  REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
231  $ work( * ), z( ldz, * )
232 * ..
233 *
234 * =====================================================================
235 *
236 * .. Parameters ..
237  REAL ZERO
238  parameter( zero = 0.0e+0 )
239 * ..
240 * .. Local Scalars ..
241  LOGICAL LQUERY
242  INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
243 * ..
244 * .. External Subroutines ..
245  EXTERNAL stgex2, 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( 'STGEXC', -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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 stgex2( 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 STGEXC
540 *
541  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine stgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
STGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition: stgex2.f:221
subroutine stgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
STGEXC
Definition: stgexc.f:220