![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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