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