LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssyequb.f
Go to the documentation of this file.
1*> \brief \b SSYEQUB
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSYEQUB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyequb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyequb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyequb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, LDA, N
25* REAL AMAX, SCOND
26* CHARACTER UPLO
27* ..
28* .. Array Arguments ..
29* REAL A( LDA, * ), S( * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SSYEQUB 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 REAL 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 REAL 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 REAL
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 REAL
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 REAL 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*> \ingroup realSYcomputational
120*
121*> \par References:
122* ================
123*>
124*> Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
125*> Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
126*> DOI 10.1023/B:NUMA.0000016606.32820.69 \n
127*> Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
128*>
129* =====================================================================
130 SUBROUTINE ssyequb( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
131*
132* -- LAPACK computational routine --
133* -- LAPACK is a software package provided by Univ. of Tennessee, --
134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*
136* .. Scalar Arguments ..
137 INTEGER INFO, LDA, N
138 REAL AMAX, SCOND
139 CHARACTER UPLO
140* ..
141* .. Array Arguments ..
142 REAL A( LDA, * ), S( * ), WORK( * )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 REAL ONE, ZERO
149 parameter( one = 1.0e0, zero = 0.0e0 )
150 INTEGER MAX_ITER
151 parameter( max_iter = 100 )
152* ..
153* .. Local Scalars ..
154 INTEGER I, J, ITER
155 REAL AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
156 $ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
157 LOGICAL UP
158* ..
159* .. External Functions ..
160 REAL SLAMCH
161 LOGICAL LSAME
162 EXTERNAL lsame, slamch
163* ..
164* .. External Subroutines ..
165 EXTERNAL slassq, xerbla
166* ..
167* .. Intrinsic Functions ..
168 INTRINSIC abs, int, log, max, min, sqrt
169* ..
170* .. Executable Statements ..
171*
172* Test the input parameters.
173*
174 info = 0
175 IF ( .NOT. ( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) ) THEN
176 info = -1
177 ELSE IF ( n .LT. 0 ) THEN
178 info = -2
179 ELSE IF ( lda .LT. max( 1, n ) ) THEN
180 info = -4
181 END IF
182 IF ( info .NE. 0 ) THEN
183 CALL xerbla( 'SSYEQUB', -info )
184 RETURN
185 END IF
186
187 up = lsame( uplo, 'U' )
188 amax = zero
189*
190* Quick return if possible.
191*
192 IF ( n .EQ. 0 ) THEN
193 scond = one
194 RETURN
195 END IF
196
197 DO i = 1, n
198 s( i ) = zero
199 END DO
200
201 amax = zero
202 IF ( up ) THEN
203 DO j = 1, n
204 DO i = 1, j-1
205 s( i ) = max( s( i ), abs( a( i, j ) ) )
206 s( j ) = max( s( j ), abs( a( i, j ) ) )
207 amax = max( amax, abs( a( i, j ) ) )
208 END DO
209 s( j ) = max( s( j ), abs( a( j, j ) ) )
210 amax = max( amax, abs( a( j, j ) ) )
211 END DO
212 ELSE
213 DO j = 1, n
214 s( j ) = max( s( j ), abs( a( j, j ) ) )
215 amax = max( amax, abs( a( j, j ) ) )
216 DO i = j+1, n
217 s( i ) = max( s( i ), abs( a( i, j ) ) )
218 s( j ) = max( s( j ), abs( a( i, j ) ) )
219 amax = max( amax, abs( a( i, j ) ) )
220 END DO
221 END DO
222 END IF
223 DO j = 1, n
224 s( j ) = 1.0e0 / s( j )
225 END DO
226
227 tol = one / sqrt( 2.0e0 * n )
228
229 DO iter = 1, max_iter
230 scale = 0.0e0
231 sumsq = 0.0e0
232* beta = |A|s
233 DO i = 1, n
234 work( i ) = zero
235 END DO
236 IF ( up ) THEN
237 DO j = 1, n
238 DO i = 1, j-1
239 work( i ) = work( i ) + abs( a( i, j ) ) * s( j )
240 work( j ) = work( j ) + abs( a( i, j ) ) * s( i )
241 END DO
242 work( j ) = work( j ) + abs( a( j, j ) ) * s( j )
243 END DO
244 ELSE
245 DO j = 1, n
246 work( j ) = work( j ) + abs( a( j, j ) ) * s( j )
247 DO i = j+1, n
248 work( i ) = work( i ) + abs( a( i, j ) ) * s( j )
249 work( j ) = work( j ) + abs( a( i, j ) ) * s( i )
250 END DO
251 END DO
252 END IF
253
254* avg = s^T beta / n
255 avg = 0.0e0
256 DO i = 1, n
257 avg = avg + s( i )*work( i )
258 END DO
259 avg = avg / n
260
261 std = 0.0e0
262 DO i = n+1, 2*n
263 work( i ) = s( i-n ) * work( i-n ) - avg
264 END DO
265 CALL slassq( n, work( n+1 ), 1, scale, sumsq )
266 std = scale * sqrt( sumsq / n )
267
268 IF ( std .LT. tol * avg ) GOTO 999
269
270 DO i = 1, n
271 t = abs( a( i, i ) )
272 si = s( i )
273 c2 = ( n-1 ) * t
274 c1 = ( n-2 ) * ( work( i ) - t*si )
275 c0 = -(t*si)*si + 2*work( i )*si - n*avg
276 d = c1*c1 - 4*c0*c2
277
278 IF ( d .LE. 0 ) THEN
279 info = -1
280 RETURN
281 END IF
282 si = -2*c0 / ( c1 + sqrt( d ) )
283
284 d = si - s( i )
285 u = zero
286 IF ( up ) THEN
287 DO j = 1, i
288 t = abs( a( j, i ) )
289 u = u + s( j )*t
290 work( j ) = work( j ) + d*t
291 END DO
292 DO j = i+1,n
293 t = abs( a( i, j ) )
294 u = u + s( j )*t
295 work( j ) = work( j ) + d*t
296 END DO
297 ELSE
298 DO j = 1, i
299 t = abs( a( i, j ) )
300 u = u + s( j )*t
301 work( j ) = work( j ) + d*t
302 END DO
303 DO j = i+1,n
304 t = abs( a( j, i ) )
305 u = u + s( j )*t
306 work( j ) = work( j ) + d*t
307 END DO
308 END IF
309
310 avg = avg + ( u + work( i ) ) * d / n
311 s( i ) = si
312 END DO
313 END DO
314
315 999 CONTINUE
316
317 smlnum = slamch( 'SAFEMIN' )
318 bignum = one / smlnum
319 smin = bignum
320 smax = zero
321 t = one / sqrt( avg )
322 base = slamch( 'B' )
323 u = one / log( base )
324 DO i = 1, n
325 s( i ) = base ** int( u * log( s( i ) * t ) )
326 smin = min( smin, s( i ) )
327 smax = max( smax, s( i ) )
328 END DO
329 scond = max( smin, smlnum ) / min( smax, bignum )
330*
331 END
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f90:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssyequb(UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO)
SSYEQUB
Definition: ssyequb.f:131