LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slaexc()

subroutine slaexc ( logical wantq,
integer n,
real, dimension( ldt, * ) t,
integer ldt,
real, dimension( ldq, * ) q,
integer ldq,
integer j1,
integer n1,
integer n2,
real, dimension( * ) work,
integer info )

SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.

Download SLAEXC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
!> an upper quasi-triangular matrix T by an orthogonal similarity
!> transformation.
!>
!> T must be in Schur canonical form, that is, block upper triangular
!> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
!> has its diagonal elements equal and its off-diagonal elements of
!> opposite sign.
!>
Parameters
[in]WANTQ
!>          WANTQ is LOGICAL
!>          = .TRUE. : accumulate the transformation in the matrix Q;
!>          = .FALSE.: do not accumulate the transformation.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix T. N >= 0.
!>
[in,out]T
!>          T is REAL array, dimension (LDT,N)
!>          On entry, the upper quasi-triangular matrix T, in Schur
!>          canonical form.
!>          On exit, the updated matrix T, again in Schur canonical form.
!>
[in]LDT
!>          LDT is INTEGER
!>          The leading dimension of the array T. LDT >= max(1,N).
!>
[in,out]Q
!>          Q is REAL array, dimension (LDQ,N)
!>          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
!>          On exit, if WANTQ is .TRUE., the updated matrix Q.
!>          If WANTQ is .FALSE., Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.
!>          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
!>
[in]J1
!>          J1 is INTEGER
!>          The index of the first row of the first block T11.
!>
[in]N1
!>          N1 is INTEGER
!>          The order of the first block T11. N1 = 0, 1 or 2.
!>
[in]N2
!>          N2 is INTEGER
!>          The order of the second block T22. N2 = 0, 1 or 2.
!>
[out]WORK
!>          WORK is REAL array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          = 1: the transformed matrix T would be too far from Schur
!>               form; the blocks are not swapped and T and Q are
!>               unchanged.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 134 of file slaexc.f.

136*
137* -- LAPACK auxiliary routine --
138* -- LAPACK is a software package provided by Univ. of Tennessee, --
139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*
141* .. Scalar Arguments ..
142 LOGICAL WANTQ
143 INTEGER INFO, J1, LDQ, LDT, N, N1, N2
144* ..
145* .. Array Arguments ..
146 REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
147* ..
148*
149* =====================================================================
150*
151* .. Parameters ..
152 REAL ZERO, ONE
153 parameter( zero = 0.0e+0, one = 1.0e+0 )
154 REAL TEN
155 parameter( ten = 1.0e+1 )
156 INTEGER LDD, LDX
157 parameter( ldd = 4, ldx = 2 )
158* ..
159* .. Local Scalars ..
160 INTEGER IERR, J2, J3, J4, K, ND
161 REAL CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
162 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
163 $ WR1, WR2, XNORM
164* ..
165* .. Local Arrays ..
166 REAL D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
167 $ X( LDX, 2 )
168* ..
169* .. External Functions ..
170 REAL SLAMCH, SLANGE
171 EXTERNAL slamch, slange
172* ..
173* .. External Subroutines ..
174 EXTERNAL slacpy, slanv2, slarfg, slarfx, slartg,
175 $ slasy2,
176 $ srot
177* ..
178* .. Intrinsic Functions ..
179 INTRINSIC abs, max
180* ..
181* .. Executable Statements ..
182*
183 info = 0
184*
185* Quick return if possible
186*
187 IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )
188 $ RETURN
189 IF( j1+n1.GT.n )
190 $ RETURN
191*
192 j2 = j1 + 1
193 j3 = j1 + 2
194 j4 = j1 + 3
195*
196 IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN
197*
198* Swap two 1-by-1 blocks.
199*
200 t11 = t( j1, j1 )
201 t22 = t( j2, j2 )
202*
203* Determine the transformation to perform the interchange.
204*
205 CALL slartg( t( j1, j2 ), t22-t11, cs, sn, temp )
206*
207* Apply transformation to the matrix T.
208*
209 IF( j3.LE.n )
210 $ CALL srot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt,
211 $ 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,
280 $ work )
281 CALL slarfx( 'R', j2, 3, u, tau, t( 1, j1 ), ldt, work )
282*
283 t( j3, j1 ) = zero
284 t( j3, j2 ) = zero
285 t( j3, j3 ) = t11
286*
287 IF( wantq ) THEN
288*
289* Accumulate transformation in the matrix Q.
290*
291 CALL slarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
292 END IF
293 GO TO 40
294*
295 20 CONTINUE
296*
297* N1 = 2, N2 = 1: generate elementary reflector H so that:
298*
299* H ( -X11 ) = ( * )
300* ( -X21 ) = ( 0 )
301* ( scale ) = ( 0 )
302*
303 u( 1 ) = -x( 1, 1 )
304 u( 2 ) = -x( 2, 1 )
305 u( 3 ) = scale
306 CALL slarfg( 3, u( 1 ), u( 2 ), 1, tau )
307 u( 1 ) = one
308 t33 = t( j3, j3 )
309*
310* Perform swap provisionally on diagonal block in D.
311*
312 CALL slarfx( 'L', 3, 3, u, tau, d, ldd, work )
313 CALL slarfx( 'R', 3, 3, u, tau, d, ldd, work )
314*
315* Test whether to reject swap.
316*
317 IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,
318 $ 1 )-t33 ) ).GT.thresh )GO TO 50
319*
320* Accept swap: apply transformation to the entire matrix T.
321*
322 CALL slarfx( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
323 CALL slarfx( 'L', 3, n-j1, u, tau, t( j1, j2 ), ldt, work )
324*
325 t( j1, j1 ) = t33
326 t( j2, j1 ) = zero
327 t( j3, j1 ) = zero
328*
329 IF( wantq ) THEN
330*
331* Accumulate transformation in the matrix Q.
332*
333 CALL slarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
334 END IF
335 GO TO 40
336*
337 30 CONTINUE
338*
339* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
340* that:
341*
342* H(2) H(1) ( -X11 -X12 ) = ( * * )
343* ( -X21 -X22 ) ( 0 * )
344* ( scale 0 ) ( 0 0 )
345* ( 0 scale ) ( 0 0 )
346*
347 u1( 1 ) = -x( 1, 1 )
348 u1( 2 ) = -x( 2, 1 )
349 u1( 3 ) = scale
350 CALL slarfg( 3, u1( 1 ), u1( 2 ), 1, tau1 )
351 u1( 1 ) = one
352*
353 temp = -tau1*( x( 1, 2 )+u1( 2 )*x( 2, 2 ) )
354 u2( 1 ) = -temp*u1( 2 ) - x( 2, 2 )
355 u2( 2 ) = -temp*u1( 3 )
356 u2( 3 ) = scale
357 CALL slarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
358 u2( 1 ) = one
359*
360* Perform swap provisionally on diagonal block in D.
361*
362 CALL slarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
363 CALL slarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
364 CALL slarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
365 CALL slarfx( 'R', 4, 3, u2, tau2, d( 1, 2 ), ldd, work )
366*
367* Test whether to reject swap.
368*
369 IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 4, 1 ) ),
370 $ abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50
371*
372* Accept swap: apply transformation to the entire matrix T.
373*
374 CALL slarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt,
375 $ work )
376 CALL slarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
377 CALL slarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt,
378 $ work )
379 CALL slarfx( 'R', j4, 3, u2, tau2, t( 1, j2 ), ldt, work )
380*
381 t( j3, j1 ) = zero
382 t( j3, j2 ) = zero
383 t( j4, j1 ) = zero
384 t( j4, j2 ) = zero
385*
386 IF( wantq ) THEN
387*
388* Accumulate transformation in the matrix Q.
389*
390 CALL slarfx( 'R', n, 3, u1, tau1, q( 1, j1 ), ldq, work )
391 CALL slarfx( 'R', n, 3, u2, tau2, q( 1, j2 ), ldq, work )
392 END IF
393*
394 40 CONTINUE
395*
396 IF( n2.EQ.2 ) THEN
397*
398* Standardize new 2-by-2 block T11
399*
400 CALL slanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),
401 $ t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )
402 CALL srot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ),
403 $ ldt,
404 $ cs, sn )
405 CALL srot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
406 IF( wantq )
407 $ CALL srot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
408 END IF
409*
410 IF( n1.EQ.2 ) THEN
411*
412* Standardize new 2-by-2 block T22
413*
414 j3 = j1 + n2
415 j4 = j3 + 1
416 CALL slanv2( t( j3, j3 ), t( j3, j4 ), t( j4, j3 ),
417 $ t( j4, j4 ), wr1, wi1, wr2, wi2, cs, sn )
418 IF( j3+2.LE.n )
419 $ CALL srot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
420 $ ldt, cs, sn )
421 CALL srot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
422 IF( wantq )
423 $ CALL srot( n, q( 1, j3 ), 1, q( 1, j4 ), 1, cs, sn )
424 END IF
425*
426 END IF
427 RETURN
428*
429* Exit with INFO = 1 if swap was rejected.
430*
431 50 info = 1
432 RETURN
433*
434* End of SLAEXC
435*
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:112
slanv2
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:125
slarfg
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:104
slarfx
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:118
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
slasy2
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:172
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: