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