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