Move LAPACK trunk into position.

1 files changed, 345 insertions, 0 deletions

diff --git a/SRC/dtrexc.f b/SRC/dtrexc.f
new file mode 100644
index 00000000..db9be753
--- /dev/null
+++ b/SRC/dtrexc.f
@@ -0,0 +1,345 @@
+      SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
+     $                   INFO )
+*
+*  -- LAPACK routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          COMPQ
+      INTEGER            IFST, ILST, INFO, LDQ, LDT, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DTREXC reorders the real Schur factorization of a real matrix
+*  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
+*  moved to row ILST.
+*
+*  The real Schur form T is reordered by an orthogonal similarity
+*  transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
+*  is updated by postmultiplying it with Z.
+*
+*  T must be in Schur canonical form (as returned by DHSEQR), 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.
+*
+*  Arguments
+*  =========
+*
+*  COMPQ   (input) CHARACTER*1
+*          = 'V':  update the matrix Q of Schur vectors;
+*          = 'N':  do not update Q.
+*
+*  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
+*          Schur canonical form.
+*          On exit, the reordered upper quasi-triangular matrix, 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 COMPQ = 'V', the matrix Q of Schur vectors.
+*          On exit, if COMPQ = 'V', Q has been postmultiplied by the
+*          orthogonal transformation matrix Z which reorders T.
+*          If COMPQ = 'N', Q is not referenced.
+*
+*  LDQ     (input) INTEGER
+*          The leading dimension of the array Q.  LDQ >= max(1,N).
+*
+*  IFST    (input/output) INTEGER
+*  ILST    (input/output) INTEGER
+*          Specify the reordering of the diagonal blocks of T.
+*          The block with row index IFST is moved to row ILST, by a
+*          sequence of transpositions between adjacent blocks.
+*          On exit, if IFST pointed on entry to the second row of a
+*          2-by-2 block, it is changed to point to the first row; ILST
+*          always points to the first row of the block in its final
+*          position (which may differ from its input value by +1 or -1).
+*          1 <= IFST <= N; 1 <= ILST <= N.
+*
+*  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*          = 1:  two adjacent blocks were too close to swap (the problem
+*                is very ill-conditioned); T may have been partially
+*                reordered, and ILST points to the first row of the
+*                current position of the block being moved.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            WANTQ
+      INTEGER            HERE, NBF, NBL, NBNEXT
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLAEXC, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Decode and test the input arguments.
+*
+      INFO = 0
+      WANTQ = LSAME( COMPQ, 'V' )
+      IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
+         INFO = -6
+      ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
+         INFO = -7
+      ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DTREXC', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.LE.1 )
+     $   RETURN
+*
+*     Determine the first row of specified block
+*     and find out it is 1 by 1 or 2 by 2.
+*
+      IF( IFST.GT.1 ) THEN
+         IF( T( IFST, IFST-1 ).NE.ZERO )
+     $      IFST = IFST - 1
+      END IF
+      NBF = 1
+      IF( IFST.LT.N ) THEN
+         IF( T( IFST+1, IFST ).NE.ZERO )
+     $      NBF = 2
+      END IF
+*
+*     Determine the first row of the final block
+*     and find out it is 1 by 1 or 2 by 2.
+*
+      IF( ILST.GT.1 ) THEN
+         IF( T( ILST, ILST-1 ).NE.ZERO )
+     $      ILST = ILST - 1
+      END IF
+      NBL = 1
+      IF( ILST.LT.N ) THEN
+         IF( T( ILST+1, ILST ).NE.ZERO )
+     $      NBL = 2
+      END IF
+*
+      IF( IFST.EQ.ILST )
+     $   RETURN
+*
+      IF( IFST.LT.ILST ) THEN
+*
+*        Update ILST
+*
+         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
+     $      ILST = ILST - 1
+         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
+     $      ILST = ILST + 1
+*
+         HERE = IFST
+*
+   10    CONTINUE
+*
+*        Swap block with next one below
+*
+         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
+*
+*           Current block either 1 by 1 or 2 by 2
+*
+            NBNEXT = 1
+            IF( HERE+NBF+1.LE.N ) THEN
+               IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
+     $            NBNEXT = 2
+            END IF
+            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
+     $                   WORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               ILST = HERE
+               RETURN
+            END IF
+            HERE = HERE + NBNEXT
+*
+*           Test if 2 by 2 block breaks into two 1 by 1 blocks
+*
+            IF( NBF.EQ.2 ) THEN
+               IF( T( HERE+1, HERE ).EQ.ZERO )
+     $            NBF = 3
+            END IF
+*
+         ELSE
+*
+*           Current block consists of two 1 by 1 blocks each of which
+*           must be swapped individually
+*
+            NBNEXT = 1
+            IF( HERE+3.LE.N ) THEN
+               IF( T( HERE+3, HERE+2 ).NE.ZERO )
+     $            NBNEXT = 2
+            END IF
+            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
+     $                   WORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               ILST = HERE
+               RETURN
+            END IF
+            IF( NBNEXT.EQ.1 ) THEN
+*
+*              Swap two 1 by 1 blocks, no problems possible
+*
+               CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
+     $                      WORK, INFO )
+               HERE = HERE + 1
+            ELSE
+*
+*              Recompute NBNEXT in case 2 by 2 split
+*
+               IF( T( HERE+2, HERE+1 ).EQ.ZERO )
+     $            NBNEXT = 1
+               IF( NBNEXT.EQ.2 ) THEN
+*
+*                 2 by 2 Block did not split
+*
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
+     $                         NBNEXT, WORK, INFO )
+                  IF( INFO.NE.0 ) THEN
+                     ILST = HERE
+                     RETURN
+                  END IF
+                  HERE = HERE + 2
+               ELSE
+*
+*                 2 by 2 Block did split
+*
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
+     $                         WORK, INFO )
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
+     $                         WORK, INFO )
+                  HERE = HERE + 2
+               END IF
+            END IF
+         END IF
+         IF( HERE.LT.ILST )
+     $      GO TO 10
+*
+      ELSE
+*
+         HERE = IFST
+   20    CONTINUE
+*
+*        Swap block with next one above
+*
+         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
+*
+*           Current block either 1 by 1 or 2 by 2
+*
+            NBNEXT = 1
+            IF( HERE.GE.3 ) THEN
+               IF( T( HERE-1, HERE-2 ).NE.ZERO )
+     $            NBNEXT = 2
+            END IF
+            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
+     $                   NBF, WORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               ILST = HERE
+               RETURN
+            END IF
+            HERE = HERE - NBNEXT
+*
+*           Test if 2 by 2 block breaks into two 1 by 1 blocks
+*
+            IF( NBF.EQ.2 ) THEN
+               IF( T( HERE+1, HERE ).EQ.ZERO )
+     $            NBF = 3
+            END IF
+*
+         ELSE
+*
+*           Current block consists of two 1 by 1 blocks each of which
+*           must be swapped individually
+*
+            NBNEXT = 1
+            IF( HERE.GE.3 ) THEN
+               IF( T( HERE-1, HERE-2 ).NE.ZERO )
+     $            NBNEXT = 2
+            END IF
+            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
+     $                   1, WORK, INFO )
+            IF( INFO.NE.0 ) THEN
+               ILST = HERE
+               RETURN
+            END IF
+            IF( NBNEXT.EQ.1 ) THEN
+*
+*              Swap two 1 by 1 blocks, no problems possible
+*
+               CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
+     $                      WORK, INFO )
+               HERE = HERE - 1
+            ELSE
+*
+*              Recompute NBNEXT in case 2 by 2 split
+*
+               IF( T( HERE, HERE-1 ).EQ.ZERO )
+     $            NBNEXT = 1
+               IF( NBNEXT.EQ.2 ) THEN
+*
+*                 2 by 2 Block did not split
+*
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
+     $                         WORK, INFO )
+                  IF( INFO.NE.0 ) THEN
+                     ILST = HERE
+                     RETURN
+                  END IF
+                  HERE = HERE - 2
+               ELSE
+*
+*                 2 by 2 Block did split
+*
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
+     $                         WORK, INFO )
+                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
+     $                         WORK, INFO )
+                  HERE = HERE - 2
+               END IF
+            END IF
+         END IF
+         IF( HERE.GT.ILST )
+     $      GO TO 20
+      END IF
+      ILST = HERE
+*
+      RETURN
+*
+*     End of DTREXC
+*
+      END