Mercurial > octave

diff libcruft/lapack/slaexc.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/lapack/slaexc.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,353 @@
+      SUBROUTINE SLAEXC( 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 ..
+      REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  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 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) 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.
+*
+*  LDT     (input)  INTEGER
+*          The leading dimension of the array T. LDT >= max(1,N).
+*
+*  Q       (input/output) 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.
+*
+*  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) REAL 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 ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      REAL               TEN
+      PARAMETER          ( TEN = 1.0E+1 )
+      INTEGER            LDD, LDX
+      PARAMETER          ( LDD = 4, LDX = 2 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            IERR, J2, J3, J4, K, ND
+      REAL               CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
+     $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
+     $                   WR1, WR2, XNORM
+*     ..
+*     .. Local Arrays ..
+      REAL               D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
+     $                   X( LDX, 2 )
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           SLAMCH, SLANGE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLACPY, SLANV2, SLARFG, SLARFX, SLARTG, SLASY2,
+     $                   SROT
+*     ..
+*     .. 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 SLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
+*
+*        Apply transformation to the matrix T.
+*
+         IF( J3.LE.N )
+     $      CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
+     $                 SN )
+         CALL SROT( 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 SROT( 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 SLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
+         DNORM = SLANGE( 'Max', ND, ND, D, LDD, WORK )
+*
+*        Compute machine-dependent threshold for test for accepting
+*        swap.
+*
+         EPS = SLAMCH( 'P' )
+         SMLNUM = SLAMCH( 'S' ) / EPS
+         THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
+*
+*        Solve T11*X - X*T22 = scale*T12 for X.
+*
+         CALL SLASY2( .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 SLARFG( 3, U( 3 ), U, 1, TAU )
+         U( 3 ) = ONE
+         T11 = T( J1, J1 )
+*
+*        Perform swap provisionally on diagonal block in D.
+*
+         CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
+         CALL SLARFX( '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 SLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
+         CALL SLARFX( '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 SLARFX( '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 SLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
+         U( 1 ) = ONE
+         T33 = T( J3, J3 )
+*
+*        Perform swap provisionally on diagonal block in D.
+*
+         CALL SLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
+         CALL SLARFX( '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 SLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
+         CALL SLARFX( '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 SLARFX( '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 SLARFG( 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 SLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
+         U2( 1 ) = ONE
+*
+*        Perform swap provisionally on diagonal block in D.
+*
+         CALL SLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
+         CALL SLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
+         CALL SLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
+         CALL SLARFX( '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 SLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
+         CALL SLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
+         CALL SLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
+         CALL SLARFX( '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 SLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
+            CALL SLARFX( '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 SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
+     $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
+            CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
+     $                 CS, SN )
+            CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
+            IF( WANTQ )
+     $         CALL SROT( 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 SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
+     $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
+            IF( J3+2.LE.N )
+     $         CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
+     $                    LDT, CS, SN )
+            CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
+            IF( WANTQ )
+     $         CALL SROT( 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 INFO = 1
+      RETURN
+*
+*     End of SLAEXC
+*
+      END