sgebal.f

Go to the documentation of this file.

*> \brief \b SGEBAL
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download SGEBAL + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgebal.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgebal.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgebal.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE SGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          JOB
*       INTEGER            IHI, ILO, INFO, LDA, N
*       ..
*       .. Array Arguments ..
*       REAL               A( LDA, * ), SCALE( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SGEBAL 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOB
*> \verbatim
*>          JOB is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*> \param[out] IHI
*> \verbatim
*>          IHI is 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.
*> \endverbatim
*>
*> \param[out] SCALE
*> \verbatim
*>          SCALE is REAL 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.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2015
*
*> \ingroup realGEcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE sgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
*
*  -- LAPACK computational routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      CHARACTER          JOB
      INTEGER            IHI, ILO, INFO, LDA, N
*     ..
*     .. Array Arguments ..
      REAL               A( lda, * ), SCALE( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter                ( zero = 0.0e+0, one = 1.0e+0 )
      REAL               SCLFAC
      parameter                ( sclfac = 2.0e+0 )
      REAL               FACTOR
      parameter                ( factor = 0.95e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOCONV
      INTEGER            I, ICA, IEXC, IRA, J, K, L, M
      REAL               C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
     $                   sfmin2
*     ..
*     .. External Functions ..
      LOGICAL            SISNAN, LSAME
      INTEGER            ISAMAX
      REAL               SLAMCH, SNRM2
      EXTERNAL           sisnan, lsame, isamax, slamch, snrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           sscal, sswap, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*
*     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( 'SGEBAL', -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 sswap( l, a( 1, j ), 1, a( 1, m ), 1 )
      CALL sswap( 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 = slamch( 'S' ) / slamch( 'P' )
      sfmax1 = one / sfmin1
      sfmin2 = sfmin1*sclfac
      sfmax2 = one / sfmin2
  140 CONTINUE
      noconv = .false.
*
      DO 200 i = k, l
*
         c = snrm2( l-k+1, a( k, i ), 1 )
         r = snrm2( l-k+1, a( i, k ), lda )
         ica = isamax( l, a( 1, i ), 1 )
         ca = abs( a( ica, i ) )
         ira = isamax( 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
            IF( sisnan( c+f+ca+r+g+ra ) ) THEN
*
*           Exit if NaN to avoid infinite loop
*
            info = -3
            CALL xerbla( 'SGEBAL', -info )
            RETURN
         END IF
         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 sscal( n-k+1, g, a( i, k ), lda )
         CALL sscal( l, f, a( 1, i ), 1 )
*
  200 CONTINUE
*
      IF( noconv )
     $   GO TO 140
*
  210 CONTINUE
      ilo = k
      ihi = l
*
      RETURN
*
*     End of SGEBAL
*
      END