![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGEBAL 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGEBAL + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebal.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebal.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebal.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGEBAL( 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 * DOUBLE PRECISION SCALE( * ) 00029 * COMPLEX*16 A( LDA, * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> ZGEBAL 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*16 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 *> \endverbatim 00087 *> 00088 *> \param[out] IHI 00089 *> \verbatim 00090 *> ILO and IHI are set to INTEGER 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 DOUBLE PRECISION 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 complex16GEcomputational 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 CBAL. 00155 *> 00156 *> Modified by Tzu-Yi Chen, Computer Science Division, University of 00157 *> California at Berkeley, USA 00158 *> \endverbatim 00159 *> 00160 * ===================================================================== 00161 SUBROUTINE ZGEBAL( 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 DOUBLE PRECISION SCALE( * ) 00174 COMPLEX*16 A( LDA, * ) 00175 * .. 00176 * 00177 * ===================================================================== 00178 * 00179 * .. Parameters .. 00180 DOUBLE PRECISION ZERO, ONE 00181 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00182 DOUBLE PRECISION SCLFAC 00183 PARAMETER ( SCLFAC = 2.0D+0 ) 00184 DOUBLE PRECISION FACTOR 00185 PARAMETER ( FACTOR = 0.95D+0 ) 00186 * .. 00187 * .. Local Scalars .. 00188 LOGICAL NOCONV 00189 INTEGER I, ICA, IEXC, IRA, J, K, L, M 00190 DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1, 00191 $ SFMIN2 00192 COMPLEX*16 CDUM 00193 * .. 00194 * .. External Functions .. 00195 LOGICAL DISNAN, LSAME 00196 INTEGER IZAMAX 00197 DOUBLE PRECISION DLAMCH 00198 EXTERNAL DISNAN, LSAME, IZAMAX, DLAMCH 00199 * .. 00200 * .. External Subroutines .. 00201 EXTERNAL XERBLA, ZDSCAL, ZSWAP 00202 * .. 00203 * .. Intrinsic Functions .. 00204 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN 00205 * .. 00206 * .. Statement Functions .. 00207 DOUBLE PRECISION CABS1 00208 * .. 00209 * .. Statement Function definitions .. 00210 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 00211 * .. 00212 * .. Executable Statements .. 00213 * 00214 * Test the input parameters 00215 * 00216 INFO = 0 00217 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 00218 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 00219 INFO = -1 00220 ELSE IF( N.LT.0 ) THEN 00221 INFO = -2 00222 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00223 INFO = -4 00224 END IF 00225 IF( INFO.NE.0 ) THEN 00226 CALL XERBLA( 'ZGEBAL', -INFO ) 00227 RETURN 00228 END IF 00229 * 00230 K = 1 00231 L = N 00232 * 00233 IF( N.EQ.0 ) 00234 $ GO TO 210 00235 * 00236 IF( LSAME( JOB, 'N' ) ) THEN 00237 DO 10 I = 1, N 00238 SCALE( I ) = ONE 00239 10 CONTINUE 00240 GO TO 210 00241 END IF 00242 * 00243 IF( LSAME( JOB, 'S' ) ) 00244 $ GO TO 120 00245 * 00246 * Permutation to isolate eigenvalues if possible 00247 * 00248 GO TO 50 00249 * 00250 * Row and column exchange. 00251 * 00252 20 CONTINUE 00253 SCALE( M ) = J 00254 IF( J.EQ.M ) 00255 $ GO TO 30 00256 * 00257 CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 00258 CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA ) 00259 * 00260 30 CONTINUE 00261 GO TO ( 40, 80 )IEXC 00262 * 00263 * Search for rows isolating an eigenvalue and push them down. 00264 * 00265 40 CONTINUE 00266 IF( L.EQ.1 ) 00267 $ GO TO 210 00268 L = L - 1 00269 * 00270 50 CONTINUE 00271 DO 70 J = L, 1, -1 00272 * 00273 DO 60 I = 1, L 00274 IF( I.EQ.J ) 00275 $ GO TO 60 00276 IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE. 00277 $ ZERO )GO TO 70 00278 60 CONTINUE 00279 * 00280 M = L 00281 IEXC = 1 00282 GO TO 20 00283 70 CONTINUE 00284 * 00285 GO TO 90 00286 * 00287 * Search for columns isolating an eigenvalue and push them left. 00288 * 00289 80 CONTINUE 00290 K = K + 1 00291 * 00292 90 CONTINUE 00293 DO 110 J = K, L 00294 * 00295 DO 100 I = K, L 00296 IF( I.EQ.J ) 00297 $ GO TO 100 00298 IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE. 00299 $ ZERO )GO TO 110 00300 100 CONTINUE 00301 * 00302 M = K 00303 IEXC = 2 00304 GO TO 20 00305 110 CONTINUE 00306 * 00307 120 CONTINUE 00308 DO 130 I = K, L 00309 SCALE( I ) = ONE 00310 130 CONTINUE 00311 * 00312 IF( LSAME( JOB, 'P' ) ) 00313 $ GO TO 210 00314 * 00315 * Balance the submatrix in rows K to L. 00316 * 00317 * Iterative loop for norm reduction 00318 * 00319 SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00320 SFMAX1 = ONE / SFMIN1 00321 SFMIN2 = SFMIN1*SCLFAC 00322 SFMAX2 = ONE / SFMIN2 00323 140 CONTINUE 00324 NOCONV = .FALSE. 00325 * 00326 DO 200 I = K, L 00327 C = ZERO 00328 R = ZERO 00329 * 00330 DO 150 J = K, L 00331 IF( J.EQ.I ) 00332 $ GO TO 150 00333 C = C + CABS1( A( J, I ) ) 00334 R = R + CABS1( A( I, J ) ) 00335 150 CONTINUE 00336 ICA = IZAMAX( L, A( 1, I ), 1 ) 00337 CA = ABS( A( ICA, I ) ) 00338 IRA = IZAMAX( N-K+1, A( I, K ), LDA ) 00339 RA = ABS( A( I, IRA+K-1 ) ) 00340 * 00341 * Guard against zero C or R due to underflow. 00342 * 00343 IF( C.EQ.ZERO .OR. R.EQ.ZERO ) 00344 $ GO TO 200 00345 G = R / SCLFAC 00346 F = ONE 00347 S = C + R 00348 160 CONTINUE 00349 IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR. 00350 $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170 00351 IF( DISNAN( C+F+CA+R+G+RA ) ) THEN 00352 * 00353 * Exit if NaN to avoid infinite loop 00354 * 00355 INFO = -3 00356 CALL XERBLA( 'ZGEBAL', -INFO ) 00357 RETURN 00358 END IF 00359 F = F*SCLFAC 00360 C = C*SCLFAC 00361 CA = CA*SCLFAC 00362 R = R / SCLFAC 00363 G = G / SCLFAC 00364 RA = RA / SCLFAC 00365 GO TO 160 00366 * 00367 170 CONTINUE 00368 G = C / SCLFAC 00369 180 CONTINUE 00370 IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR. 00371 $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190 00372 F = F / SCLFAC 00373 C = C / SCLFAC 00374 G = G / SCLFAC 00375 CA = CA / SCLFAC 00376 R = R*SCLFAC 00377 RA = RA*SCLFAC 00378 GO TO 180 00379 * 00380 * Now balance. 00381 * 00382 190 CONTINUE 00383 IF( ( C+R ).GE.FACTOR*S ) 00384 $ GO TO 200 00385 IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN 00386 IF( F*SCALE( I ).LE.SFMIN1 ) 00387 $ GO TO 200 00388 END IF 00389 IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN 00390 IF( SCALE( I ).GE.SFMAX1 / F ) 00391 $ GO TO 200 00392 END IF 00393 G = ONE / F 00394 SCALE( I ) = SCALE( I )*F 00395 NOCONV = .TRUE. 00396 * 00397 CALL ZDSCAL( N-K+1, G, A( I, K ), LDA ) 00398 CALL ZDSCAL( L, F, A( 1, I ), 1 ) 00399 * 00400 200 CONTINUE 00401 * 00402 IF( NOCONV ) 00403 $ GO TO 140 00404 * 00405 210 CONTINUE 00406 ILO = K 00407 IHI = L 00408 * 00409 RETURN 00410 * 00411 * End of ZGEBAL 00412 * 00413 END