LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
csyequb.f
Go to the documentation of this file.
1*> \brief \b CSYEQUB
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CSYEQUB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csyequb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csyequb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csyequb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSYEQUB( 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* COMPLEX A( LDA, * ), WORK( * )
30* REAL S( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CSYEQUB computes row and column scalings intended to equilibrate a
40*> symmetric matrix A (with respect to the Euclidean norm) and reduce
41*> its condition number. The scale factors S are computed by the BIN
42*> algorithm (see references) so that the scaled matrix B with elements
43*> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
44*> the 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*> = 'U': Upper triangle of A is stored;
55*> = 'L': Lower triangle of A is stored.
56*> \endverbatim
57*>
58*> \param[in] N
59*> \verbatim
60*> N is INTEGER
61*> The order of the matrix A. N >= 0.
62*> \endverbatim
63*>
64*> \param[in] A
65*> \verbatim
66*> A is COMPLEX array, dimension (LDA,N)
67*> The N-by-N symmetric matrix whose scaling factors are to be
68*> computed.
69*> \endverbatim
70*>
71*> \param[in] LDA
72*> \verbatim
73*> LDA is INTEGER
74*> The leading dimension of the array A. LDA >= max(1,N).
75*> \endverbatim
76*>
77*> \param[out] S
78*> \verbatim
79*> S is REAL array, dimension (N)
80*> If INFO = 0, S contains the scale factors for A.
81*> \endverbatim
82*>
83*> \param[out] SCOND
84*> \verbatim
85*> SCOND is REAL
86*> If INFO = 0, S contains the ratio of the smallest S(i) to
87*> the largest S(i). If SCOND >= 0.1 and AMAX is neither too
88*> large nor too small, it is not worth scaling by S.
89*> \endverbatim
90*>
91*> \param[out] AMAX
92*> \verbatim
93*> AMAX is REAL
94*> Largest absolute value of any matrix element. If AMAX is
95*> very close to overflow or very close to underflow, the
96*> matrix should be scaled.
97*> \endverbatim
98*>
99*> \param[out] WORK
100*> \verbatim
101*> WORK is COMPLEX array, dimension (2*N)
102*> \endverbatim
103*>
104*> \param[out] INFO
105*> \verbatim
106*> INFO is INTEGER
107*> = 0: successful exit
108*> < 0: if INFO = -i, the i-th argument had an illegal value
109*> > 0: if INFO = i, the i-th diagonal element is nonpositive.
110*> \endverbatim
111*
112* Authors:
113* ========
114*
115*> \author Univ. of Tennessee
116*> \author Univ. of California Berkeley
117*> \author Univ. of Colorado Denver
118*> \author NAG Ltd.
119*
120*> \ingroup heequb
121*
122*> \par References:
123* ================
124*>
125*> Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
126*> Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
127*> DOI 10.1023/B:NUMA.0000016606.32820.69 \n
128*> Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
129*>
130* =====================================================================
131 SUBROUTINE csyequb( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
132*
133* -- LAPACK computational routine --
134* -- LAPACK is a software package provided by Univ. of Tennessee, --
135* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137* .. Scalar Arguments ..
138 INTEGER INFO, LDA, N
139 REAL AMAX, SCOND
140 CHARACTER UPLO
141* ..
142* .. Array Arguments ..
143 COMPLEX A( LDA, * ), WORK( * )
144 REAL S( * )
145* ..
146*
147* =====================================================================
148*
149* .. Parameters ..
150 REAL ONE, ZERO
151 parameter( one = 1.0e0, zero = 0.0e0 )
152 INTEGER MAX_ITER
153 parameter( max_iter = 100 )
154* ..
155* .. Local Scalars ..
156 INTEGER I, J, ITER
157 REAL AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
158 $ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
159 LOGICAL UP
160 COMPLEX ZDUM
161* ..
162* .. External Functions ..
163 REAL SLAMCH
164 LOGICAL LSAME
165 EXTERNAL lsame, slamch
166* ..
167* .. External Subroutines ..
168 EXTERNAL classq, xerbla
169* ..
170* .. Intrinsic Functions ..
171 INTRINSIC abs, aimag, int, log, max, min, real, sqrt
172* ..
173* .. Statement Functions ..
174 REAL CABS1
175* ..
176* .. Statement Function Definitions ..
177 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
178* ..
179* .. Executable Statements ..
180*
181* Test the input parameters.
182*
183 info = 0
184 IF ( .NOT. ( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) ) THEN
185 info = -1
186 ELSE IF ( n .LT. 0 ) THEN
187 info = -2
188 ELSE IF ( lda .LT. max( 1, n ) ) THEN
189 info = -4
190 END IF
191 IF ( info .NE. 0 ) THEN
192 CALL xerbla( 'CSYEQUB', -info )
193 RETURN
194 END IF
195
196 up = lsame( uplo, 'U' )
197 amax = zero
198*
199* Quick return if possible.
200*
201 IF ( n .EQ. 0 ) THEN
202 scond = one
203 RETURN
204 END IF
205
206 DO i = 1, n
207 s( i ) = zero
208 END DO
209
210 amax = zero
211 IF ( up ) THEN
212 DO j = 1, n
213 DO i = 1, j-1
214 s( i ) = max( s( i ), cabs1( a( i, j ) ) )
215 s( j ) = max( s( j ), cabs1( a( i, j ) ) )
216 amax = max( amax, cabs1( a( i, j ) ) )
217 END DO
218 s( j ) = max( s( j ), cabs1( a( j, j ) ) )
219 amax = max( amax, cabs1( a( j, j ) ) )
220 END DO
221 ELSE
222 DO j = 1, n
223 s( j ) = max( s( j ), cabs1( a( j, j ) ) )
224 amax = max( amax, cabs1( a( j, j ) ) )
225 DO i = j+1, n
226 s( i ) = max( s( i ), cabs1( a( i, j ) ) )
227 s( j ) = max( s( j ), cabs1( a( i, j ) ) )
228 amax = max( amax, cabs1( a( i, j ) ) )
229 END DO
230 END DO
231 END IF
232 DO j = 1, n
233 s( j ) = 1.0 / s( j )
234 END DO
235
236 tol = one / sqrt( 2.0e0 * n )
237
238 DO iter = 1, max_iter
239 scale = 0.0e0
240 sumsq = 0.0e0
241* beta = |A|s
242 DO i = 1, n
243 work( i ) = zero
244 END DO
245 IF ( up ) THEN
246 DO j = 1, n
247 DO i = 1, j-1
248 work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
249 work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
250 END DO
251 work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
252 END DO
253 ELSE
254 DO j = 1, n
255 work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
256 DO i = j+1, n
257 work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
258 work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
259 END DO
260 END DO
261 END IF
262
263* avg = s^T beta / n
264 avg = 0.0e0
265 DO i = 1, n
266 avg = avg + real( s( i )*work( i ) )
267 END DO
268 avg = avg / n
269
270 std = 0.0e0
271 DO i = n+1, 2*n
272 work( i ) = s( i-n ) * work( i-n ) - avg
273 END DO
274 CALL classq( n, work( 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 = cabs1( a( i, i ) )
281 si = s( i )
282 c2 = ( n-1 ) * t
283 c1 = real( n-2 ) * ( real( work( i ) ) - t*si )
284 c0 = -(t*si)*si + 2 * real( 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 = cabs1( 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 = cabs1( 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 = cabs1( 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 = cabs1( 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 + real( work( i ) ) ) * d / n
320 s( i ) = si
321 END DO
322 END DO
323
324 999 CONTINUE
325
326 smlnum = slamch( 'SAFEMIN' )
327 bignum = one / smlnum
328 smin = bignum
329 smax = zero
330 t = one / sqrt( avg )
331 base = slamch( 'B' )
332 u = one / log( base )
333 DO i = 1, n
334 s( i ) = base ** int( u * log( s( i ) * t ) )
335 smin = min( smin, s( i ) )
336 smax = max( smax, s( i ) )
337 END DO
338 scond = max( smin, smlnum ) / min( smax, bignum )
339*
340 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine csyequb(uplo, n, a, lda, s, scond, amax, work, info)
CSYEQUB
Definition csyequb.f:132
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:124