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

◆ 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 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 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.

Definition at line 134 of file dlaexc.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 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
147* ..
148*
149* =====================================================================
150*
151* .. Parameters ..
152 DOUBLE PRECISION ZERO, ONE
153 parameter( zero = 0.0d+0, one = 1.0d+0 )
154 DOUBLE PRECISION TEN
155 parameter( ten = 1.0d+1 )
156 INTEGER LDD, LDX
157 parameter( ldd = 4, ldx = 2 )
158* ..
159* .. Local Scalars ..
160 INTEGER IERR, J2, J3, J4, K, ND
161 DOUBLE PRECISION 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 DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
167 $ X( LDX, 2 )
168* ..
169* .. External Functions ..
170 DOUBLE PRECISION DLAMCH, DLANGE
171 EXTERNAL dlamch, dlange
172* ..
173* .. External Subroutines ..
174 EXTERNAL dlacpy, dlanv2, dlarfg, dlarfx, dlartg,
175 $ dlasy2,
176 $ drot
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 dlartg( t( j1, j2 ), t22-t11, cs, sn, temp )
206*
207* Apply transformation to the matrix T.
208*
209 IF( j3.LE.n )
210 $ CALL drot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt,
211 $ cs,
212 $ sn )
213 CALL drot( 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 drot( 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 dlacpy( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
234 dnorm = dlange( 'Max', nd, nd, d, ldd, work )
235*
236* Compute machine-dependent threshold for test for accepting
237* swap.
238*
239 eps = dlamch( 'P' )
240 smlnum = dlamch( 'S' ) / eps
241 thresh = max( ten*eps*dnorm, smlnum )
242*
243* Solve T11*X - X*T22 = scale*T12 for X.
244*
245 CALL dlasy2( .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 dlarfg( 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 dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
270 CALL dlarfx( '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 dlarfx( 'L', 3, n-j1+1, u, tau, t( j1, j1 ), ldt,
280 $ work )
281 CALL dlarfx( '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 dlarfx( '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 dlarfg( 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 dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
313 CALL dlarfx( '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 dlarfx( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
323 CALL dlarfx( '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 dlarfx( '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 dlarfg( 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 dlarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
358 u2( 1 ) = one
359*
360* Perform swap provisionally on diagonal block in D.
361*
362 CALL dlarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
363 CALL dlarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
364 CALL dlarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
365 CALL dlarfx( '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 dlarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt,
375 $ work )
376 CALL dlarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
377 CALL dlarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt,
378 $ 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 ),
403 $ ldt,
404 $ cs, sn )
405 CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
406 IF( wantq )
407 $ CALL drot( 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 dlanv2( 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 drot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
420 $ ldt, cs, sn )
421 CALL drot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
422 IF( wantq )
423 $ CALL drot( 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 CONTINUE
432 info = 1
433 RETURN
434*
435* End of DLAEXC
436*
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:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:112
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:125
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:104
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:118
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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:172
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: