LAPACK 3.12.0
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 136 of file dlaexc.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 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
149* ..
150*
151* =====================================================================
152*
153* .. Parameters ..
154 DOUBLE PRECISION ZERO, ONE
155 parameter( zero = 0.0d+0, one = 1.0d+0 )
156 DOUBLE PRECISION TEN
157 parameter( ten = 1.0d+1 )
158 INTEGER LDD, LDX
159 parameter( ldd = 4, ldx = 2 )
160* ..
161* .. Local Scalars ..
162 INTEGER IERR, J2, J3, J4, K, ND
163 DOUBLE PRECISION 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 DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
169 $ X( LDX, 2 )
170* ..
171* .. External Functions ..
172 DOUBLE PRECISION DLAMCH, DLANGE
173 EXTERNAL dlamch, dlange
174* ..
175* .. External Subroutines ..
176 EXTERNAL dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2,
177 $ drot
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 dlartg( t( j1, j2 ), t22-t11, cs, sn, temp )
207*
208* Apply transformation to the matrix T.
209*
210 IF( j3.LE.n )
211 $ CALL drot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, 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, work )
280 CALL dlarfx( '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 dlarfx( '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 dlarfg( 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 dlarfx( 'L', 3, 3, u, tau, d, ldd, work )
312 CALL dlarfx( '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 dlarfx( 'R', j3, 3, u, tau, t( 1, j1 ), ldt, work )
322 CALL dlarfx( '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 dlarfx( '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 dlarfg( 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 dlarfg( 3, u2( 1 ), u2( 2 ), 1, tau2 )
357 u2( 1 ) = one
358*
359* Perform swap provisionally on diagonal block in D.
360*
361 CALL dlarfx( 'L', 3, 4, u1, tau1, d, ldd, work )
362 CALL dlarfx( 'R', 4, 3, u1, tau1, d, ldd, work )
363 CALL dlarfx( 'L', 3, 4, u2, tau2, d( 2, 1 ), ldd, work )
364 CALL dlarfx( '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 dlarfx( 'L', 3, n-j1+1, u1, tau1, t( j1, j1 ), ldt, work )
374 CALL dlarfx( 'R', j4, 3, u1, tau1, t( 1, j1 ), ldt, work )
375 CALL dlarfx( 'L', 3, n-j1+1, u2, tau2, t( j2, j1 ), ldt, work )
376 CALL dlarfx( '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 dlarfx( 'R', n, 3, u1, tau1, q( 1, j1 ), ldq, work )
388 CALL dlarfx( '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 dlanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),
398 $ t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )
399 CALL drot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ), ldt,
400 $ cs, sn )
401 CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
402 IF( wantq )
403 $ CALL drot( 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 dlanv2( 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 drot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
416 $ ldt, cs, sn )
417 CALL drot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
418 IF( wantq )
419 $ CALL drot( 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 CONTINUE
428 info = 1
429 RETURN
430*
431* End of DLAEXC
432*
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:103
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:114
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:127
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
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:120
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:174
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: