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