LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dlaexc()

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
December 2016

Definition at line 140 of file dlaexc.f.

140 *
141 * -- LAPACK auxiliary routine (version 3.7.0) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * December 2016
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 *
Here is the call graph for this function:
Here is the caller graph for this function:
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
dlarfg
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
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
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
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
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:94
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70
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
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