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

◆ 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*
slacpy
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
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
slanv2
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
slarfg
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
slarfx
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
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
slasy2
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
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: