Annotation of rpl/lapack/lapack/zgebal.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZGEBAL
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGEBAL + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebal.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebal.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebal.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * CHARACTER JOB
! 25: * INTEGER IHI, ILO, INFO, LDA, N
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * DOUBLE PRECISION SCALE( * )
! 29: * COMPLEX*16 A( LDA, * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> ZGEBAL balances a general complex matrix A. This involves, first,
! 39: *> permuting A by a similarity transformation to isolate eigenvalues
! 40: *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
! 41: *> diagonal; and second, applying a diagonal similarity transformation
! 42: *> to rows and columns ILO to IHI to make the rows and columns as
! 43: *> close in norm as possible. Both steps are optional.
! 44: *>
! 45: *> Balancing may reduce the 1-norm of the matrix, and improve the
! 46: *> accuracy of the computed eigenvalues and/or eigenvectors.
! 47: *> \endverbatim
! 48: *
! 49: * Arguments:
! 50: * ==========
! 51: *
! 52: *> \param[in] JOB
! 53: *> \verbatim
! 54: *> JOB is CHARACTER*1
! 55: *> Specifies the operations to be performed on A:
! 56: *> = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
! 57: *> for i = 1,...,N;
! 58: *> = 'P': permute only;
! 59: *> = 'S': scale only;
! 60: *> = 'B': both permute and scale.
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] N
! 64: *> \verbatim
! 65: *> N is INTEGER
! 66: *> The order of the matrix A. N >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in,out] A
! 70: *> \verbatim
! 71: *> A is COMPLEX*16 array, dimension (LDA,N)
! 72: *> On entry, the input matrix A.
! 73: *> On exit, A is overwritten by the balanced matrix.
! 74: *> If JOB = 'N', A is not referenced.
! 75: *> See Further Details.
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in] LDA
! 79: *> \verbatim
! 80: *> LDA is INTEGER
! 81: *> The leading dimension of the array A. LDA >= max(1,N).
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[out] ILO
! 85: *> \verbatim
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[out] IHI
! 89: *> \verbatim
! 90: *> ILO and IHI are set to INTEGER such that on exit
! 91: *> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
! 92: *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[out] SCALE
! 96: *> \verbatim
! 97: *> SCALE is DOUBLE PRECISION array, dimension (N)
! 98: *> Details of the permutations and scaling factors applied to
! 99: *> A. If P(j) is the index of the row and column interchanged
! 100: *> with row and column j and D(j) is the scaling factor
! 101: *> applied to row and column j, then
! 102: *> SCALE(j) = P(j) for j = 1,...,ILO-1
! 103: *> = D(j) for j = ILO,...,IHI
! 104: *> = P(j) for j = IHI+1,...,N.
! 105: *> The order in which the interchanges are made is N to IHI+1,
! 106: *> then 1 to ILO-1.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[out] INFO
! 110: *> \verbatim
! 111: *> INFO is INTEGER
! 112: *> = 0: successful exit.
! 113: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 114: *> \endverbatim
! 115: *
! 116: * Authors:
! 117: * ========
! 118: *
! 119: *> \author Univ. of Tennessee
! 120: *> \author Univ. of California Berkeley
! 121: *> \author Univ. of Colorado Denver
! 122: *> \author NAG Ltd.
! 123: *
! 124: *> \date November 2011
! 125: *
! 126: *> \ingroup complex16GEcomputational
! 127: *
! 128: *> \par Further Details:
! 129: * =====================
! 130: *>
! 131: *> \verbatim
! 132: *>
! 133: *> The permutations consist of row and column interchanges which put
! 134: *> the matrix in the form
! 135: *>
! 136: *> ( T1 X Y )
! 137: *> P A P = ( 0 B Z )
! 138: *> ( 0 0 T2 )
! 139: *>
! 140: *> where T1 and T2 are upper triangular matrices whose eigenvalues lie
! 141: *> along the diagonal. The column indices ILO and IHI mark the starting
! 142: *> and ending columns of the submatrix B. Balancing consists of applying
! 143: *> a diagonal similarity transformation inv(D) * B * D to make the
! 144: *> 1-norms of each row of B and its corresponding column nearly equal.
! 145: *> The output matrix is
! 146: *>
! 147: *> ( T1 X*D Y )
! 148: *> ( 0 inv(D)*B*D inv(D)*Z ).
! 149: *> ( 0 0 T2 )
! 150: *>
! 151: *> Information about the permutations P and the diagonal matrix D is
! 152: *> returned in the vector SCALE.
! 153: *>
! 154: *> This subroutine is based on the EISPACK routine CBAL.
! 155: *>
! 156: *> Modified by Tzu-Yi Chen, Computer Science Division, University of
! 157: *> California at Berkeley, USA
! 158: *> \endverbatim
! 159: *>
! 160: * =====================================================================
1.1 bertrand 161: SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
162: *
1.9 ! bertrand 163: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 164: * -- LAPACK is a software package provided by Univ. of Tennessee, --
165: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 166: * November 2011
1.1 bertrand 167: *
168: * .. Scalar Arguments ..
169: CHARACTER JOB
170: INTEGER IHI, ILO, INFO, LDA, N
171: * ..
172: * .. Array Arguments ..
173: DOUBLE PRECISION SCALE( * )
174: COMPLEX*16 A( LDA, * )
175: * ..
176: *
177: * =====================================================================
178: *
179: * .. Parameters ..
180: DOUBLE PRECISION ZERO, ONE
181: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
182: DOUBLE PRECISION SCLFAC
183: PARAMETER ( SCLFAC = 2.0D+0 )
184: DOUBLE PRECISION FACTOR
185: PARAMETER ( FACTOR = 0.95D+0 )
186: * ..
187: * .. Local Scalars ..
188: LOGICAL NOCONV
189: INTEGER I, ICA, IEXC, IRA, J, K, L, M
190: DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
191: $ SFMIN2
192: COMPLEX*16 CDUM
193: * ..
194: * .. External Functions ..
1.5 bertrand 195: LOGICAL DISNAN, LSAME
1.1 bertrand 196: INTEGER IZAMAX
197: DOUBLE PRECISION DLAMCH
1.5 bertrand 198: EXTERNAL DISNAN, LSAME, IZAMAX, DLAMCH
1.1 bertrand 199: * ..
200: * .. External Subroutines ..
201: EXTERNAL XERBLA, ZDSCAL, ZSWAP
202: * ..
203: * .. Intrinsic Functions ..
204: INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
205: * ..
206: * .. Statement Functions ..
207: DOUBLE PRECISION CABS1
208: * ..
209: * .. Statement Function definitions ..
210: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
211: * ..
212: * .. Executable Statements ..
213: *
214: * Test the input parameters
215: *
216: INFO = 0
217: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
218: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
219: INFO = -1
220: ELSE IF( N.LT.0 ) THEN
221: INFO = -2
222: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
223: INFO = -4
224: END IF
225: IF( INFO.NE.0 ) THEN
226: CALL XERBLA( 'ZGEBAL', -INFO )
227: RETURN
228: END IF
229: *
230: K = 1
231: L = N
232: *
233: IF( N.EQ.0 )
234: $ GO TO 210
235: *
236: IF( LSAME( JOB, 'N' ) ) THEN
237: DO 10 I = 1, N
238: SCALE( I ) = ONE
239: 10 CONTINUE
240: GO TO 210
241: END IF
242: *
243: IF( LSAME( JOB, 'S' ) )
244: $ GO TO 120
245: *
246: * Permutation to isolate eigenvalues if possible
247: *
248: GO TO 50
249: *
250: * Row and column exchange.
251: *
252: 20 CONTINUE
253: SCALE( M ) = J
254: IF( J.EQ.M )
255: $ GO TO 30
256: *
257: CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
258: CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
259: *
260: 30 CONTINUE
261: GO TO ( 40, 80 )IEXC
262: *
263: * Search for rows isolating an eigenvalue and push them down.
264: *
265: 40 CONTINUE
266: IF( L.EQ.1 )
267: $ GO TO 210
268: L = L - 1
269: *
270: 50 CONTINUE
271: DO 70 J = L, 1, -1
272: *
273: DO 60 I = 1, L
274: IF( I.EQ.J )
275: $ GO TO 60
276: IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
277: $ ZERO )GO TO 70
278: 60 CONTINUE
279: *
280: M = L
281: IEXC = 1
282: GO TO 20
283: 70 CONTINUE
284: *
285: GO TO 90
286: *
287: * Search for columns isolating an eigenvalue and push them left.
288: *
289: 80 CONTINUE
290: K = K + 1
291: *
292: 90 CONTINUE
293: DO 110 J = K, L
294: *
295: DO 100 I = K, L
296: IF( I.EQ.J )
297: $ GO TO 100
298: IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
299: $ ZERO )GO TO 110
300: 100 CONTINUE
301: *
302: M = K
303: IEXC = 2
304: GO TO 20
305: 110 CONTINUE
306: *
307: 120 CONTINUE
308: DO 130 I = K, L
309: SCALE( I ) = ONE
310: 130 CONTINUE
311: *
312: IF( LSAME( JOB, 'P' ) )
313: $ GO TO 210
314: *
315: * Balance the submatrix in rows K to L.
316: *
317: * Iterative loop for norm reduction
318: *
319: SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
320: SFMAX1 = ONE / SFMIN1
321: SFMIN2 = SFMIN1*SCLFAC
322: SFMAX2 = ONE / SFMIN2
323: 140 CONTINUE
324: NOCONV = .FALSE.
325: *
326: DO 200 I = K, L
327: C = ZERO
328: R = ZERO
329: *
330: DO 150 J = K, L
331: IF( J.EQ.I )
332: $ GO TO 150
333: C = C + CABS1( A( J, I ) )
334: R = R + CABS1( A( I, J ) )
335: 150 CONTINUE
336: ICA = IZAMAX( L, A( 1, I ), 1 )
337: CA = ABS( A( ICA, I ) )
338: IRA = IZAMAX( N-K+1, A( I, K ), LDA )
339: RA = ABS( A( I, IRA+K-1 ) )
340: *
341: * Guard against zero C or R due to underflow.
342: *
343: IF( C.EQ.ZERO .OR. R.EQ.ZERO )
344: $ GO TO 200
345: G = R / SCLFAC
346: F = ONE
347: S = C + R
348: 160 CONTINUE
349: IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
350: $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
1.5 bertrand 351: IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
352: *
353: * Exit if NaN to avoid infinite loop
354: *
355: INFO = -3
356: CALL XERBLA( 'ZGEBAL', -INFO )
357: RETURN
358: END IF
1.1 bertrand 359: F = F*SCLFAC
360: C = C*SCLFAC
361: CA = CA*SCLFAC
362: R = R / SCLFAC
363: G = G / SCLFAC
364: RA = RA / SCLFAC
365: GO TO 160
366: *
367: 170 CONTINUE
368: G = C / SCLFAC
369: 180 CONTINUE
370: IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
371: $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
372: F = F / SCLFAC
373: C = C / SCLFAC
374: G = G / SCLFAC
375: CA = CA / SCLFAC
376: R = R*SCLFAC
377: RA = RA*SCLFAC
378: GO TO 180
379: *
380: * Now balance.
381: *
382: 190 CONTINUE
383: IF( ( C+R ).GE.FACTOR*S )
384: $ GO TO 200
385: IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
386: IF( F*SCALE( I ).LE.SFMIN1 )
387: $ GO TO 200
388: END IF
389: IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
390: IF( SCALE( I ).GE.SFMAX1 / F )
391: $ GO TO 200
392: END IF
393: G = ONE / F
394: SCALE( I ) = SCALE( I )*F
395: NOCONV = .TRUE.
396: *
397: CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
398: CALL ZDSCAL( L, F, A( 1, I ), 1 )
399: *
400: 200 CONTINUE
401: *
402: IF( NOCONV )
403: $ GO TO 140
404: *
405: 210 CONTINUE
406: ILO = K
407: IHI = L
408: *
409: RETURN
410: *
411: * End of ZGEBAL
412: *
413: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgebal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack