LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dsyequb.f
Go to the documentation of this file.
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