LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlaexc ( logical  WANTQ,
integer  N,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
integer  J1,
integer  N1,
integer  N2,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose: 
 DLAEXC 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 elemnts 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.
Date
September 2012

Definition at line 140 of file dlaexc.f.

140 *
141 * -- LAPACK auxiliary routine (version 3.4.2) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * September 2012
145 *
146 * .. Scalar Arguments ..
147  LOGICAL wantq
148  INTEGER info, j1, ldq, ldt, n, n1, n2
149 * ..
150 * .. Array Arguments ..
151  DOUBLE PRECISION q( ldq, * ), t( ldt, * ), work( * )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Parameters ..
157  DOUBLE PRECISION zero, one
158  parameter ( zero = 0.0d+0, one = 1.0d+0 )
159  DOUBLE PRECISION ten
160  parameter ( ten = 1.0d+1 )
161  INTEGER ldd, ldx
162  parameter ( ldd = 4, ldx = 2 )
163 * ..
164 * .. Local Scalars ..
165  INTEGER ierr, j2, j3, j4, k, nd
166  DOUBLE PRECISION cs, dnorm, eps, scale, smlnum, sn, t11, t22,
167  $ t33, tau, tau1, tau2, temp, thresh, wi1, wi2,
168  $ wr1, wr2, xnorm
169 * ..
170 * .. Local Arrays ..
171  DOUBLE PRECISION d( ldd, 4 ), u( 3 ), u1( 3 ), u2( 3 ),
172  $ x( ldx, 2 )
173 * ..
174 * .. External Functions ..
175  DOUBLE PRECISION dlamch, dlange
176  EXTERNAL dlamch, dlange
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2,
180  $ drot
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max
184 * ..
185 * .. Executable Statements ..
186 *
187  info = 0
188 *
189 * Quick return if possible
190 *
191  IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )
192  $ RETURN
193  IF( j1+n1.GT.n )
194  $ RETURN
195 *
196  j2 = j1 + 1
197  j3 = j1 + 2
198  j4 = j1 + 3
199 *
200  IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN
201 *
202 * Swap two 1-by-1 blocks.
203 *
204  t11 = t( j1, j1 )
205  t22 = t( j2, j2 )
206 *
207 * Determine the transformation to perform the interchange.
208 *
209  CALL dlartg( t( j1, j2 ), t22-t11, cs, sn, temp )
210 *
211 * Apply transformation to the matrix T.
212 *
213  IF( j3.LE.n )
214  $ CALL drot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,
215  $ sn )
216  CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
217 *
218  t( j1, j1 ) = t22
219  t( j2, j2 ) = t11
220 *
221  IF( wantq ) THEN
222 *
223 * Accumulate transformation in the matrix Q.
224 *
225  CALL drot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
226  END IF
227 *
228  ELSE
229 *
230 * Swapping involves at least one 2-by-2 block.
231 *
232 * Copy the diagonal block of order N1+N2 to the local array D
233 * and compute its norm.
234 *
235  nd = n1 + n2
236  CALL dlacpy( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
237  dnorm = dlange( 'Max', nd, nd, d, ldd, work )
238 *
239 * Compute machine-dependent threshold for test for accepting
240 * swap.
241 *
242  eps = dlamch( 'P' )
243  smlnum = dlamch( 'S' ) / eps
244  thresh = max( ten*eps*dnorm, smlnum )
245 *
246 * Solve T11*X - X*T22 = scale*T12 for X.
247 *
248  CALL dlasy2( .false., .false., -1, n1, n2, d, ldd,
249  $ d( n1+1, n1+1 ), ldd, d( 1, n1+1 ), ldd, scale, x,
250  $ ldx, xnorm, ierr )
251 *
252 * Swap the adjacent diagonal blocks.
253 *
254  k = n1 + n1 + n2 - 3
255  GO TO ( 10, 20, 30 )k
256 *
257  10 CONTINUE
258 *
259 * N1 = 1, N2 = 2: generate elementary reflector H so that:
260 *
261 * ( scale, X11, X12 ) H = ( 0, 0, * )
262 *
263  u( 1 ) = scale
264  u( 2 ) = x( 1, 1 )
265  u( 3 ) = x( 1, 2 )
266  CALL dlarfg( 3, u( 3 ), u, 1, tau )
267  u( 3 ) = one
268  t11 = t( j1, j1 )
269 *
270 * Perform swap provisionally on diagonal block in D.
271 *
272  CALL dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
273  CALL dlarfx( 'R', 3, 3, u, tau, d, ldd, work )
274 *
275 * Test whether to reject swap.
276 *
277  IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 3,
278  $ 3 )-t11 ) ).GT.thresh )GO TO 50
279 *
280 * Accept swap: apply transformation to the entire matrix T.
281 *
282  CALL dlarfx( 'L', 3, n-j1+1, u, tau, t( j1, j1 ), ldt, work )
283  CALL dlarfx( 'R', j2, 3, u, tau, t( 1, j1 ), ldt, work )
284 *
285  t( j3, j1 ) = zero
286  t( j3, j2 ) = zero
287  t( j3, j3 ) = t11
288 *
289  IF( wantq ) THEN
290 *
291 * Accumulate transformation in the matrix Q.
292 *
293  CALL dlarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
294  END IF
295  GO TO 40
296 *
297  20 CONTINUE
298 *
299 * N1 = 2, N2 = 1: generate elementary reflector H so that:
300 *
301 * H ( -X11 ) = ( * )
302 * ( -X21 ) = ( 0 )
303 * ( scale ) = ( 0 )
304 *
305  u( 1 ) = -x( 1, 1 )
306  u( 2 ) = -x( 2, 1 )
307  u( 3 ) = scale
308  CALL dlarfg( 3, u( 1 ), u( 2 ), 1, tau )
309  u( 1 ) = one
310  t33 = t( j3, j3 )
311 *
312 * Perform swap provisionally on diagonal block in D.
313 *
314  CALL dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
315  CALL dlarfx( 'R', 3, 3, u, tau, d, ldd, work )
316 *
317 * Test whether to reject swap.
318 *
319  IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,
320  $ 1 )-t33 ) ).GT.thresh )GO TO 50
321 *
322 * Accept swap: apply transformation to the entire matrix T.
323 *
324  CALL dlarfx( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
325  CALL dlarfx( 'L', 3, n-j1, u, tau, t( j1, j2 ), ldt, work )
326 *
327  t( j1, j1 ) = t33
328  t( j2, j1 ) = zero
329  t( j3, j1 ) = zero
330 *
331  IF( wantq ) THEN
332 *
333 * Accumulate transformation in the matrix Q.
334 *
335  CALL dlarfx( 'R', n, 3, u, tau, q( 1, j1 ), ldq, work )
336  END IF
337  GO TO 40
338 *
339  30 CONTINUE
340 *
341 * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
342 * that:
343 *
344 * H(2) H(1) ( -X11 -X12 ) = ( * * )
345 * ( -X21 -X22 ) ( 0 * )
346 * ( scale 0 ) ( 0 0 )
347 * ( 0 scale ) ( 0 0 )
348 *
349  u1( 1 ) = -x( 1, 1 )
350  u1( 2 ) = -x( 2, 1 )
351  u1( 3 ) = scale
352  CALL dlarfg( 3, u1( 1 ), u1( 2 ), 1, tau1 )
353  u1( 1 ) = one
354 *
355  temp = -tau1*( x( 1, 2 )+u1( 2 )*x( 2, 2 ) )
356  u2( 1 ) = -temp*u1( 2 ) - x( 2, 2 )
357  u2( 2 ) = -temp*u1( 3 )
358  u2( 3 ) = scale
359  CALL dlarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
360  u2( 1 ) = one
361 *
362 * Perform swap provisionally on diagonal block in D.
363 *
364  CALL dlarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
365  CALL dlarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
366  CALL dlarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
367  CALL dlarfx( 'R', 4, 3, u2, tau2, d( 1, 2 ), ldd, work )
368 *
369 * Test whether to reject swap.
370 *
371  IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 4, 1 ) ),
372  $ abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50
373 *
374 * Accept swap: apply transformation to the entire matrix T.
375 *
376  CALL dlarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt, work )
377  CALL dlarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
378  CALL dlarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt, work )
379  CALL dlarfx( '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 dlarfx( 'R', n, 3, u1, tau1, q( 1, j1 ), ldq, work )
391  CALL dlarfx( '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 dlanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),
401  $ t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )
402  CALL drot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ), ldt,
403  $ cs, sn )
404  CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
405  IF( wantq )
406  $ CALL drot( n, q( 1, j1 ), 1, q( 1, j2 ), 1, cs, sn )
407  END IF
408 *
409  IF( n1.EQ.2 ) THEN
410 *
411 * Standardize new 2-by-2 block T22
412 *
413  j3 = j1 + n2
414  j4 = j3 + 1
415  CALL dlanv2( t( j3, j3 ), t( j3, j4 ), t( j4, j3 ),
416  $ t( j4, j4 ), wr1, wi1, wr2, wi2, cs, sn )
417  IF( j3+2.LE.n )
418  $ CALL drot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
419  $ ldt, cs, sn )
420  CALL drot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
421  IF( wantq )
422  $ CALL drot( n, q( 1, j3 ), 1, q( 1, j4 ), 1, cs, sn )
423  END IF
424 *
425  END IF
426  RETURN
427 *
428 * Exit with INFO = 1 if swap was rejected.
429 *
430  50 CONTINUE
431  info = 1
432  RETURN
433 *
434 * End of DLAEXC
435 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
dlasy2
subroutine dlasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition: dlasy2.f:176
dlarfx
subroutine dlarfx(SIDE, M, N, V, TAU, C, LDC, WORK)
DLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the ...
Definition: dlarfx.f:122
dlange
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
dlarfg
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
dlanv2
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form...
Definition: dlanv2.f:129

Here is the call graph for this function:

Here is the caller graph for this function: