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