1: *> \brief \b DSYEQUB
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYEQUB + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyequb.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyequb.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyequb.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, N
25: * DOUBLE PRECISION AMAX, SCOND
26: * CHARACTER UPLO
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), S( * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DSYEQUB computes row and column scalings intended to equilibrate a
39: *> symmetric matrix A and reduce its condition number
40: *> (with respect to the two-norm). S contains the scale factors,
41: *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
42: *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This
43: *> choice of S puts the condition number of B within a factor N of the
44: *> smallest possible condition number over all possible diagonal
45: *> scalings.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> Specifies whether the details of the factorization are stored
55: *> as an upper or lower triangular matrix.
56: *> = 'U': Upper triangular, form is A = U*D*U**T;
57: *> = 'L': Lower triangular, form is A = L*D*L**T.
58: *> \endverbatim
59: *>
60: *> \param[in] N
61: *> \verbatim
62: *> N is INTEGER
63: *> The order of the matrix A. N >= 0.
64: *> \endverbatim
65: *>
66: *> \param[in] A
67: *> \verbatim
68: *> A is DOUBLE PRECISION array, dimension (LDA,N)
69: *> The N-by-N symmetric matrix whose scaling
70: *> factors are to be computed. Only the diagonal elements of A
71: *> are referenced.
72: *> \endverbatim
73: *>
74: *> \param[in] LDA
75: *> \verbatim
76: *> LDA is INTEGER
77: *> The leading dimension of the array A. LDA >= max(1,N).
78: *> \endverbatim
79: *>
80: *> \param[out] S
81: *> \verbatim
82: *> S is DOUBLE PRECISION array, dimension (N)
83: *> If INFO = 0, S contains the scale factors for A.
84: *> \endverbatim
85: *>
86: *> \param[out] SCOND
87: *> \verbatim
88: *> SCOND is DOUBLE PRECISION
89: *> If INFO = 0, S contains the ratio of the smallest S(i) to
90: *> the largest S(i). If SCOND >= 0.1 and AMAX is neither too
91: *> large nor too small, it is not worth scaling by S.
92: *> \endverbatim
93: *>
94: *> \param[out] AMAX
95: *> \verbatim
96: *> AMAX is DOUBLE PRECISION
97: *> Absolute value of largest matrix element. If AMAX is very
98: *> close to overflow or very close to underflow, the matrix
99: *> should be scaled.
100: *> \endverbatim
101: *>
102: *> \param[out] WORK
103: *> \verbatim
104: *> WORK is DOUBLE PRECISION array, dimension (3*N)
105: *> \endverbatim
106: *>
107: *> \param[out] INFO
108: *> \verbatim
109: *> INFO is INTEGER
110: *> = 0: successful exit
111: *> < 0: if INFO = -i, the i-th argument had an illegal value
112: *> > 0: if INFO = i, the i-th diagonal element is nonpositive.
113: *> \endverbatim
114: *
115: * Authors:
116: * ========
117: *
118: *> \author Univ. of Tennessee
119: *> \author Univ. of California Berkeley
120: *> \author Univ. of Colorado Denver
121: *> \author NAG Ltd.
122: *
123: *> \date November 2011
124: *
125: *> \ingroup doubleSYcomputational
126: *
127: *> \par References:
128: * ================
129: *>
130: *> Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
131: *> Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
132: *> DOI 10.1023/B:NUMA.0000016606.32820.69 \n
133: *> Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf
134: *>
135: * =====================================================================
136: SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
137: *
138: * -- LAPACK computational routine (version 3.4.0) --
139: * -- LAPACK is a software package provided by Univ. of Tennessee, --
140: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
141: * November 2011
142: *
143: * .. Scalar Arguments ..
144: INTEGER INFO, LDA, N
145: DOUBLE PRECISION AMAX, SCOND
146: CHARACTER UPLO
147: * ..
148: * .. Array Arguments ..
149: DOUBLE PRECISION A( LDA, * ), S( * ), WORK( * )
150: * ..
151: *
152: * =====================================================================
153: *
154: * .. Parameters ..
155: DOUBLE PRECISION ONE, ZERO
156: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
157: INTEGER MAX_ITER
158: PARAMETER ( MAX_ITER = 100 )
159: * ..
160: * .. Local Scalars ..
161: INTEGER I, J, ITER
162: DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
163: $ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
164: LOGICAL UP
165: * ..
166: * .. External Functions ..
167: DOUBLE PRECISION DLAMCH
168: LOGICAL LSAME
169: EXTERNAL DLAMCH, LSAME
170: * ..
171: * .. External Subroutines ..
172: EXTERNAL DLASSQ
173: * ..
174: * .. Intrinsic Functions ..
175: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
176: * ..
177: * .. Executable Statements ..
178: *
179: * Test input parameters.
180: *
181: INFO = 0
182: IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
183: INFO = -1
184: ELSE IF ( N .LT. 0 ) THEN
185: INFO = -2
186: ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
187: INFO = -4
188: END IF
189: IF ( INFO .NE. 0 ) THEN
190: CALL XERBLA( 'DSYEQUB', -INFO )
191: RETURN
192: END IF
193:
194: UP = LSAME( UPLO, 'U' )
195: AMAX = ZERO
196: *
197: * Quick return if possible.
198: *
199: IF ( N .EQ. 0 ) THEN
200: SCOND = ONE
201: RETURN
202: END IF
203:
204: DO I = 1, N
205: S( I ) = ZERO
206: END DO
207:
208: AMAX = ZERO
209: IF ( UP ) THEN
210: DO J = 1, N
211: DO I = 1, J-1
212: S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
213: S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
214: AMAX = MAX( AMAX, ABS( A(I, J) ) )
215: END DO
216: S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
217: AMAX = MAX( AMAX, ABS( A( J, J ) ) )
218: END DO
219: ELSE
220: DO J = 1, N
221: S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
222: AMAX = MAX( AMAX, ABS( A( J, J ) ) )
223: DO I = J+1, N
224: S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
225: S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
226: AMAX = MAX( AMAX, ABS( A( I, J ) ) )
227: END DO
228: END DO
229: END IF
230: DO J = 1, N
231: S( J ) = 1.0D+0 / S( J )
232: END DO
233:
234: TOL = ONE / SQRT(2.0D0 * N)
235:
236: DO ITER = 1, MAX_ITER
237: SCALE = 0.0D+0
238: SUMSQ = 0.0D+0
239: * BETA = |A|S
240: DO I = 1, N
241: WORK(I) = ZERO
242: END DO
243: IF ( UP ) THEN
244: DO J = 1, N
245: DO I = 1, J-1
246: T = ABS( A( I, J ) )
247: WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
248: WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
249: END DO
250: WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
251: END DO
252: ELSE
253: DO J = 1, N
254: WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
255: DO I = J+1, N
256: T = ABS( A( I, J ) )
257: WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
258: WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
259: END DO
260: END DO
261: END IF
262:
263: * avg = s^T beta / n
264: AVG = 0.0D+0
265: DO I = 1, N
266: AVG = AVG + S( I )*WORK( I )
267: END DO
268: AVG = AVG / N
269:
270: STD = 0.0D+0
271: DO I = 2*N+1, 3*N
272: WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG
273: END DO
274: CALL DLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ )
275: STD = SCALE * SQRT( SUMSQ / N )
276:
277: IF ( STD .LT. TOL * AVG ) GOTO 999
278:
279: DO I = 1, N
280: T = ABS( A( I, I ) )
281: SI = S( I )
282: C2 = ( N-1 ) * T
283: C1 = ( N-2 ) * ( WORK( I ) - T*SI )
284: C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
285: D = C1*C1 - 4*C0*C2
286:
287: IF ( D .LE. 0 ) THEN
288: INFO = -1
289: RETURN
290: END IF
291: SI = -2*C0 / ( C1 + SQRT( D ) )
292:
293: D = SI - S( I )
294: U = ZERO
295: IF ( UP ) THEN
296: DO J = 1, I
297: T = ABS( A( J, I ) )
298: U = U + S( J )*T
299: WORK( J ) = WORK( J ) + D*T
300: END DO
301: DO J = I+1,N
302: T = ABS( A( I, J ) )
303: U = U + S( J )*T
304: WORK( J ) = WORK( J ) + D*T
305: END DO
306: ELSE
307: DO J = 1, I
308: T = ABS( A( I, J ) )
309: U = U + S( J )*T
310: WORK( J ) = WORK( J ) + D*T
311: END DO
312: DO J = I+1,N
313: T = ABS( A( J, I ) )
314: U = U + S( J )*T
315: WORK( J ) = WORK( J ) + D*T
316: END DO
317: END IF
318:
319: AVG = AVG + ( U + WORK( I ) ) * D / N
320: S( I ) = SI
321:
322: END DO
323:
324: END DO
325:
326: 999 CONTINUE
327:
328: SMLNUM = DLAMCH( 'SAFEMIN' )
329: BIGNUM = ONE / SMLNUM
330: SMIN = BIGNUM
331: SMAX = ZERO
332: T = ONE / SQRT(AVG)
333: BASE = DLAMCH( 'B' )
334: U = ONE / LOG( BASE )
335: DO I = 1, N
336: S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
337: SMIN = MIN( SMIN, S( I ) )
338: SMAX = MAX( SMAX, S( I ) )
339: END DO
340: SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
341: *
342: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsyequb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack