Move LAPACK trunk into position.

1 files changed, 322 insertions, 0 deletions

diff --git a/SRC/dgebal.f b/SRC/dgebal.f
new file mode 100644
index 00000000..1796577b
--- /dev/null
+++ b/SRC/dgebal.f
@@ -0,0 +1,322 @@
+      SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
+*
+*  -- LAPACK routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOB
+      INTEGER            IHI, ILO, INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), SCALE( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DGEBAL balances a general real matrix A.  This involves, first,
+*  permuting A by a similarity transformation to isolate eigenvalues
+*  in the first 1 to ILO-1 and last IHI+1 to N elements on the
+*  diagonal; and second, applying a diagonal similarity transformation
+*  to rows and columns ILO to IHI to make the rows and columns as
+*  close in norm as possible.  Both steps are optional.
+*
+*  Balancing may reduce the 1-norm of the matrix, and improve the
+*  accuracy of the computed eigenvalues and/or eigenvectors.
+*
+*  Arguments
+*  =========
+*
+*  JOB     (input) CHARACTER*1
+*          Specifies the operations to be performed on A:
+*          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
+*                  for i = 1,...,N;
+*          = 'P':  permute only;
+*          = 'S':  scale only;
+*          = 'B':  both permute and scale.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+*          On entry, the input matrix A.
+*          On exit,  A is overwritten by the balanced matrix.
+*          If JOB = 'N', A is not referenced.
+*          See Further Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  ILO     (output) INTEGER
+*  IHI     (output) INTEGER
+*          ILO and IHI are set to integers such that on exit
+*          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
+*          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
+*
+*  SCALE   (output) DOUBLE PRECISION array, dimension (N)
+*          Details of the permutations and scaling factors applied to
+*          A.  If P(j) is the index of the row and column interchanged
+*          with row and column j and D(j) is the scaling factor
+*          applied to row and column j, then
+*          SCALE(j) = P(j)    for j = 1,...,ILO-1
+*                   = D(j)    for j = ILO,...,IHI
+*                   = P(j)    for j = IHI+1,...,N.
+*          The order in which the interchanges are made is N to IHI+1,
+*          then 1 to ILO-1.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit.
+*          < 0:  if INFO = -i, the i-th argument had an illegal value.
+*
+*  Further Details
+*  ===============
+*
+*  The permutations consist of row and column interchanges which put
+*  the matrix in the form
+*
+*             ( T1   X   Y  )
+*     P A P = (  0   B   Z  )
+*             (  0   0   T2 )
+*
+*  where T1 and T2 are upper triangular matrices whose eigenvalues lie
+*  along the diagonal.  The column indices ILO and IHI mark the starting
+*  and ending columns of the submatrix B. Balancing consists of applying
+*  a diagonal similarity transformation inv(D) * B * D to make the
+*  1-norms of each row of B and its corresponding column nearly equal.
+*  The output matrix is
+*
+*     ( T1     X*D          Y    )
+*     (  0  inv(D)*B*D  inv(D)*Z ).
+*     (  0      0           T2   )
+*
+*  Information about the permutations P and the diagonal matrix D is
+*  returned in the vector SCALE.
+*
+*  This subroutine is based on the EISPACK routine BALANC.
+*
+*  Modified by Tzu-Yi Chen, Computer Science Division, University of
+*    California at Berkeley, USA
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      DOUBLE PRECISION   SCLFAC
+      PARAMETER          ( SCLFAC = 2.0D+0 )
+      DOUBLE PRECISION   FACTOR
+      PARAMETER          ( FACTOR = 0.95D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            NOCONV
+      INTEGER            I, ICA, IEXC, IRA, J, K, L, M
+      DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
+     $                   SFMIN2
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            IDAMAX
+      DOUBLE PRECISION   DLAMCH
+      EXTERNAL           LSAME, IDAMAX, DLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DSCAL, DSWAP, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters
+*
+      INFO = 0
+      IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
+     $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DGEBAL', -INFO )
+         RETURN
+      END IF
+*
+      K = 1
+      L = N
+*
+      IF( N.EQ.0 )
+     $   GO TO 210
+*
+      IF( LSAME( JOB, 'N' ) ) THEN
+         DO 10 I = 1, N
+            SCALE( I ) = ONE
+   10    CONTINUE
+         GO TO 210
+      END IF
+*
+      IF( LSAME( JOB, 'S' ) )
+     $   GO TO 120
+*
+*     Permutation to isolate eigenvalues if possible
+*
+      GO TO 50
+*
+*     Row and column exchange.
+*
+   20 CONTINUE
+      SCALE( M ) = J
+      IF( J.EQ.M )
+     $   GO TO 30
+*
+      CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
+      CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
+*
+   30 CONTINUE
+      GO TO ( 40, 80 )IEXC
+*
+*     Search for rows isolating an eigenvalue and push them down.
+*
+   40 CONTINUE
+      IF( L.EQ.1 )
+     $   GO TO 210
+      L = L - 1
+*
+   50 CONTINUE
+      DO 70 J = L, 1, -1
+*
+         DO 60 I = 1, L
+            IF( I.EQ.J )
+     $         GO TO 60
+            IF( A( J, I ).NE.ZERO )
+     $         GO TO 70
+   60    CONTINUE
+*
+         M = L
+         IEXC = 1
+         GO TO 20
+   70 CONTINUE
+*
+      GO TO 90
+*
+*     Search for columns isolating an eigenvalue and push them left.
+*
+   80 CONTINUE
+      K = K + 1
+*
+   90 CONTINUE
+      DO 110 J = K, L
+*
+         DO 100 I = K, L
+            IF( I.EQ.J )
+     $         GO TO 100
+            IF( A( I, J ).NE.ZERO )
+     $         GO TO 110
+  100    CONTINUE
+*
+         M = K
+         IEXC = 2
+         GO TO 20
+  110 CONTINUE
+*
+  120 CONTINUE
+      DO 130 I = K, L
+         SCALE( I ) = ONE
+  130 CONTINUE
+*
+      IF( LSAME( JOB, 'P' ) )
+     $   GO TO 210
+*
+*     Balance the submatrix in rows K to L.
+*
+*     Iterative loop for norm reduction
+*
+      SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
+      SFMAX1 = ONE / SFMIN1
+      SFMIN2 = SFMIN1*SCLFAC
+      SFMAX2 = ONE / SFMIN2
+  140 CONTINUE
+      NOCONV = .FALSE.
+*
+      DO 200 I = K, L
+         C = ZERO
+         R = ZERO
+*
+         DO 150 J = K, L
+            IF( J.EQ.I )
+     $         GO TO 150
+            C = C + ABS( A( J, I ) )
+            R = R + ABS( A( I, J ) )
+  150    CONTINUE
+         ICA = IDAMAX( L, A( 1, I ), 1 )
+         CA = ABS( A( ICA, I ) )
+         IRA = IDAMAX( N-K+1, A( I, K ), LDA )
+         RA = ABS( A( I, IRA+K-1 ) )
+*
+*        Guard against zero C or R due to underflow.
+*
+         IF( C.EQ.ZERO .OR. R.EQ.ZERO )
+     $      GO TO 200
+         G = R / SCLFAC
+         F = ONE
+         S = C + R
+  160    CONTINUE
+         IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
+     $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
+         F = F*SCLFAC
+         C = C*SCLFAC
+         CA = CA*SCLFAC
+         R = R / SCLFAC
+         G = G / SCLFAC
+         RA = RA / SCLFAC
+         GO TO 160
+*
+  170    CONTINUE
+         G = C / SCLFAC
+  180    CONTINUE
+         IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
+     $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
+         F = F / SCLFAC
+         C = C / SCLFAC
+         G = G / SCLFAC
+         CA = CA / SCLFAC
+         R = R*SCLFAC
+         RA = RA*SCLFAC
+         GO TO 180
+*
+*        Now balance.
+*
+  190    CONTINUE
+         IF( ( C+R ).GE.FACTOR*S )
+     $      GO TO 200
+         IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
+            IF( F*SCALE( I ).LE.SFMIN1 )
+     $         GO TO 200
+         END IF
+         IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
+            IF( SCALE( I ).GE.SFMAX1 / F )
+     $         GO TO 200
+         END IF
+         G = ONE / F
+         SCALE( I ) = SCALE( I )*F
+         NOCONV = .TRUE.
+*
+         CALL DSCAL( N-K+1, G, A( I, K ), LDA )
+         CALL DSCAL( L, F, A( 1, I ), 1 )
+*
+  200 CONTINUE
+*
+      IF( NOCONV )
+     $   GO TO 140
+*
+  210 CONTINUE
+      ILO = K
+      IHI = L
+*
+      RETURN
+*
+*     End of DGEBAL
+*
+      END