LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slaexc.f
Go to the documentation of this file.
1*> \brief \b SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLAEXC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaexc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaexc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaexc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* LOGICAL WANTQ
26* INTEGER INFO, J1, LDQ, LDT, N, N1, N2
27* ..
28* .. Array Arguments ..
29* REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
39*> an upper quasi-triangular matrix T by an orthogonal similarity
40*> transformation.
41*>
42*> T must be in Schur canonical form, that is, block upper triangular
43*> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
44*> has its diagonal elements equal and its off-diagonal elements of
45*> opposite sign.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] WANTQ
52*> \verbatim
53*> WANTQ is LOGICAL
54*> = .TRUE. : accumulate the transformation in the matrix Q;
55*> = .FALSE.: do not accumulate the transformation.
56*> \endverbatim
57*>
58*> \param[in] N
59*> \verbatim
60*> N is INTEGER
61*> The order of the matrix T. N >= 0.
62*> \endverbatim
63*>
64*> \param[in,out] T
65*> \verbatim
66*> T is REAL array, dimension (LDT,N)
67*> On entry, the upper quasi-triangular matrix T, in Schur
68*> canonical form.
69*> On exit, the updated matrix T, again in Schur canonical form.
70*> \endverbatim
71*>
72*> \param[in] LDT
73*> \verbatim
74*> LDT is INTEGER
75*> The leading dimension of the array T. LDT >= max(1,N).
76*> \endverbatim
77*>
78*> \param[in,out] Q
79*> \verbatim
80*> Q is REAL array, dimension (LDQ,N)
81*> On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
82*> On exit, if WANTQ is .TRUE., the updated matrix Q.
83*> If WANTQ is .FALSE., Q is not referenced.
84*> \endverbatim
85*>
86*> \param[in] LDQ
87*> \verbatim
88*> LDQ is INTEGER
89*> The leading dimension of the array Q.
90*> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
91*> \endverbatim
92*>
93*> \param[in] J1
94*> \verbatim
95*> J1 is INTEGER
96*> The index of the first row of the first block T11.
97*> \endverbatim
98*>
99*> \param[in] N1
100*> \verbatim
101*> N1 is INTEGER
102*> The order of the first block T11. N1 = 0, 1 or 2.
103*> \endverbatim
104*>
105*> \param[in] N2
106*> \verbatim
107*> N2 is INTEGER
108*> The order of the second block T22. N2 = 0, 1 or 2.
109*> \endverbatim
110*>
111*> \param[out] WORK
112*> \verbatim
113*> WORK is REAL array, dimension (N)
114*> \endverbatim
115*>
116*> \param[out] INFO
117*> \verbatim
118*> INFO is INTEGER
119*> = 0: successful exit
120*> = 1: the transformed matrix T would be too far from Schur
121*> form; the blocks are not swapped and T and Q are
122*> unchanged.
123*> \endverbatim
124*
125* Authors:
126* ========
127*
128*> \author Univ. of Tennessee
129*> \author Univ. of California Berkeley
130*> \author Univ. of Colorado Denver
131*> \author NAG Ltd.
132*
133*> \ingroup laexc
134*
135* =====================================================================
136 SUBROUTINE slaexc( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
137 $ INFO )
138*
139* -- LAPACK auxiliary routine --
140* -- LAPACK is a software package provided by Univ. of Tennessee, --
141* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142*
143* .. Scalar Arguments ..
144 LOGICAL WANTQ
145 INTEGER INFO, J1, LDQ, LDT, N, N1, N2
146* ..
147* .. Array Arguments ..
148 REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
149* ..
150*
151* =====================================================================
152*
153* .. Parameters ..
154 REAL ZERO, ONE
155 parameter( zero = 0.0e+0, one = 1.0e+0 )
156 REAL TEN
157 parameter( ten = 1.0e+1 )
158 INTEGER LDD, LDX
159 parameter( ldd = 4, ldx = 2 )
160* ..
161* .. Local Scalars ..
162 INTEGER IERR, J2, J3, J4, K, ND
163 REAL CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
164 $ t33, tau, tau1, tau2, temp, thresh, wi1, wi2,
165 $ wr1, wr2, xnorm
166* ..
167* .. Local Arrays ..
168 REAL D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
169 $ x( ldx, 2 )
170* ..
171* .. External Functions ..
172 REAL SLAMCH, SLANGE
173 EXTERNAL slamch, slange
174* ..
175* .. External Subroutines ..
176 EXTERNAL slacpy, slanv2, slarfg, slarfx, slartg, slasy2,
177 $ srot
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC abs, max
181* ..
182* .. Executable Statements ..
183*
184 info = 0
185*
186* Quick return if possible
187*
188 IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )
189 $ RETURN
190 IF( j1+n1.GT.n )
191 $ RETURN
192*
193 j2 = j1 + 1
194 j3 = j1 + 2
195 j4 = j1 + 3
196*
197 IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN
198*
199* Swap two 1-by-1 blocks.
200*
201 t11 = t( j1, j1 )
202 t22 = t( j2, j2 )
203*
204* Determine the transformation to perform the interchange.
205*
206 CALL slartg( t( j1, j2 ), t22-t11, cs, sn, temp )
207*
208* Apply transformation to the matrix T.
209*
210 IF( j3.LE.n )
211 $ CALL srot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,
212 $ sn )
213 CALL srot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
214*
215 t( j1, j1 ) = t22
216 t( j2, j2 ) = t11
217*
218 IF( wantq ) THEN
219*
220* Accumulate transformation in the matrix Q.
221*
222 CALL srot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
223 END IF
224*
225 ELSE
226*
227* Swapping involves at least one 2-by-2 block.
228*
229* Copy the diagonal block of order N1+N2 to the local array D
230* and compute its norm.
231*
232 nd = n1 + n2
233 CALL slacpy( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
234 dnorm = slange( 'Max', nd, nd, d, ldd, work )
235*
236* Compute machine-dependent threshold for test for accepting
237* swap.
238*
239 eps = slamch( 'P' )
240 smlnum = slamch( 'S' ) / eps
241 thresh = max( ten*eps*dnorm, smlnum )
242*
243* Solve T11*X - X*T22 = scale*T12 for X.
244*
245 CALL slasy2( .false., .false., -1, n1, n2, d, ldd,
246 $ d( n1+1, n1+1 ), ldd, d( 1, n1+1 ), ldd, scale, x,
247 $ ldx, xnorm, ierr )
248*
249* Swap the adjacent diagonal blocks.
250*
251 k = n1 + n1 + n2 - 3
252 GO TO ( 10, 20, 30 )k
253*
254 10 CONTINUE
255*
256* N1 = 1, N2 = 2: generate elementary reflector H so that:
257*
258* ( scale, X11, X12 ) H = ( 0, 0, * )
259*
260 u( 1 ) = scale
261 u( 2 ) = x( 1, 1 )
262 u( 3 ) = x( 1, 2 )
263 CALL slarfg( 3, u( 3 ), u, 1, tau )
264 u( 3 ) = one
265 t11 = t( j1, j1 )
266*
267* Perform swap provisionally on diagonal block in D.
268*
269 CALL slarfx( 'L', 3, 3, u, tau, d, ldd, work )
270 CALL slarfx( 'R', 3, 3, u, tau, d, ldd, work )
271*
272* Test whether to reject swap.
273*
274 IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 3,
275 $ 3 )-t11 ) ).GT.thresh )GO TO 50
276*
277* Accept swap: apply transformation to the entire matrix T.
278*
279 CALL slarfx( 'L', 3, n-j1+1, u, tau, t( j1, j1 ), ldt, work )
280 CALL slarfx( 'R', j2, 3, u, tau, t( 1, j1 ), ldt, work )
281*
282 t( j3, j1 ) = zero
283 t( j3, j2 ) = zero
284 t( j3, j3 ) = t11
285*
286 IF( wantq ) THEN
287*
288* Accumulate transformation in the matrix Q.
289*
290 CALL slarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
291 END IF
292 GO TO 40
293*
294 20 CONTINUE
295*
296* N1 = 2, N2 = 1: generate elementary reflector H so that:
297*
298* H ( -X11 ) = ( * )
299* ( -X21 ) = ( 0 )
300* ( scale ) = ( 0 )
301*
302 u( 1 ) = -x( 1, 1 )
303 u( 2 ) = -x( 2, 1 )
304 u( 3 ) = scale
305 CALL slarfg( 3, u( 1 ), u( 2 ), 1, tau )
306 u( 1 ) = one
307 t33 = t( j3, j3 )
308*
309* Perform swap provisionally on diagonal block in D.
310*
311 CALL slarfx( 'L', 3, 3, u, tau, d, ldd, work )
312 CALL slarfx( 'R', 3, 3, u, tau, d, ldd, work )
313*
314* Test whether to reject swap.
315*
316 IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,
317 $ 1 )-t33 ) ).GT.thresh )GO TO 50
318*
319* Accept swap: apply transformation to the entire matrix T.
320*
321 CALL slarfx( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
322 CALL slarfx( 'L', 3, n-j1, u, tau, t( j1, j2 ), ldt, work )
323*
324 t( j1, j1 ) = t33
325 t( j2, j1 ) = zero
326 t( j3, j1 ) = zero
327*
328 IF( wantq ) THEN
329*
330* Accumulate transformation in the matrix Q.
331*
332 CALL slarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
333 END IF
334 GO TO 40
335*
336 30 CONTINUE
337*
338* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
339* that:
340*
341* H(2) H(1) ( -X11 -X12 ) = ( * * )
342* ( -X21 -X22 ) ( 0 * )
343* ( scale 0 ) ( 0 0 )
344* ( 0 scale ) ( 0 0 )
345*
346 u1( 1 ) = -x( 1, 1 )
347 u1( 2 ) = -x( 2, 1 )
348 u1( 3 ) = scale
349 CALL slarfg( 3, u1( 1 ), u1( 2 ), 1, tau1 )
350 u1( 1 ) = one
351*
352 temp = -tau1*( x( 1, 2 )+u1( 2 )*x( 2, 2 ) )
353 u2( 1 ) = -temp*u1( 2 ) - x( 2, 2 )
354 u2( 2 ) = -temp*u1( 3 )
355 u2( 3 ) = scale
356 CALL slarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
357 u2( 1 ) = one
358*
359* Perform swap provisionally on diagonal block in D.
360*
361 CALL slarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
362 CALL slarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
363 CALL slarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
364 CALL slarfx( 'R', 4, 3, u2, tau2, d( 1, 2 ), ldd, work )
365*
366* Test whether to reject swap.
367*
368 IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 4, 1 ) ),
369 $ abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50
370*
371* Accept swap: apply transformation to the entire matrix T.
372*
373 CALL slarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt, work )
374 CALL slarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
375 CALL slarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt, work )
376 CALL slarfx( 'R', j4, 3, u2, tau2, t( 1, j2 ), ldt, work )
377*
378 t( j3, j1 ) = zero
379 t( j3, j2 ) = zero
380 t( j4, j1 ) = zero
381 t( j4, j2 ) = zero
382*
383 IF( wantq ) THEN
384*
385* Accumulate transformation in the matrix Q.
386*
387 CALL slarfx( 'R', n, 3, u1, tau1, q( 1, j1 ), ldq, work )
388 CALL slarfx( 'R', n, 3, u2, tau2, q( 1, j2 ), ldq, work )
389 END IF
390*
391 40 CONTINUE
392*
393 IF( n2.EQ.2 ) THEN
394*
395* Standardize new 2-by-2 block T11
396*
397 CALL slanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),
398 $ t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )
399 CALL srot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ), ldt,
400 $ cs, sn )
401 CALL srot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
402 IF( wantq )
403 $ CALL srot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
404 END IF
405*
406 IF( n1.EQ.2 ) THEN
407*
408* Standardize new 2-by-2 block T22
409*
410 j3 = j1 + n2
411 j4 = j3 + 1
412 CALL slanv2( t( j3, j3 ), t( j3, j4 ), t( j4, j3 ),
413 $ t( j4, j4 ), wr1, wi1, wr2, wi2, cs, sn )
414 IF( j3+2.LE.n )
415 $ CALL srot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
416 $ ldt, cs, sn )
417 CALL srot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
418 IF( wantq )
419 $ CALL srot( n, q( 1, j3 ), 1, q( 1, j4 ), 1, cs, sn )
420 END IF
421*
422 END IF
423 RETURN
424*
425* Exit with INFO = 1 if swap was rejected.
426*
427 50 info = 1
428 RETURN
429*
430* End of SLAEXC
431*
432 END
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine slaexc(wantq, n, t, ldt, q, ldq, j1, n1, n2, work, info)
SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form...
Definition slaexc.f:138
subroutine slanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition slanv2.f:127
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
subroutine slarfx(side, m, n, v, tau, c, ldc, work)
SLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the ...
Definition slarfx.f:120
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
subroutine slasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition slasy2.f:174
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92