Annotation of rpl/lapack/lapack/dsyequb.f, revision 1.5
1.5 ! bertrand 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: * =====================================================================
1.1 bertrand 136: SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
137: *
1.5 ! bertrand 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
1.1 bertrand 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