LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
strexc.f
Go to the documentation of this file.
1 *> \brief \b STREXC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STREXC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strexc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strexc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strexc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ
26 * INTEGER IFST, ILST, INFO, LDQ, LDT, N
27 * ..
28 * .. Array Arguments ..
29 * REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> STREXC reorders the real Schur factorization of a real matrix
39 *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
40 *> moved to row ILST.
41 *>
42 *> The real Schur form T is reordered by an orthogonal similarity
43 *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
44 *> is updated by postmultiplying it with Z.
45 *>
46 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
47 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
48 *> 2-by-2 diagonal block has its diagonal elements equal and its
49 *> off-diagonal elements of opposite sign.
50 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] COMPQ
56 *> \verbatim
57 *> COMPQ is CHARACTER*1
58 *> = 'V': update the matrix Q of Schur vectors;
59 *> = 'N': do not update Q.
60 *> \endverbatim
61 *>
62 *> \param[in] N
63 *> \verbatim
64 *> N is INTEGER
65 *> The order of the matrix T. N >= 0.
66 *> \endverbatim
67 *>
68 *> \param[in,out] T
69 *> \verbatim
70 *> T is REAL array, dimension (LDT,N)
71 *> On entry, the upper quasi-triangular matrix T, in Schur
72 *> Schur canonical form.
73 *> On exit, the reordered upper quasi-triangular matrix, again
74 *> in Schur canonical form.
75 *> \endverbatim
76 *>
77 *> \param[in] LDT
78 *> \verbatim
79 *> LDT is INTEGER
80 *> The leading dimension of the array T. LDT >= max(1,N).
81 *> \endverbatim
82 *>
83 *> \param[in,out] Q
84 *> \verbatim
85 *> Q is REAL array, dimension (LDQ,N)
86 *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
87 *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
88 *> orthogonal transformation matrix Z which reorders T.
89 *> If COMPQ = 'N', Q is not referenced.
90 *> \endverbatim
91 *>
92 *> \param[in] LDQ
93 *> \verbatim
94 *> LDQ is INTEGER
95 *> The leading dimension of the array Q. LDQ >= max(1,N).
96 *> \endverbatim
97 *>
98 *> \param[in,out] IFST
99 *> \verbatim
100 *> IFST is INTEGER
101 *> \endverbatim
102 *>
103 *> \param[in,out] ILST
104 *> \verbatim
105 *> ILST is INTEGER
106 *>
107 *> Specify the reordering of the diagonal blocks of T.
108 *> The block with row index IFST is moved to row ILST, by a
109 *> sequence of transpositions between adjacent blocks.
110 *> On exit, if IFST pointed on entry to the second row of a
111 *> 2-by-2 block, it is changed to point to the first row; ILST
112 *> always points to the first row of the block in its final
113 *> position (which may differ from its input value by +1 or -1).
114 *> 1 <= IFST <= N; 1 <= ILST <= N.
115 *> \endverbatim
116 *>
117 *> \param[out] WORK
118 *> \verbatim
119 *> WORK is REAL array, dimension (N)
120 *> \endverbatim
121 *>
122 *> \param[out] INFO
123 *> \verbatim
124 *> INFO is INTEGER
125 *> = 0: successful exit
126 *> < 0: if INFO = -i, the i-th argument had an illegal value
127 *> = 1: two adjacent blocks were too close to swap (the problem
128 *> is very ill-conditioned); T may have been partially
129 *> reordered, and ILST points to the first row of the
130 *> current position of the block being moved.
131 *> \endverbatim
132 *
133 * Authors:
134 * ========
135 *
136 *> \author Univ. of Tennessee
137 *> \author Univ. of California Berkeley
138 *> \author Univ. of Colorado Denver
139 *> \author NAG Ltd.
140 *
141 *> \date November 2011
142 *
143 *> \ingroup realOTHERcomputational
144 *
145 * =====================================================================
146  SUBROUTINE strexc( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
147  $ info )
148 *
149 * -- LAPACK computational routine (version 3.4.0) --
150 * -- LAPACK is a software package provided by Univ. of Tennessee, --
151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152 * November 2011
153 *
154 * .. Scalar Arguments ..
155  CHARACTER compq
156  INTEGER ifst, ilst, info, ldq, ldt, n
157 * ..
158 * .. Array Arguments ..
159  REAL q( ldq, * ), t( ldt, * ), work( * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * .. Parameters ..
165  REAL zero
166  parameter( zero = 0.0e+0 )
167 * ..
168 * .. Local Scalars ..
169  LOGICAL wantq
170  INTEGER here, nbf, nbl, nbnext
171 * ..
172 * .. External Functions ..
173  LOGICAL lsame
174  EXTERNAL lsame
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL slaexc, xerbla
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC max
181 * ..
182 * .. Executable Statements ..
183 *
184 * Decode and test the input arguments.
185 *
186  info = 0
187  wantq = lsame( compq, 'V' )
188  IF( .NOT.wantq .AND. .NOT.lsame( compq, 'N' ) ) THEN
189  info = -1
190  ELSE IF( n.LT.0 ) THEN
191  info = -2
192  ELSE IF( ldt.LT.max( 1, n ) ) THEN
193  info = -4
194  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.max( 1, n ) ) ) THEN
195  info = -6
196  ELSE IF( ifst.LT.1 .OR. ifst.GT.n ) THEN
197  info = -7
198  ELSE IF( ilst.LT.1 .OR. ilst.GT.n ) THEN
199  info = -8
200  END IF
201  IF( info.NE.0 ) THEN
202  CALL xerbla( 'STREXC', -info )
203  return
204  END IF
205 *
206 * Quick return if possible
207 *
208  IF( n.LE.1 )
209  $ return
210 *
211 * Determine the first row of specified block
212 * and find out it is 1 by 1 or 2 by 2.
213 *
214  IF( ifst.GT.1 ) THEN
215  IF( t( ifst, ifst-1 ).NE.zero )
216  $ ifst = ifst - 1
217  END IF
218  nbf = 1
219  IF( ifst.LT.n ) THEN
220  IF( t( ifst+1, ifst ).NE.zero )
221  $ nbf = 2
222  END IF
223 *
224 * Determine the first row of the final block
225 * and find out it is 1 by 1 or 2 by 2.
226 *
227  IF( ilst.GT.1 ) THEN
228  IF( t( ilst, ilst-1 ).NE.zero )
229  $ ilst = ilst - 1
230  END IF
231  nbl = 1
232  IF( ilst.LT.n ) THEN
233  IF( t( ilst+1, ilst ).NE.zero )
234  $ nbl = 2
235  END IF
236 *
237  IF( ifst.EQ.ilst )
238  $ return
239 *
240  IF( ifst.LT.ilst ) THEN
241 *
242 * Update ILST
243 *
244  IF( nbf.EQ.2 .AND. nbl.EQ.1 )
245  $ ilst = ilst - 1
246  IF( nbf.EQ.1 .AND. nbl.EQ.2 )
247  $ ilst = ilst + 1
248 *
249  here = ifst
250 *
251  10 continue
252 *
253 * Swap block with next one below
254 *
255  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
256 *
257 * Current block either 1 by 1 or 2 by 2
258 *
259  nbnext = 1
260  IF( here+nbf+1.LE.n ) THEN
261  IF( t( here+nbf+1, here+nbf ).NE.zero )
262  $ nbnext = 2
263  END IF
264  CALL slaexc( wantq, n, t, ldt, q, ldq, here, nbf, nbnext,
265  $ work, info )
266  IF( info.NE.0 ) THEN
267  ilst = here
268  return
269  END IF
270  here = here + nbnext
271 *
272 * Test if 2 by 2 block breaks into two 1 by 1 blocks
273 *
274  IF( nbf.EQ.2 ) THEN
275  IF( t( here+1, here ).EQ.zero )
276  $ nbf = 3
277  END IF
278 *
279  ELSE
280 *
281 * Current block consists of two 1 by 1 blocks each of which
282 * must be swapped individually
283 *
284  nbnext = 1
285  IF( here+3.LE.n ) THEN
286  IF( t( here+3, here+2 ).NE.zero )
287  $ nbnext = 2
288  END IF
289  CALL slaexc( wantq, n, t, ldt, q, ldq, here+1, 1, nbnext,
290  $ work, info )
291  IF( info.NE.0 ) THEN
292  ilst = here
293  return
294  END IF
295  IF( nbnext.EQ.1 ) THEN
296 *
297 * Swap two 1 by 1 blocks, no problems possible
298 *
299  CALL slaexc( wantq, n, t, ldt, q, ldq, here, 1, nbnext,
300  $ work, info )
301  here = here + 1
302  ELSE
303 *
304 * Recompute NBNEXT in case 2 by 2 split
305 *
306  IF( t( here+2, here+1 ).EQ.zero )
307  $ nbnext = 1
308  IF( nbnext.EQ.2 ) THEN
309 *
310 * 2 by 2 Block did not split
311 *
312  CALL slaexc( wantq, n, t, ldt, q, ldq, here, 1,
313  $ nbnext, work, info )
314  IF( info.NE.0 ) THEN
315  ilst = here
316  return
317  END IF
318  here = here + 2
319  ELSE
320 *
321 * 2 by 2 Block did split
322 *
323  CALL slaexc( wantq, n, t, ldt, q, ldq, here, 1, 1,
324  $ work, info )
325  CALL slaexc( wantq, n, t, ldt, q, ldq, here+1, 1, 1,
326  $ work, info )
327  here = here + 2
328  END IF
329  END IF
330  END IF
331  IF( here.LT.ilst )
332  $ go to 10
333 *
334  ELSE
335 *
336  here = ifst
337  20 continue
338 *
339 * Swap block with next one above
340 *
341  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
342 *
343 * Current block either 1 by 1 or 2 by 2
344 *
345  nbnext = 1
346  IF( here.GE.3 ) THEN
347  IF( t( here-1, here-2 ).NE.zero )
348  $ nbnext = 2
349  END IF
350  CALL slaexc( wantq, n, t, ldt, q, ldq, here-nbnext, nbnext,
351  $ nbf, work, info )
352  IF( info.NE.0 ) THEN
353  ilst = here
354  return
355  END IF
356  here = here - nbnext
357 *
358 * Test if 2 by 2 block breaks into two 1 by 1 blocks
359 *
360  IF( nbf.EQ.2 ) THEN
361  IF( t( here+1, here ).EQ.zero )
362  $ nbf = 3
363  END IF
364 *
365  ELSE
366 *
367 * Current block consists of two 1 by 1 blocks each of which
368 * must be swapped individually
369 *
370  nbnext = 1
371  IF( here.GE.3 ) THEN
372  IF( t( here-1, here-2 ).NE.zero )
373  $ nbnext = 2
374  END IF
375  CALL slaexc( wantq, n, t, ldt, q, ldq, here-nbnext, nbnext,
376  $ 1, work, info )
377  IF( info.NE.0 ) THEN
378  ilst = here
379  return
380  END IF
381  IF( nbnext.EQ.1 ) THEN
382 *
383 * Swap two 1 by 1 blocks, no problems possible
384 *
385  CALL slaexc( wantq, n, t, ldt, q, ldq, here, nbnext, 1,
386  $ work, info )
387  here = here - 1
388  ELSE
389 *
390 * Recompute NBNEXT in case 2 by 2 split
391 *
392  IF( t( here, here-1 ).EQ.zero )
393  $ nbnext = 1
394  IF( nbnext.EQ.2 ) THEN
395 *
396 * 2 by 2 Block did not split
397 *
398  CALL slaexc( wantq, n, t, ldt, q, ldq, here-1, 2, 1,
399  $ work, info )
400  IF( info.NE.0 ) THEN
401  ilst = here
402  return
403  END IF
404  here = here - 2
405  ELSE
406 *
407 * 2 by 2 Block did split
408 *
409  CALL slaexc( wantq, n, t, ldt, q, ldq, here, 1, 1,
410  $ work, info )
411  CALL slaexc( wantq, n, t, ldt, q, ldq, here-1, 1, 1,
412  $ work, info )
413  here = here - 2
414  END IF
415  END IF
416  END IF
417  IF( here.GT.ilst )
418  $ go to 20
419  END IF
420  ilst = here
421 *
422  return
423 *
424 * End of STREXC
425 *
426  END