summaryrefslogtreecommitdiff
path: root/2.3-1/src/fortran/lapack/dlaexc.f
diff options
context:
space:
mode:
Diffstat (limited to '2.3-1/src/fortran/lapack/dlaexc.f')
-rw-r--r--2.3-1/src/fortran/lapack/dlaexc.f354
1 files changed, 354 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dlaexc.f b/2.3-1/src/fortran/lapack/dlaexc.f
new file mode 100644
index 00000000..18e7d247
--- /dev/null
+++ b/2.3-1/src/fortran/lapack/dlaexc.f
@@ -0,0 +1,354 @@
+ SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
+ $ INFO )
+*
+* -- LAPACK auxiliary routine (version 3.1) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ LOGICAL WANTQ
+ INTEGER INFO, J1, LDQ, LDT, N, N1, N2
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
+* ..
+*
+* 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.
+*
+* Arguments
+* =========
+*
+* WANTQ (input) LOGICAL
+* = .TRUE. : accumulate the transformation in the matrix Q;
+* = .FALSE.: do not accumulate the transformation.
+*
+* N (input) INTEGER
+* The order of the matrix T. N >= 0.
+*
+* T (input/output) 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.
+*
+* LDT (input) INTEGER
+* The leading dimension of the array T. LDT >= max(1,N).
+*
+* Q (input/output) 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.
+*
+* LDQ (input) INTEGER
+* The leading dimension of the array Q.
+* LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
+*
+* J1 (input) INTEGER
+* The index of the first row of the first block T11.
+*
+* N1 (input) INTEGER
+* The order of the first block T11. N1 = 0, 1 or 2.
+*
+* N2 (input) INTEGER
+* The order of the second block T22. N2 = 0, 1 or 2.
+*
+* WORK (workspace) DOUBLE PRECISION array, dimension (N)
+*
+* INFO (output) 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.
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ DOUBLE PRECISION TEN
+ PARAMETER ( TEN = 1.0D+1 )
+ INTEGER LDD, LDX
+ PARAMETER ( LDD = 4, LDX = 2 )
+* ..
+* .. Local Scalars ..
+ INTEGER IERR, J2, J3, J4, K, ND
+ DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
+ $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
+ $ WR1, WR2, XNORM
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
+ $ X( LDX, 2 )
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION DLAMCH, DLANGE
+ EXTERNAL DLAMCH, DLANGE
+* ..
+* .. External Subroutines ..
+ EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
+ $ DROT
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX
+* ..
+* .. Executable Statements ..
+*
+ INFO = 0
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
+ $ RETURN
+ IF( J1+N1.GT.N )
+ $ RETURN
+*
+ J2 = J1 + 1
+ J3 = J1 + 2
+ J4 = J1 + 3
+*
+ IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
+*
+* Swap two 1-by-1 blocks.
+*
+ T11 = T( J1, J1 )
+ T22 = T( J2, J2 )
+*
+* Determine the transformation to perform the interchange.
+*
+ CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
+*
+* Apply transformation to the matrix T.
+*
+ IF( J3.LE.N )
+ $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
+ $ SN )
+ CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
+*
+ T( J1, J1 ) = T22
+ T( J2, J2 ) = T11
+*
+ IF( WANTQ ) THEN
+*
+* Accumulate transformation in the matrix Q.
+*
+ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
+ END IF
+*
+ ELSE
+*
+* Swapping involves at least one 2-by-2 block.
+*
+* Copy the diagonal block of order N1+N2 to the local array D
+* and compute its norm.
+*
+ ND = N1 + N2
+ CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
+ DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
+*
+* Compute machine-dependent threshold for test for accepting
+* swap.
+*
+ EPS = DLAMCH( 'P' )
+ SMLNUM = DLAMCH( 'S' ) / EPS
+ THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
+*
+* Solve T11*X - X*T22 = scale*T12 for X.
+*
+ CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
+ $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
+ $ LDX, XNORM, IERR )
+*
+* Swap the adjacent diagonal blocks.
+*
+ K = N1 + N1 + N2 - 3
+ GO TO ( 10, 20, 30 )K
+*
+ 10 CONTINUE
+*
+* N1 = 1, N2 = 2: generate elementary reflector H so that:
+*
+* ( scale, X11, X12 ) H = ( 0, 0, * )
+*
+ U( 1 ) = SCALE
+ U( 2 ) = X( 1, 1 )
+ U( 3 ) = X( 1, 2 )
+ CALL DLARFG( 3, U( 3 ), U, 1, TAU )
+ U( 3 ) = ONE
+ T11 = T( J1, J1 )
+*
+* Perform swap provisionally on diagonal block in D.
+*
+ CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
+ CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
+*
+* Test whether to reject swap.
+*
+ IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
+ $ 3 )-T11 ) ).GT.THRESH )GO TO 50
+*
+* Accept swap: apply transformation to the entire matrix T.
+*
+ CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
+ CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
+*
+ T( J3, J1 ) = ZERO
+ T( J3, J2 ) = ZERO
+ T( J3, J3 ) = T11
+*
+ IF( WANTQ ) THEN
+*
+* Accumulate transformation in the matrix Q.
+*
+ CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
+ END IF
+ GO TO 40
+*
+ 20 CONTINUE
+*
+* N1 = 2, N2 = 1: generate elementary reflector H so that:
+*
+* H ( -X11 ) = ( * )
+* ( -X21 ) = ( 0 )
+* ( scale ) = ( 0 )
+*
+ U( 1 ) = -X( 1, 1 )
+ U( 2 ) = -X( 2, 1 )
+ U( 3 ) = SCALE
+ CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
+ U( 1 ) = ONE
+ T33 = T( J3, J3 )
+*
+* Perform swap provisionally on diagonal block in D.
+*
+ CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
+ CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
+*
+* Test whether to reject swap.
+*
+ IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
+ $ 1 )-T33 ) ).GT.THRESH )GO TO 50
+*
+* Accept swap: apply transformation to the entire matrix T.
+*
+ CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
+ CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
+*
+ T( J1, J1 ) = T33
+ T( J2, J1 ) = ZERO
+ T( J3, J1 ) = ZERO
+*
+ IF( WANTQ ) THEN
+*
+* Accumulate transformation in the matrix Q.
+*
+ CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
+ END IF
+ GO TO 40
+*
+ 30 CONTINUE
+*
+* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
+* that:
+*
+* H(2) H(1) ( -X11 -X12 ) = ( * * )
+* ( -X21 -X22 ) ( 0 * )
+* ( scale 0 ) ( 0 0 )
+* ( 0 scale ) ( 0 0 )
+*
+ U1( 1 ) = -X( 1, 1 )
+ U1( 2 ) = -X( 2, 1 )
+ U1( 3 ) = SCALE
+ CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
+ U1( 1 ) = ONE
+*
+ TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
+ U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
+ U2( 2 ) = -TEMP*U1( 3 )
+ U2( 3 ) = SCALE
+ CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
+ U2( 1 ) = ONE
+*
+* Perform swap provisionally on diagonal block in D.
+*
+ CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
+ CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
+ CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
+ CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
+*
+* Test whether to reject swap.
+*
+ IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
+ $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
+*
+* Accept swap: apply transformation to the entire matrix T.
+*
+ CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
+ CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
+ CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
+ CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
+*
+ T( J3, J1 ) = ZERO
+ T( J3, J2 ) = ZERO
+ T( J4, J1 ) = ZERO
+ T( J4, J2 ) = ZERO
+*
+ IF( WANTQ ) THEN
+*
+* Accumulate transformation in the matrix Q.
+*
+ CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
+ CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
+ END IF
+*
+ 40 CONTINUE
+*
+ IF( N2.EQ.2 ) THEN
+*
+* Standardize new 2-by-2 block T11
+*
+ CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
+ $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
+ CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
+ $ CS, SN )
+ CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
+ IF( WANTQ )
+ $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
+ END IF
+*
+ IF( N1.EQ.2 ) THEN
+*
+* Standardize new 2-by-2 block T22
+*
+ J3 = J1 + N2
+ J4 = J3 + 1
+ CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
+ $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
+ IF( J3+2.LE.N )
+ $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
+ $ LDT, CS, SN )
+ CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
+ IF( WANTQ )
+ $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
+ END IF
+*
+ END IF
+ RETURN
+*
+* Exit with INFO = 1 if swap was rejected.
+*
+ 50 CONTINUE
+ INFO = 1
+ RETURN
+*
+* End of DLAEXC
+*
+ END
generated by cgit (git 2.34.1) at 2024-12-20 10:59:20 +0000