1: SUBROUTINE ZGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
2: $ AMAX, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, KL, KU, LDAB, M, N
11: DOUBLE PRECISION AMAX, COLCND, ROWCND
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION C( * ), R( * )
15: COMPLEX*16 AB( LDAB, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZGBEQU computes row and column scalings intended to equilibrate an
22: * M-by-N band matrix A and reduce its condition number. R returns the
23: * row scale factors and C the column scale factors, chosen to try to
24: * make the largest element in each row and column of the matrix B with
25: * elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
26: *
27: * R(i) and C(j) are restricted to be between SMLNUM = smallest safe
28: * number and BIGNUM = largest safe number. Use of these scaling
29: * factors is not guaranteed to reduce the condition number of A but
30: * works well in practice.
31: *
32: * Arguments
33: * =========
34: *
35: * M (input) INTEGER
36: * The number of rows of the matrix A. M >= 0.
37: *
38: * N (input) INTEGER
39: * The number of columns of the matrix A. N >= 0.
40: *
41: * KL (input) INTEGER
42: * The number of subdiagonals within the band of A. KL >= 0.
43: *
44: * KU (input) INTEGER
45: * The number of superdiagonals within the band of A. KU >= 0.
46: *
47: * AB (input) COMPLEX*16 array, dimension (LDAB,N)
48: * The band matrix A, stored in rows 1 to KL+KU+1. The j-th
49: * column of A is stored in the j-th column of the array AB as
50: * follows:
51: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
52: *
53: * LDAB (input) INTEGER
54: * The leading dimension of the array AB. LDAB >= KL+KU+1.
55: *
56: * R (output) DOUBLE PRECISION array, dimension (M)
57: * If INFO = 0, or INFO > M, R contains the row scale factors
58: * for A.
59: *
60: * C (output) DOUBLE PRECISION array, dimension (N)
61: * If INFO = 0, C contains the column scale factors for A.
62: *
63: * ROWCND (output) DOUBLE PRECISION
64: * If INFO = 0 or INFO > M, ROWCND contains the ratio of the
65: * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
66: * AMAX is neither too large nor too small, it is not worth
67: * scaling by R.
68: *
69: * COLCND (output) DOUBLE PRECISION
70: * If INFO = 0, COLCND contains the ratio of the smallest
71: * C(i) to the largest C(i). If COLCND >= 0.1, it is not
72: * worth scaling by C.
73: *
74: * AMAX (output) DOUBLE PRECISION
75: * Absolute value of largest matrix element. If AMAX is very
76: * close to overflow or very close to underflow, the matrix
77: * should be scaled.
78: *
79: * INFO (output) INTEGER
80: * = 0: successful exit
81: * < 0: if INFO = -i, the i-th argument had an illegal value
82: * > 0: if INFO = i, and i is
83: * <= M: the i-th row of A is exactly zero
84: * > M: the (i-M)-th column of A is exactly zero
85: *
86: * =====================================================================
87: *
88: * .. Parameters ..
89: DOUBLE PRECISION ONE, ZERO
90: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
91: * ..
92: * .. Local Scalars ..
93: INTEGER I, J, KD
94: DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM
95: COMPLEX*16 ZDUM
96: * ..
97: * .. External Functions ..
98: DOUBLE PRECISION DLAMCH
99: EXTERNAL DLAMCH
100: * ..
101: * .. External Subroutines ..
102: EXTERNAL XERBLA
103: * ..
104: * .. Intrinsic Functions ..
105: INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
106: * ..
107: * .. Statement Functions ..
108: DOUBLE PRECISION CABS1
109: * ..
110: * .. Statement Function definitions ..
111: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
112: * ..
113: * .. Executable Statements ..
114: *
115: * Test the input parameters
116: *
117: INFO = 0
118: IF( M.LT.0 ) THEN
119: INFO = -1
120: ELSE IF( N.LT.0 ) THEN
121: INFO = -2
122: ELSE IF( KL.LT.0 ) THEN
123: INFO = -3
124: ELSE IF( KU.LT.0 ) THEN
125: INFO = -4
126: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
127: INFO = -6
128: END IF
129: IF( INFO.NE.0 ) THEN
130: CALL XERBLA( 'ZGBEQU', -INFO )
131: RETURN
132: END IF
133: *
134: * Quick return if possible
135: *
136: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
137: ROWCND = ONE
138: COLCND = ONE
139: AMAX = ZERO
140: RETURN
141: END IF
142: *
143: * Get machine constants.
144: *
145: SMLNUM = DLAMCH( 'S' )
146: BIGNUM = ONE / SMLNUM
147: *
148: * Compute row scale factors.
149: *
150: DO 10 I = 1, M
151: R( I ) = ZERO
152: 10 CONTINUE
153: *
154: * Find the maximum element in each row.
155: *
156: KD = KU + 1
157: DO 30 J = 1, N
158: DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
159: R( I ) = MAX( R( I ), CABS1( AB( KD+I-J, J ) ) )
160: 20 CONTINUE
161: 30 CONTINUE
162: *
163: * Find the maximum and minimum scale factors.
164: *
165: RCMIN = BIGNUM
166: RCMAX = ZERO
167: DO 40 I = 1, M
168: RCMAX = MAX( RCMAX, R( I ) )
169: RCMIN = MIN( RCMIN, R( I ) )
170: 40 CONTINUE
171: AMAX = RCMAX
172: *
173: IF( RCMIN.EQ.ZERO ) THEN
174: *
175: * Find the first zero scale factor and return an error code.
176: *
177: DO 50 I = 1, M
178: IF( R( I ).EQ.ZERO ) THEN
179: INFO = I
180: RETURN
181: END IF
182: 50 CONTINUE
183: ELSE
184: *
185: * Invert the scale factors.
186: *
187: DO 60 I = 1, M
188: R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
189: 60 CONTINUE
190: *
191: * Compute ROWCND = min(R(I)) / max(R(I))
192: *
193: ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
194: END IF
195: *
196: * Compute column scale factors
197: *
198: DO 70 J = 1, N
199: C( J ) = ZERO
200: 70 CONTINUE
201: *
202: * Find the maximum element in each column,
203: * assuming the row scaling computed above.
204: *
205: KD = KU + 1
206: DO 90 J = 1, N
207: DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
208: C( J ) = MAX( C( J ), CABS1( AB( KD+I-J, J ) )*R( I ) )
209: 80 CONTINUE
210: 90 CONTINUE
211: *
212: * Find the maximum and minimum scale factors.
213: *
214: RCMIN = BIGNUM
215: RCMAX = ZERO
216: DO 100 J = 1, N
217: RCMIN = MIN( RCMIN, C( J ) )
218: RCMAX = MAX( RCMAX, C( J ) )
219: 100 CONTINUE
220: *
221: IF( RCMIN.EQ.ZERO ) THEN
222: *
223: * Find the first zero scale factor and return an error code.
224: *
225: DO 110 J = 1, N
226: IF( C( J ).EQ.ZERO ) THEN
227: INFO = M + J
228: RETURN
229: END IF
230: 110 CONTINUE
231: ELSE
232: *
233: * Invert the scale factors.
234: *
235: DO 120 J = 1, N
236: C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
237: 120 CONTINUE
238: *
239: * Compute COLCND = min(C(J)) / max(C(J))
240: *
241: COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
242: END IF
243: *
244: RETURN
245: *
246: * End of ZGBEQU
247: *
248: END
CVSweb interface <joel.bertrand@systella.fr>