LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 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 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 Functions ..
245 REAL SROUNDUP_LWORK
246 EXTERNAL sroundup_lwork
247* ..
248* .. External Subroutines ..
249 EXTERNAL stgex2, xerbla
250* ..
251* .. Intrinsic Functions ..
252 INTRINSIC max
253* ..
254* .. Executable Statements ..
255*
256* Decode and test input arguments.
257*
258 info = 0
259 lquery = ( lwork.EQ.-1 )
260 IF( n.LT.0 ) THEN
261 info = -3
262 ELSE IF( lda.LT.max( 1, n ) ) THEN
263 info = -5
264 ELSE IF( ldb.LT.max( 1, n ) ) THEN
265 info = -7
266 ELSE IF( ldq.LT.1 .OR. wantq .AND. ( ldq.LT.max( 1, n ) ) ) THEN
267 info = -9
268 ELSE IF( ldz.LT.1 .OR. wantz .AND. ( ldz.LT.max( 1, n ) ) ) THEN
269 info = -11
270 ELSE IF( ifst.LT.1 .OR. ifst.GT.n ) THEN
271 info = -12
272 ELSE IF( ilst.LT.1 .OR. ilst.GT.n ) THEN
273 info = -13
274 END IF
275*
276 IF( info.EQ.0 ) THEN
277 IF( n.LE.1 ) THEN
278 lwmin = 1
279 ELSE
280 lwmin = 4*n + 16
281 END IF
282 work(1) = lwmin
283*
284 IF (lwork.LT.lwmin .AND. .NOT.lquery) THEN
285 info = -15
286 END IF
287 END IF
288*
289 IF( info.NE.0 ) THEN
290 CALL xerbla( 'STGEXC', -info )
291 RETURN
292 ELSE IF( lquery ) THEN
293 RETURN
294 END IF
295*
296* Quick return if possible
297*
298 IF( n.LE.1 )
299 $ RETURN
300*
301* Determine the first row of the specified block and find out
302* if it is 1-by-1 or 2-by-2.
303*
304 IF( ifst.GT.1 ) THEN
305 IF( a( ifst, ifst-1 ).NE.zero )
306 $ ifst = ifst - 1
307 END IF
308 nbf = 1
309 IF( ifst.LT.n ) THEN
310 IF( a( ifst+1, ifst ).NE.zero )
311 $ nbf = 2
312 END IF
313*
314* Determine the first row of the final block
315* and find out if it is 1-by-1 or 2-by-2.
316*
317 IF( ilst.GT.1 ) THEN
318 IF( a( ilst, ilst-1 ).NE.zero )
319 $ ilst = ilst - 1
320 END IF
321 nbl = 1
322 IF( ilst.LT.n ) THEN
323 IF( a( ilst+1, ilst ).NE.zero )
324 $ nbl = 2
325 END IF
326 IF( ifst.EQ.ilst )
327 $ RETURN
328*
329 IF( ifst.LT.ilst ) THEN
330*
331* Update ILST.
332*
333 IF( nbf.EQ.2 .AND. nbl.EQ.1 )
334 $ ilst = ilst - 1
335 IF( nbf.EQ.1 .AND. nbl.EQ.2 )
336 $ ilst = ilst + 1
337*
338 here = ifst
339*
340 10 CONTINUE
341*
342* Swap with next one below.
343*
344 IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
345*
346* Current block either 1-by-1 or 2-by-2.
347*
348 nbnext = 1
349 IF( here+nbf+1.LE.n ) THEN
350 IF( a( here+nbf+1, here+nbf ).NE.zero )
351 $ nbnext = 2
352 END IF
353 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
354 $ ldz, here, nbf, nbnext, work, lwork, info )
355 IF( info.NE.0 ) THEN
356 ilst = here
357 RETURN
358 END IF
359 here = here + nbnext
360*
361* Test if 2-by-2 block breaks into two 1-by-1 blocks.
362*
363 IF( nbf.EQ.2 ) THEN
364 IF( a( here+1, here ).EQ.zero )
365 $ nbf = 3
366 END IF
367*
368 ELSE
369*
370* Current block consists of two 1-by-1 blocks, each of which
371* must be swapped individually.
372*
373 nbnext = 1
374 IF( here+3.LE.n ) THEN
375 IF( a( here+3, here+2 ).NE.zero )
376 $ nbnext = 2
377 END IF
378 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
379 $ ldz, here+1, 1, nbnext, work, lwork, info )
380 IF( info.NE.0 ) THEN
381 ilst = here
382 RETURN
383 END IF
384 IF( nbnext.EQ.1 ) THEN
385*
386* Swap two 1-by-1 blocks.
387*
388 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
389 $ ldz, here, 1, 1, work, lwork, info )
390 IF( info.NE.0 ) THEN
391 ilst = here
392 RETURN
393 END IF
394 here = here + 1
395*
396 ELSE
397*
398* Recompute NBNEXT in case of 2-by-2 split.
399*
400 IF( a( here+2, here+1 ).EQ.zero )
401 $ nbnext = 1
402 IF( nbnext.EQ.2 ) THEN
403*
404* 2-by-2 block did not split.
405*
406 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
407 $ z, ldz, here, 1, nbnext, work, lwork,
408 $ info )
409 IF( info.NE.0 ) THEN
410 ilst = here
411 RETURN
412 END IF
413 here = here + 2
414 ELSE
415*
416* 2-by-2 block did split.
417*
418 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
419 $ z, ldz, here, 1, 1, work, lwork, info )
420 IF( info.NE.0 ) THEN
421 ilst = here
422 RETURN
423 END IF
424 here = here + 1
425 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
426 $ z, ldz, here, 1, 1, work, lwork, info )
427 IF( info.NE.0 ) THEN
428 ilst = here
429 RETURN
430 END IF
431 here = here + 1
432 END IF
433*
434 END IF
435 END IF
436 IF( here.LT.ilst )
437 $ GO TO 10
438 ELSE
439 here = ifst
440*
441 20 CONTINUE
442*
443* Swap with next one below.
444*
445 IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
446*
447* Current block either 1-by-1 or 2-by-2.
448*
449 nbnext = 1
450 IF( here.GE.3 ) THEN
451 IF( a( here-1, here-2 ).NE.zero )
452 $ nbnext = 2
453 END IF
454 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
455 $ ldz, here-nbnext, nbnext, nbf, work, lwork,
456 $ info )
457 IF( info.NE.0 ) THEN
458 ilst = here
459 RETURN
460 END IF
461 here = here - nbnext
462*
463* Test if 2-by-2 block breaks into two 1-by-1 blocks.
464*
465 IF( nbf.EQ.2 ) THEN
466 IF( a( here+1, here ).EQ.zero )
467 $ nbf = 3
468 END IF
469*
470 ELSE
471*
472* Current block consists of two 1-by-1 blocks, each of which
473* must be swapped individually.
474*
475 nbnext = 1
476 IF( here.GE.3 ) THEN
477 IF( a( here-1, here-2 ).NE.zero )
478 $ nbnext = 2
479 END IF
480 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
481 $ ldz, here-nbnext, nbnext, 1, work, lwork,
482 $ info )
483 IF( info.NE.0 ) THEN
484 ilst = here
485 RETURN
486 END IF
487 IF( nbnext.EQ.1 ) THEN
488*
489* Swap two 1-by-1 blocks.
490*
491 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
492 $ ldz, here, nbnext, 1, work, lwork, info )
493 IF( info.NE.0 ) THEN
494 ilst = here
495 RETURN
496 END IF
497 here = here - 1
498 ELSE
499*
500* Recompute NBNEXT in case of 2-by-2 split.
501*
502 IF( a( here, here-1 ).EQ.zero )
503 $ nbnext = 1
504 IF( nbnext.EQ.2 ) THEN
505*
506* 2-by-2 block did not split.
507*
508 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
509 $ z, ldz, here-1, 2, 1, work, lwork, info )
510 IF( info.NE.0 ) THEN
511 ilst = here
512 RETURN
513 END IF
514 here = here - 2
515 ELSE
516*
517* 2-by-2 block did split.
518*
519 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
520 $ z, ldz, here, 1, 1, work, lwork, info )
521 IF( info.NE.0 ) THEN
522 ilst = here
523 RETURN
524 END IF
525 here = here - 1
526 CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
527 $ z, ldz, here, 1, 1, work, lwork, info )
528 IF( info.NE.0 ) THEN
529 ilst = here
530 RETURN
531 END IF
532 here = here - 1
533 END IF
534 END IF
535 END IF
536 IF( here.GT.ilst )
537 $ GO TO 20
538 END IF
539 ilst = here
540 work( 1 ) = sroundup_lwork(lwmin)
541 RETURN
542*
543* End of STGEXC
544*
545 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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