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: *> ILO is INTEGER
87: *> \endverbatim
88: *>
89: *> \param[out] IHI
90: *> \verbatim
91: *> IHI is INTEGER
92: *> ILO and IHI are set to INTEGER such that on exit
93: *> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
94: *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
95: *> \endverbatim
96: *>
97: *> \param[out] SCALE
98: *> \verbatim
99: *> SCALE is DOUBLE PRECISION array, dimension (N)
100: *> Details of the permutations and scaling factors applied to
101: *> A. If P(j) is the index of the row and column interchanged
102: *> with row and column j and D(j) is the scaling factor
103: *> applied to row and column j, then
104: *> SCALE(j) = P(j) for j = 1,...,ILO-1
105: *> = D(j) for j = ILO,...,IHI
106: *> = P(j) for j = IHI+1,...,N.
107: *> The order in which the interchanges are made is N to IHI+1,
108: *> then 1 to ILO-1.
109: *> \endverbatim
110: *>
111: *> \param[out] INFO
112: *> \verbatim
113: *> INFO is INTEGER
114: *> = 0: successful exit.
115: *> < 0: if INFO = -i, the i-th argument had an illegal value.
116: *> \endverbatim
117: *
118: * Authors:
119: * ========
120: *
121: *> \author Univ. of Tennessee
122: *> \author Univ. of California Berkeley
123: *> \author Univ. of Colorado Denver
124: *> \author NAG Ltd.
125: *
126: *> \date June 2017
127: *
128: *> \ingroup complex16GEcomputational
129: *
130: *> \par Further Details:
131: * =====================
132: *>
133: *> \verbatim
134: *>
135: *> The permutations consist of row and column interchanges which put
136: *> the matrix in the form
137: *>
138: *> ( T1 X Y )
139: *> P A P = ( 0 B Z )
140: *> ( 0 0 T2 )
141: *>
142: *> where T1 and T2 are upper triangular matrices whose eigenvalues lie
143: *> along the diagonal. The column indices ILO and IHI mark the starting
144: *> and ending columns of the submatrix B. Balancing consists of applying
145: *> a diagonal similarity transformation inv(D) * B * D to make the
146: *> 1-norms of each row of B and its corresponding column nearly equal.
147: *> The output matrix is
148: *>
149: *> ( T1 X*D Y )
150: *> ( 0 inv(D)*B*D inv(D)*Z ).
151: *> ( 0 0 T2 )
152: *>
153: *> Information about the permutations P and the diagonal matrix D is
154: *> returned in the vector SCALE.
155: *>
156: *> This subroutine is based on the EISPACK routine CBAL.
157: *>
158: *> Modified by Tzu-Yi Chen, Computer Science Division, University of
159: *> California at Berkeley, USA
160: *> \endverbatim
161: *>
162: * =====================================================================
163: SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
164: *
165: * -- LAPACK computational routine (version 3.7.1) --
166: * -- LAPACK is a software package provided by Univ. of Tennessee, --
167: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168: * June 2017
169: *
170: * .. Scalar Arguments ..
171: CHARACTER JOB
172: INTEGER IHI, ILO, INFO, LDA, N
173: * ..
174: * .. Array Arguments ..
175: DOUBLE PRECISION SCALE( * )
176: COMPLEX*16 A( LDA, * )
177: * ..
178: *
179: * =====================================================================
180: *
181: * .. Parameters ..
182: DOUBLE PRECISION ZERO, ONE
183: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
184: DOUBLE PRECISION SCLFAC
185: PARAMETER ( SCLFAC = 2.0D+0 )
186: DOUBLE PRECISION FACTOR
187: PARAMETER ( FACTOR = 0.95D+0 )
188: * ..
189: * .. Local Scalars ..
190: LOGICAL NOCONV
191: INTEGER I, ICA, IEXC, IRA, J, K, L, M
192: DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
193: $ SFMIN2
194: * ..
195: * .. External Functions ..
196: LOGICAL DISNAN, LSAME
197: INTEGER IZAMAX
198: DOUBLE PRECISION DLAMCH, DZNRM2
199: EXTERNAL DISNAN, LSAME, IZAMAX, DLAMCH, DZNRM2
200: * ..
201: * .. External Subroutines ..
202: EXTERNAL XERBLA, ZDSCAL, ZSWAP
203: * ..
204: * .. Intrinsic Functions ..
205: INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
206: *
207: * Test the input parameters
208: *
209: INFO = 0
210: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
211: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
212: INFO = -1
213: ELSE IF( N.LT.0 ) THEN
214: INFO = -2
215: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
216: INFO = -4
217: END IF
218: IF( INFO.NE.0 ) THEN
219: CALL XERBLA( 'ZGEBAL', -INFO )
220: RETURN
221: END IF
222: *
223: K = 1
224: L = N
225: *
226: IF( N.EQ.0 )
227: $ GO TO 210
228: *
229: IF( LSAME( JOB, 'N' ) ) THEN
230: DO 10 I = 1, N
231: SCALE( I ) = ONE
232: 10 CONTINUE
233: GO TO 210
234: END IF
235: *
236: IF( LSAME( JOB, 'S' ) )
237: $ GO TO 120
238: *
239: * Permutation to isolate eigenvalues if possible
240: *
241: GO TO 50
242: *
243: * Row and column exchange.
244: *
245: 20 CONTINUE
246: SCALE( M ) = J
247: IF( J.EQ.M )
248: $ GO TO 30
249: *
250: CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
251: CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
252: *
253: 30 CONTINUE
254: GO TO ( 40, 80 )IEXC
255: *
256: * Search for rows isolating an eigenvalue and push them down.
257: *
258: 40 CONTINUE
259: IF( L.EQ.1 )
260: $ GO TO 210
261: L = L - 1
262: *
263: 50 CONTINUE
264: DO 70 J = L, 1, -1
265: *
266: DO 60 I = 1, L
267: IF( I.EQ.J )
268: $ GO TO 60
269: IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
270: $ ZERO )GO TO 70
271: 60 CONTINUE
272: *
273: M = L
274: IEXC = 1
275: GO TO 20
276: 70 CONTINUE
277: *
278: GO TO 90
279: *
280: * Search for columns isolating an eigenvalue and push them left.
281: *
282: 80 CONTINUE
283: K = K + 1
284: *
285: 90 CONTINUE
286: DO 110 J = K, L
287: *
288: DO 100 I = K, L
289: IF( I.EQ.J )
290: $ GO TO 100
291: IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
292: $ ZERO )GO TO 110
293: 100 CONTINUE
294: *
295: M = K
296: IEXC = 2
297: GO TO 20
298: 110 CONTINUE
299: *
300: 120 CONTINUE
301: DO 130 I = K, L
302: SCALE( I ) = ONE
303: 130 CONTINUE
304: *
305: IF( LSAME( JOB, 'P' ) )
306: $ GO TO 210
307: *
308: * Balance the submatrix in rows K to L.
309: *
310: * Iterative loop for norm reduction
311: *
312: SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
313: SFMAX1 = ONE / SFMIN1
314: SFMIN2 = SFMIN1*SCLFAC
315: SFMAX2 = ONE / SFMIN2
316: 140 CONTINUE
317: NOCONV = .FALSE.
318: *
319: DO 200 I = K, L
320: *
321: C = DZNRM2( L-K+1, A( K, I ), 1 )
322: R = DZNRM2( L-K+1, A( I, K ), LDA )
323: ICA = IZAMAX( L, A( 1, I ), 1 )
324: CA = ABS( A( ICA, I ) )
325: IRA = IZAMAX( N-K+1, A( I, K ), LDA )
326: RA = ABS( A( I, IRA+K-1 ) )
327: *
328: * Guard against zero C or R due to underflow.
329: *
330: IF( C.EQ.ZERO .OR. R.EQ.ZERO )
331: $ GO TO 200
332: G = R / SCLFAC
333: F = ONE
334: S = C + R
335: 160 CONTINUE
336: IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
337: $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
338: IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
339: *
340: * Exit if NaN to avoid infinite loop
341: *
342: INFO = -3
343: CALL XERBLA( 'ZGEBAL', -INFO )
344: RETURN
345: END IF
346: F = F*SCLFAC
347: C = C*SCLFAC
348: CA = CA*SCLFAC
349: R = R / SCLFAC
350: G = G / SCLFAC
351: RA = RA / SCLFAC
352: GO TO 160
353: *
354: 170 CONTINUE
355: G = C / SCLFAC
356: 180 CONTINUE
357: IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
358: $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
359: F = F / SCLFAC
360: C = C / SCLFAC
361: G = G / SCLFAC
362: CA = CA / SCLFAC
363: R = R*SCLFAC
364: RA = RA*SCLFAC
365: GO TO 180
366: *
367: * Now balance.
368: *
369: 190 CONTINUE
370: IF( ( C+R ).GE.FACTOR*S )
371: $ GO TO 200
372: IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
373: IF( F*SCALE( I ).LE.SFMIN1 )
374: $ GO TO 200
375: END IF
376: IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
377: IF( SCALE( I ).GE.SFMAX1 / F )
378: $ GO TO 200
379: END IF
380: G = ONE / F
381: SCALE( I ) = SCALE( I )*F
382: NOCONV = .TRUE.
383: *
384: CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
385: CALL ZDSCAL( L, F, A( 1, I ), 1 )
386: *
387: 200 CONTINUE
388: *
389: IF( NOCONV )
390: $ GO TO 140
391: *
392: 210 CONTINUE
393: ILO = K
394: IHI = L
395: *
396: RETURN
397: *
398: * End of ZGEBAL
399: *
400: 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