LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 (with respect to the Euclidean norm) and reduce
40 *> its condition number. The scale factors S are computed by the BIN
41 *> algorithm (see references) so that the scaled matrix B with elements
42 *> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
43 *> the smallest possible condition number over all possible diagonal
44 *> scalings.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] UPLO
51 *> \verbatim
52 *> UPLO is CHARACTER*1
53 *> = 'U': Upper triangle of A is stored;
54 *> = 'L': Lower triangle of A is stored.
55 *> \endverbatim
56 *>
57 *> \param[in] N
58 *> \verbatim
59 *> N is INTEGER
60 *> The order of the matrix A. N >= 0.
61 *> \endverbatim
62 *>
63 *> \param[in] A
64 *> \verbatim
65 *> A is DOUBLE PRECISION array, dimension (LDA,N)
66 *> The N-by-N symmetric matrix whose scaling factors are to be
67 *> computed.
68 *> \endverbatim
69 *>
70 *> \param[in] LDA
71 *> \verbatim
72 *> LDA is INTEGER
73 *> The leading dimension of the array A. LDA >= max(1,N).
74 *> \endverbatim
75 *>
76 *> \param[out] S
77 *> \verbatim
78 *> S is DOUBLE PRECISION array, dimension (N)
79 *> If INFO = 0, S contains the scale factors for A.
80 *> \endverbatim
81 *>
82 *> \param[out] SCOND
83 *> \verbatim
84 *> SCOND is DOUBLE PRECISION
85 *> If INFO = 0, S contains the ratio of the smallest S(i) to
86 *> the largest S(i). If SCOND >= 0.1 and AMAX is neither too
87 *> large nor too small, it is not worth scaling by S.
88 *> \endverbatim
89 *>
90 *> \param[out] AMAX
91 *> \verbatim
92 *> AMAX is DOUBLE PRECISION
93 *> Largest absolute value of any matrix element. If AMAX is
94 *> very close to overflow or very close to underflow, the
95 *> matrix should be scaled.
96 *> \endverbatim
97 *>
98 *> \param[out] WORK
99 *> \verbatim
100 *> WORK is DOUBLE PRECISION array, dimension (2*N)
101 *> \endverbatim
102 *>
103 *> \param[out] INFO
104 *> \verbatim
105 *> INFO is INTEGER
106 *> = 0: successful exit
107 *> < 0: if INFO = -i, the i-th argument had an illegal value
108 *> > 0: if INFO = i, the i-th diagonal element is nonpositive.
109 *> \endverbatim
110 *
111 * Authors:
112 * ========
113 *
114 *> \author Univ. of Tennessee
115 *> \author Univ. of California Berkeley
116 *> \author Univ. of Colorado Denver
117 *> \author NAG Ltd.
118 *
119 *> \date November 2017
120 *
121 *> \ingroup doubleSYcomputational
122 *
123 *> \par References:
124 * ================
125 *>
126 *> Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
127 *> Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
128 *> DOI 10.1023/B:NUMA.0000016606.32820.69 \n
129 *> Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
130 *>
131 * =====================================================================
132  SUBROUTINE dsyequb( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
133 *
134 * -- LAPACK computational routine (version 3.8.0) --
135 * -- LAPACK is a software package provided by Univ. of Tennessee, --
136 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137 * November 2017
138 *
139 * .. Scalar Arguments ..
140  INTEGER INFO, LDA, N
141  DOUBLE PRECISION AMAX, SCOND
142  CHARACTER UPLO
143 * ..
144 * .. Array Arguments ..
145  DOUBLE PRECISION A( LDA, * ), S( * ), WORK( * )
146 * ..
147 *
148 * =====================================================================
149 *
150 * .. Parameters ..
151  DOUBLE PRECISION ONE, ZERO
152  parameter( one = 1.0d0, zero = 0.0d0 )
153  INTEGER MAX_ITER
154  parameter( max_iter = 100 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER I, J, ITER
158  DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
159  $ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
160  LOGICAL UP
161 * ..
162 * .. External Functions ..
163  DOUBLE PRECISION DLAMCH
164  LOGICAL LSAME
165  EXTERNAL dlamch, lsame
166 * ..
167 * .. External Subroutines ..
168  EXTERNAL dlassq, xerbla
169 * ..
170 * .. Intrinsic Functions ..
171  INTRINSIC abs, int, log, max, min, sqrt
172 * ..
173 * .. Executable Statements ..
174 *
175 * Test the input parameters.
176 *
177  info = 0
178  IF ( .NOT. ( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) ) THEN
179  info = -1
180  ELSE IF ( n .LT. 0 ) THEN
181  info = -2
182  ELSE IF ( lda .LT. max( 1, n ) ) THEN
183  info = -4
184  END IF
185  IF ( info .NE. 0 ) THEN
186  CALL xerbla( 'DSYEQUB', -info )
187  RETURN
188  END IF
189 
190  up = lsame( uplo, 'U' )
191  amax = zero
192 *
193 * Quick return if possible.
194 *
195  IF ( n .EQ. 0 ) THEN
196  scond = one
197  RETURN
198  END IF
199 
200  DO i = 1, n
201  s( i ) = zero
202  END DO
203 
204  amax = zero
205  IF ( up ) THEN
206  DO j = 1, n
207  DO i = 1, j-1
208  s( i ) = max( s( i ), abs( a( i, j ) ) )
209  s( j ) = max( s( j ), abs( a( i, j ) ) )
210  amax = max( amax, abs( a( i, j ) ) )
211  END DO
212  s( j ) = max( s( j ), abs( a( j, j ) ) )
213  amax = max( amax, abs( a( j, j ) ) )
214  END DO
215  ELSE
216  DO j = 1, n
217  s( j ) = max( s( j ), abs( a( j, j ) ) )
218  amax = max( amax, abs( a( j, j ) ) )
219  DO i = j+1, n
220  s( i ) = max( s( i ), abs( a( i, j ) ) )
221  s( j ) = max( s( j ), abs( a( i, j ) ) )
222  amax = max( amax, abs( a( i, j ) ) )
223  END DO
224  END DO
225  END IF
226  DO j = 1, n
227  s( j ) = 1.0d0 / s( j )
228  END DO
229 
230  tol = one / sqrt( 2.0d0 * n )
231 
232  DO iter = 1, max_iter
233  scale = 0.0d0
234  sumsq = 0.0d0
235 * beta = |A|s
236  DO i = 1, n
237  work( i ) = zero
238  END DO
239  IF ( up ) THEN
240  DO j = 1, n
241  DO i = 1, j-1
242  work( i ) = work( i ) + abs( a( i, j ) ) * s( j )
243  work( j ) = work( j ) + abs( a( i, j ) ) * s( i )
244  END DO
245  work( j ) = work( j ) + abs( a( j, j ) ) * s( j )
246  END DO
247  ELSE
248  DO j = 1, n
249  work( j ) = work( j ) + abs( a( j, j ) ) * s( j )
250  DO i = j+1, n
251  work( i ) = work( i ) + abs( a( i, j ) ) * s( j )
252  work( j ) = work( j ) + abs( a( i, j ) ) * s( i )
253  END DO
254  END DO
255  END IF
256 
257 * avg = s^T beta / n
258  avg = 0.0d0
259  DO i = 1, n
260  avg = avg + s( i )*work( i )
261  END DO
262  avg = avg / n
263 
264  std = 0.0d0
265  DO i = n+1, 2*n
266  work( i ) = s( i-n ) * work( i-n ) - avg
267  END DO
268  CALL dlassq( n, work( n+1 ), 1, scale, sumsq )
269  std = scale * sqrt( sumsq / n )
270 
271  IF ( std .LT. tol * avg ) GOTO 999
272 
273  DO i = 1, n
274  t = abs( a( i, i ) )
275  si = s( i )
276  c2 = ( n-1 ) * t
277  c1 = ( n-2 ) * ( work( i ) - t*si )
278  c0 = -(t*si)*si + 2*work( i )*si - n*avg
279  d = c1*c1 - 4*c0*c2
280 
281  IF ( d .LE. 0 ) THEN
282  info = -1
283  RETURN
284  END IF
285  si = -2*c0 / ( c1 + sqrt( d ) )
286 
287  d = si - s( i )
288  u = zero
289  IF ( up ) THEN
290  DO j = 1, i
291  t = abs( a( j, i ) )
292  u = u + s( j )*t
293  work( j ) = work( j ) + d*t
294  END DO
295  DO j = i+1,n
296  t = abs( a( i, j ) )
297  u = u + s( j )*t
298  work( j ) = work( j ) + d*t
299  END DO
300  ELSE
301  DO j = 1, i
302  t = abs( a( i, j ) )
303  u = u + s( j )*t
304  work( j ) = work( j ) + d*t
305  END DO
306  DO j = i+1,n
307  t = abs( a( j, i ) )
308  u = u + s( j )*t
309  work( j ) = work( j ) + d*t
310  END DO
311  END IF
312 
313  avg = avg + ( u + work( i ) ) * d / n
314  s( i ) = si
315  END DO
316  END DO
317 
318  999 CONTINUE
319 
320  smlnum = dlamch( 'SAFEMIN' )
321  bignum = one / smlnum
322  smin = bignum
323  smax = zero
324  t = one / sqrt( avg )
325  base = dlamch( 'B' )
326  u = one / log( base )
327  DO i = 1, n
328  s( i ) = base ** int( u * log( s( i ) * t ) )
329  smin = min( smin, s( i ) )
330  smax = max( smax, s( i ) )
331  END DO
332  scond = max( smin, smlnum ) / min( smax, bignum )
333 *
334  END
dsyequb
subroutine dsyequb(UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO)
DSYEQUB
Definition: dsyequb.f:133
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlassq
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105