LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgebal.f
Go to the documentation of this file.
00001 *> \brief \b SGEBAL
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGEBAL + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgebal.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgebal.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgebal.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          JOB
00025 *       INTEGER            IHI, ILO, INFO, LDA, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               A( LDA, * ), SCALE( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> SGEBAL balances a general real matrix A.  This involves, first,
00038 *> permuting A by a similarity transformation to isolate eigenvalues
00039 *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
00040 *> diagonal; and second, applying a diagonal similarity transformation
00041 *> to rows and columns ILO to IHI to make the rows and columns as
00042 *> close in norm as possible.  Both steps are optional.
00043 *>
00044 *> Balancing may reduce the 1-norm of the matrix, and improve the
00045 *> accuracy of the computed eigenvalues and/or eigenvectors.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] JOB
00052 *> \verbatim
00053 *>          JOB is CHARACTER*1
00054 *>          Specifies the operations to be performed on A:
00055 *>          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
00056 *>                  for i = 1,...,N;
00057 *>          = 'P':  permute only;
00058 *>          = 'S':  scale only;
00059 *>          = 'B':  both permute and scale.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.  N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in,out] A
00069 *> \verbatim
00070 *>          A is REAL array, dimension (LDA,N)
00071 *>          On entry, the input matrix A.
00072 *>          On exit,  A is overwritten by the balanced matrix.
00073 *>          If JOB = 'N', A is not referenced.
00074 *>          See Further Details.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDA
00078 *> \verbatim
00079 *>          LDA is INTEGER
00080 *>          The leading dimension of the array A.  LDA >= max(1,N).
00081 *> \endverbatim
00082 *>
00083 *> \param[out] ILO
00084 *> \verbatim
00085 *>          ILO is INTEGER
00086 *> \endverbatim
00087 *> \param[out] IHI
00088 *> \verbatim
00089 *>          IHI is INTEGER
00090 *>          ILO and IHI are set to integers such that on exit
00091 *>          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
00092 *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
00093 *> \endverbatim
00094 *>
00095 *> \param[out] SCALE
00096 *> \verbatim
00097 *>          SCALE is REAL array, dimension (N)
00098 *>          Details of the permutations and scaling factors applied to
00099 *>          A.  If P(j) is the index of the row and column interchanged
00100 *>          with row and column j and D(j) is the scaling factor
00101 *>          applied to row and column j, then
00102 *>          SCALE(j) = P(j)    for j = 1,...,ILO-1
00103 *>                   = D(j)    for j = ILO,...,IHI
00104 *>                   = P(j)    for j = IHI+1,...,N.
00105 *>          The order in which the interchanges are made is N to IHI+1,
00106 *>          then 1 to ILO-1.
00107 *> \endverbatim
00108 *>
00109 *> \param[out] INFO
00110 *> \verbatim
00111 *>          INFO is INTEGER
00112 *>          = 0:  successful exit.
00113 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00114 *> \endverbatim
00115 *
00116 *  Authors:
00117 *  ========
00118 *
00119 *> \author Univ. of Tennessee 
00120 *> \author Univ. of California Berkeley 
00121 *> \author Univ. of Colorado Denver 
00122 *> \author NAG Ltd. 
00123 *
00124 *> \date November 2011
00125 *
00126 *> \ingroup realGEcomputational
00127 *
00128 *> \par Further Details:
00129 *  =====================
00130 *>
00131 *> \verbatim
00132 *>
00133 *>  The permutations consist of row and column interchanges which put
00134 *>  the matrix in the form
00135 *>
00136 *>             ( T1   X   Y  )
00137 *>     P A P = (  0   B   Z  )
00138 *>             (  0   0   T2 )
00139 *>
00140 *>  where T1 and T2 are upper triangular matrices whose eigenvalues lie
00141 *>  along the diagonal.  The column indices ILO and IHI mark the starting
00142 *>  and ending columns of the submatrix B. Balancing consists of applying
00143 *>  a diagonal similarity transformation inv(D) * B * D to make the
00144 *>  1-norms of each row of B and its corresponding column nearly equal.
00145 *>  The output matrix is
00146 *>
00147 *>     ( T1     X*D          Y    )
00148 *>     (  0  inv(D)*B*D  inv(D)*Z ).
00149 *>     (  0      0           T2   )
00150 *>
00151 *>  Information about the permutations P and the diagonal matrix D is
00152 *>  returned in the vector SCALE.
00153 *>
00154 *>  This subroutine is based on the EISPACK routine BALANC.
00155 *>
00156 *>  Modified by Tzu-Yi Chen, Computer Science Division, University of
00157 *>    California at Berkeley, USA
00158 *> \endverbatim
00159 *>
00160 *  =====================================================================
00161       SUBROUTINE SGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
00162 *
00163 *  -- LAPACK computational routine (version 3.4.0) --
00164 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00165 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00166 *     November 2011
00167 *
00168 *     .. Scalar Arguments ..
00169       CHARACTER          JOB
00170       INTEGER            IHI, ILO, INFO, LDA, N
00171 *     ..
00172 *     .. Array Arguments ..
00173       REAL               A( LDA, * ), SCALE( * )
00174 *     ..
00175 *
00176 *  =====================================================================
00177 *
00178 *     .. Parameters ..
00179       REAL               ZERO, ONE
00180       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00181       REAL               SCLFAC
00182       PARAMETER          ( SCLFAC = 2.0E+0 )
00183       REAL               FACTOR
00184       PARAMETER          ( FACTOR = 0.95E+0 )
00185 *     ..
00186 *     .. Local Scalars ..
00187       LOGICAL            NOCONV
00188       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
00189       REAL               C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
00190      $                   SFMIN2
00191 *     ..
00192 *     .. External Functions ..
00193       LOGICAL            SISNAN, LSAME
00194       INTEGER            ISAMAX
00195       REAL               SLAMCH
00196       EXTERNAL           SISNAN, LSAME, ISAMAX, SLAMCH
00197 *     ..
00198 *     .. External Subroutines ..
00199       EXTERNAL           SSCAL, SSWAP, XERBLA
00200 *     ..
00201 *     .. Intrinsic Functions ..
00202       INTRINSIC          ABS, MAX, MIN
00203 *     ..
00204 *     .. Executable Statements ..
00205 *
00206 *     Test the input parameters
00207 *
00208       INFO = 0
00209       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
00210      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
00211          INFO = -1
00212       ELSE IF( N.LT.0 ) THEN
00213          INFO = -2
00214       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00215          INFO = -4
00216       END IF
00217       IF( INFO.NE.0 ) THEN
00218          CALL XERBLA( 'SGEBAL', -INFO )
00219          RETURN
00220       END IF
00221 *
00222       K = 1
00223       L = N
00224 *
00225       IF( N.EQ.0 )
00226      $   GO TO 210
00227 *
00228       IF( LSAME( JOB, 'N' ) ) THEN
00229          DO 10 I = 1, N
00230             SCALE( I ) = ONE
00231    10    CONTINUE
00232          GO TO 210
00233       END IF
00234 *
00235       IF( LSAME( JOB, 'S' ) )
00236      $   GO TO 120
00237 *
00238 *     Permutation to isolate eigenvalues if possible
00239 *
00240       GO TO 50
00241 *
00242 *     Row and column exchange.
00243 *
00244    20 CONTINUE
00245       SCALE( M ) = J
00246       IF( J.EQ.M )
00247      $   GO TO 30
00248 *
00249       CALL SSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
00250       CALL SSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
00251 *
00252    30 CONTINUE
00253       GO TO ( 40, 80 )IEXC
00254 *
00255 *     Search for rows isolating an eigenvalue and push them down.
00256 *
00257    40 CONTINUE
00258       IF( L.EQ.1 )
00259      $   GO TO 210
00260       L = L - 1
00261 *
00262    50 CONTINUE
00263       DO 70 J = L, 1, -1
00264 *
00265          DO 60 I = 1, L
00266             IF( I.EQ.J )
00267      $         GO TO 60
00268             IF( A( J, I ).NE.ZERO )
00269      $         GO TO 70
00270    60    CONTINUE
00271 *
00272          M = L
00273          IEXC = 1
00274          GO TO 20
00275    70 CONTINUE
00276 *
00277       GO TO 90
00278 *
00279 *     Search for columns isolating an eigenvalue and push them left.
00280 *
00281    80 CONTINUE
00282       K = K + 1
00283 *
00284    90 CONTINUE
00285       DO 110 J = K, L
00286 *
00287          DO 100 I = K, L
00288             IF( I.EQ.J )
00289      $         GO TO 100
00290             IF( A( I, J ).NE.ZERO )
00291      $         GO TO 110
00292   100    CONTINUE
00293 *
00294          M = K
00295          IEXC = 2
00296          GO TO 20
00297   110 CONTINUE
00298 *
00299   120 CONTINUE
00300       DO 130 I = K, L
00301          SCALE( I ) = ONE
00302   130 CONTINUE
00303 *
00304       IF( LSAME( JOB, 'P' ) )
00305      $   GO TO 210
00306 *
00307 *     Balance the submatrix in rows K to L.
00308 *
00309 *     Iterative loop for norm reduction
00310 *
00311       SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
00312       SFMAX1 = ONE / SFMIN1
00313       SFMIN2 = SFMIN1*SCLFAC
00314       SFMAX2 = ONE / SFMIN2
00315   140 CONTINUE
00316       NOCONV = .FALSE.
00317 *
00318       DO 200 I = K, L
00319          C = ZERO
00320          R = ZERO
00321 *
00322          DO 150 J = K, L
00323             IF( J.EQ.I )
00324      $         GO TO 150
00325             C = C + ABS( A( J, I ) )
00326             R = R + ABS( A( I, J ) )
00327   150    CONTINUE
00328          ICA = ISAMAX( L, A( 1, I ), 1 )
00329          CA = ABS( A( ICA, I ) )
00330          IRA = ISAMAX( N-K+1, A( I, K ), LDA )
00331          RA = ABS( A( I, IRA+K-1 ) )
00332 *
00333 *        Guard against zero C or R due to underflow.
00334 *
00335          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
00336      $      GO TO 200
00337          G = R / SCLFAC
00338          F = ONE
00339          S = C + R
00340   160    CONTINUE
00341          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
00342      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
00343          F = F*SCLFAC
00344          C = C*SCLFAC
00345          CA = CA*SCLFAC
00346          R = R / SCLFAC
00347          G = G / SCLFAC
00348          RA = RA / SCLFAC
00349          GO TO 160
00350 *
00351   170    CONTINUE
00352          G = C / SCLFAC
00353   180    CONTINUE
00354          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
00355      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
00356             IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
00357 *
00358 *           Exit if NaN to avoid infinite loop
00359 *
00360             INFO = -3
00361             CALL XERBLA( 'SGEBAL', -INFO )
00362             RETURN
00363          END IF
00364          F = F / SCLFAC
00365          C = C / SCLFAC
00366          G = G / SCLFAC
00367          CA = CA / SCLFAC
00368          R = R*SCLFAC
00369          RA = RA*SCLFAC
00370          GO TO 180
00371 *
00372 *        Now balance.
00373 *
00374   190    CONTINUE
00375          IF( ( C+R ).GE.FACTOR*S )
00376      $      GO TO 200
00377          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
00378             IF( F*SCALE( I ).LE.SFMIN1 )
00379      $         GO TO 200
00380          END IF
00381          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
00382             IF( SCALE( I ).GE.SFMAX1 / F )
00383      $         GO TO 200
00384          END IF
00385          G = ONE / F
00386          SCALE( I ) = SCALE( I )*F
00387          NOCONV = .TRUE.
00388 *
00389          CALL SSCAL( N-K+1, G, A( I, K ), LDA )
00390          CALL SSCAL( L, F, A( 1, I ), 1 )
00391 *
00392   200 CONTINUE
00393 *
00394       IF( NOCONV )
00395      $   GO TO 140
00396 *
00397   210 CONTINUE
00398       ILO = K
00399       IHI = L
00400 *
00401       RETURN
00402 *
00403 *     End of SGEBAL
00404 *
00405       END
 All Files Functions