Refactored GEBAL family

Author: eprovst

SRC/cgebal.f

@@ -85,6 +85,7 @@
 *> \verbatim
 *>          ILO is INTEGER
 *> \endverbatim
+*>
 *> \param[out] IHI
 *> \verbatim
 *>          IHI is INTEGER
@@ -154,6 +155,9 @@
 *>
 *>  Modified by Tzu-Yi Chen, Computer Science Division, University of
 *>    California at Berkeley, USA
+*>
+*>  Refactored by Evert Provoost, Department of Computer Science,
+*>    KU Leuven, Belgium
 *> \endverbatim
 *>
 *  =====================================================================
@@ -183,8 +187,8 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
       PARAMETER          ( FACTOR = 0.95E+0 )
 *     ..
 *     .. Local Scalars ..
-      LOGICAL            NOCONV
-      INTEGER            I, ICA, IEXC, IRA, J, K, L, M
+      LOGICAL            NOCONV, CANSWAP
+      INTEGER            I, ICA, IRA, J, K, L
       REAL               C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
      $                   SFMIN2
 *     ..
@@ -195,10 +199,10 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
       EXTERNAL           SISNAN, LSAME, ICAMAX, SLAMCH, SCNRM2
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CSSCAL, CSWAP, XERBLA
+      EXTERNAL           XERBLA, CSSCAL, CSWAP
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
+      INTRINSIC          ABS, REAL, AIMAG, MAX, MIN
 *
 *     Test the input parameters
 *
@@ -216,176 +220,196 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
          RETURN
       END IF
 *
-      K = 1
-      L = N
+*     Quick returns.
 *
-      IF( N.EQ.0 )
-     $   GO TO 210
+      IF( N.EQ.0 ) THEN
+         ILO = 1
+         IHI = 0
+         RETURN
+      END IF
 *
       IF( LSAME( JOB, 'N' ) ) THEN
-         DO 10 I = 1, N
+         DO I = 1, N
             SCALE( I ) = ONE
-   10    CONTINUE
-         GO TO 210
+         END DO
+         ILO = 1
+         IHI = N
+         RETURN
       END IF
 *
-      IF( LSAME( JOB, 'S' ) )
-     $   GO TO 120
+*     Permutation to isolate eigenvalues if possible.
+*
+      K = 1
+      L = N
+*
+      IF( .NOT.LSAME( JOB, 'S' ) ) THEN
+*
+*        Row and column exchange.
+*
+         NOCONV = .TRUE.
+         DO WHILE( NOCONV )
+*
+*           Search for rows isolating an eigenvalue and push them down.
+*
+            NOCONV = .FALSE.
+            DO I = L, 1, -1
+               CANSWAP = .TRUE.
+               DO J = 1, L
+                  IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
+     $                AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
+                     CANSWAP = .FALSE.
+                     EXIT
+                  END IF
+               END DO
+*
+               IF( CANSWAP ) THEN
+                  SCALE( L ) = I
+                  IF( I.NE.L ) THEN
+                     CALL CSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
+                     CALL CSWAP( N-K+1, A( I, K ), LDA, A( L, K ), LDA )
+                  END IF
+                  NOCONV = .TRUE.
+*
+                  IF( L.EQ.1 ) THEN
+                     ILO = 1
+                     IHI = 1
+                     RETURN
+                  END IF
+*
+                  L = L - 1
+               END IF
+            END DO
+*
+         END DO
+
+         NOCONV = .TRUE.
+         DO WHILE( NOCONV )
+*
+*           Search for columns isolating an eigenvalue and push them left.
+*
+            NOCONV = .FALSE.
+            DO J = K, L
+               CANSWAP = .TRUE.
+               DO I = K, L
+                  IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
+     $                AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
+                     CANSWAP = .FALSE.
+                     EXIT
+                  END IF
+               END DO
+*
+               IF( CANSWAP ) THEN
+                  SCALE( K ) = J
+                  IF( J.NE.K ) THEN
+                     CALL CSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
+                     CALL CSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
+                  END IF
+                  NOCONV = .TRUE.
+*
+                  K = K + 1
+               END IF
+            END DO
+*
+         END DO
+*
+      END IF
 *
-*     Permutation to isolate eigenvalues if possible
+*     Initialize SCALE for non-permuted submatrix.
 *
-      GO TO 50
+      DO I = K, L
+         SCALE( I ) = ONE
+      END DO
 *
-*     Row and column exchange.
+*     If we only had to permute, we are done.
 *
-   20 CONTINUE
-      SCALE( M ) = J
-      IF( J.EQ.M )
-     $   GO TO 30
+      IF( LSAME( JOB, 'P' ) ) THEN
+         ILO = K
+         IHI = L
+         RETURN
+      END IF
 *
-      CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
-      CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
+*     Balance the submatrix in rows K to L.
 *
-   30 CONTINUE
-      GO TO ( 40, 80 )IEXC
+*     Iterative loop for norm reduction.
 *
-*     Search for rows isolating an eigenvalue and push them down.
+      SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
+      SFMAX1 = ONE / SFMIN1
+      SFMIN2 = SFMIN1*SCLFAC
+      SFMAX2 = ONE / SFMIN2
 *
-   40 CONTINUE
-      IF( L.EQ.1 )
-     $   GO TO 210
-      L = L - 1
+      NOCONV = .TRUE.
+      DO WHILE( NOCONV )
+         NOCONV = .FALSE.
 *
-   50 CONTINUE
-      DO 70 J = L, 1, -1
+         DO I = K, L
 *
-         DO 60 I = 1, L
-            IF( I.EQ.J )
-     $         GO TO 60
-            IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE.
-     $          ZERO )GO TO 70
-   60    CONTINUE
+            C = SCNRM2( L-K+1, A( K, I ), 1 )
+            R = SCNRM2( L-K+1, A( I, K ), LDA )
+            ICA = ICAMAX( L, A( 1, I ), 1 )
+            CA = ABS( A( ICA, I ) )
+            IRA = ICAMAX( N-K+1, A( I, K ), LDA )
+            RA = ABS( A( I, IRA+K-1 ) )
 *
-         M = L
-         IEXC = 1
-         GO TO 20
-   70 CONTINUE
+*           Guard against zero C or R due to underflow.
 *
-      GO TO 90
+            IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
 *
-*     Search for columns isolating an eigenvalue and push them left.
+            G = R / SCLFAC
+            F = ONE
+            S = C + R
 *
-   80 CONTINUE
-      K = K + 1
+            DO WHILE( C.LT.G .AND. MAX( F, C, CA ).LT.SFMAX2 .AND.
+     $                MIN( R, G, RA ).GT.SFMIN2 )
 *
-   90 CONTINUE
-      DO 110 J = K, L
+               IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
 *
-         DO 100 I = K, L
-            IF( I.EQ.J )
-     $         GO TO 100
-            IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE.
-     $          ZERO )GO TO 110
-  100    CONTINUE
+*                 Exit if NaN to avoid infinite loop
 *
-         M = K
-         IEXC = 2
-         GO TO 20
-  110 CONTINUE
+                  INFO = -3
+                  CALL XERBLA( 'CGEBAL', -INFO )
+                  RETURN
+               END IF
 *
-  120 CONTINUE
-      DO 130 I = K, L
-         SCALE( I ) = ONE
-  130 CONTINUE
+               F = F*SCLFAC
+               C = C*SCLFAC
+               CA = CA*SCLFAC
+               R = R / SCLFAC
+               G = G / SCLFAC
+               RA = RA / SCLFAC
+            END DO
 *
-      IF( LSAME( JOB, 'P' ) )
-     $   GO TO 210
+            G = C / SCLFAC
 *
-*     Balance the submatrix in rows K to L.
+            DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
+     $                MIN( F, C, G, CA ).GT.SFMIN2 )
+               F = F / SCLFAC
+               C = C / SCLFAC
+               G = G / SCLFAC
+               CA = CA / SCLFAC
+               R = R*SCLFAC
+               RA = RA*SCLFAC
+            END DO
 *
-*     Iterative loop for norm reduction
+*           Now balance.
 *
-      SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
-      SFMAX1 = ONE / SFMIN1
-      SFMIN2 = SFMIN1*SCLFAC
-      SFMAX2 = ONE / SFMIN2
-  140 CONTINUE
-      NOCONV = .FALSE.
-*
-      DO 200 I = K, L
-*
-         C = SCNRM2( L-K+1, A( K, I ), 1 )
-         R = SCNRM2( L-K+1, A( I , K ), LDA )
-         ICA = ICAMAX( L, A( 1, I ), 1 )
-         CA = ABS( A( ICA, I ) )
-         IRA = ICAMAX( 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
-            IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
-*
-*           Exit if NaN to avoid infinite loop
-*
-            INFO = -3
-            CALL XERBLA( 'CGEBAL', -INFO )
-            RETURN
-         END IF
-         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.
+            IF( ( C+R ).GE.FACTOR*S ) CYCLE
+            IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
+               IF( F*SCALE( I ).LE.SFMIN1 ) CYCLE
+            END IF
+            IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
+               IF( SCALE( I ).GE.SFMAX1 / F ) CYCLE
+            END IF
+            G = ONE / F
+            SCALE( I ) = SCALE( I )*F
+            NOCONV = .TRUE.
 *
-         CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
-         CALL CSSCAL( L, F, A( 1, I ), 1 )
+            CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
+            CALL CSSCAL( L, F, A( 1, I ), 1 )
 *
-  200 CONTINUE
+         END DO
 *
-      IF( NOCONV )
-     $   GO TO 140
+      END DO
 *
-  210 CONTINUE
       ILO = K
       IHI = L
 *