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