LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 136 of file slaexc.f.

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 *
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 slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:113
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:114
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 slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
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 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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: