#include "ESMF_LapackBlas.inc" !> \brief \b DGEBAL ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DGEBAL + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgebal.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgebal.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgebal.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) ! ! .. Scalar Arguments .. ! CHARACTER JOB ! INTEGER IHI, ILO, INFO, LDA, N ! .. ! .. Array Arguments .. ! DOUBLE PRECISION A( LDA, * ), SCALE( * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> 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. !> \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 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. !> \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 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. !> \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 June 2017 ! !> \ingroup doubleGEcomputational ! !> \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 DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) ! ! -- LAPACK computational routine (version 3.7.1) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! June 2017 ! ! .. Scalar Arguments .. CHARACTER JOB INTEGER IHI, ILO, INFO, LDA, N ! .. ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), SCALE( * ) ! .. ! ! ===================================================================== ! ! .. 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 DISNAN, LSAME INTEGER IDAMAX DOUBLE PRECISION DLAMCH, DNRM2 EXTERNAL DISNAN, LSAME, IDAMAX, DLAMCH, DNRM2 ! .. ! .. External Subroutines .. EXTERNAL DSCAL, DSWAP, 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( '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 = DNRM2( L-K+1, A( K, I ), 1 ) R = DNRM2( L-K+1, A( I, K ), LDA ) 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 IF( DISNAN( C+F+CA+R+G+RA ) ) THEN ! ! Exit if NaN to avoid infinite loop ! INFO = -3 CALL XERBLA( 'DGEBAL', -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. ! 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