LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ zsyequb()

 subroutine zsyequb ( character UPLO, integer N, complex*16, dimension( lda, * ) A, integer LDA, double precision, dimension( * ) S, double precision SCOND, double precision AMAX, complex*16, dimension( * ) WORK, integer INFO )

ZSYEQUB

Download ZSYEQUB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
``` ZSYEQUB computes row and column scalings intended to equilibrate a
symmetric matrix A (with respect to the Euclidean norm) and reduce
its condition number. The scale factors S are computed by the BIN
algorithm (see references) so that the scaled matrix B with elements
B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
the smallest possible condition number over all possible diagonal
scalings.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] A ``` A is COMPLEX*16 array, dimension (LDA,N) The N-by-N symmetric matrix whose scaling factors are to be computed.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] S ``` S is DOUBLE PRECISION array, dimension (N) If INFO = 0, S contains the scale factors for A.``` [out] SCOND ``` SCOND is DOUBLE PRECISION If INFO = 0, S contains the ratio of the smallest S(i) to the largest S(i). If SCOND >= 0.1 and AMAX is neither too large nor too small, it is not worth scaling by S.``` [out] AMAX ``` AMAX is DOUBLE PRECISION Largest absolute value of any matrix element. If AMAX is very close to overflow or very close to underflow, the matrix should be scaled.``` [out] WORK ` WORK is COMPLEX*16 array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the i-th diagonal element is nonpositive.```
Date
November 2017
References:
Livne, O.E. and Golub, G.H., "Scaling by Binormalization",
Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004.
DOI 10.1023/B:NUMA.0000016606.32820.69
Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679

Definition at line 134 of file zsyequb.f.

134 *
135 * -- LAPACK computational routine (version 3.8.0) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * November 2017
139 *
140 * .. Scalar Arguments ..
141  INTEGER info, lda, n
142  DOUBLE PRECISION amax, scond
143  CHARACTER uplo
144 * ..
145 * .. Array Arguments ..
146  COMPLEX*16 a( lda, * ), work( * )
147  DOUBLE PRECISION s( * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  DOUBLE PRECISION one, zero
154  parameter( one = 1.0d0, zero = 0.0d0 )
155  INTEGER max_iter
156  parameter( max_iter = 100 )
157 * ..
158 * .. Local Scalars ..
159  INTEGER i, j, iter
160  DOUBLE PRECISION avg, std, tol, c0, c1, c2, t, u, si, d, base,
161  \$ smin, smax, smlnum, bignum, scale, sumsq
162  LOGICAL up
163  COMPLEX*16 zdum
164 * ..
165 * .. External Functions ..
166  DOUBLE PRECISION dlamch
167  LOGICAL lsame
168  EXTERNAL dlamch, lsame
169 * ..
170 * .. External Subroutines ..
171  EXTERNAL zlassq, xerbla
172 * ..
173 * .. Intrinsic Functions ..
174  INTRINSIC abs, dble, dimag, int, log, max, min, sqrt
175 * ..
176 * .. Statement Functions ..
177  DOUBLE PRECISION cabs1
178 * ..
179 * .. Statement Function Definitions ..
180  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
181 * ..
182 * .. Executable Statements ..
183 *
184 * Test the input parameters.
185 *
186  info = 0
187  IF ( .NOT. ( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) ) THEN
188  info = -1
189  ELSE IF ( n .LT. 0 ) THEN
190  info = -2
191  ELSE IF ( lda .LT. max( 1, n ) ) THEN
192  info = -4
193  END IF
194  IF ( info .NE. 0 ) THEN
195  CALL xerbla( 'ZSYEQUB', -info )
196  RETURN
197  END IF
198
199  up = lsame( uplo, 'U' )
200  amax = zero
201 *
202 * Quick return if possible.
203 *
204  IF ( n .EQ. 0 ) THEN
205  scond = one
206  RETURN
207  END IF
208
209  DO i = 1, n
210  s( i ) = zero
211  END DO
212
213  amax = zero
214  IF ( up ) THEN
215  DO j = 1, n
216  DO i = 1, j-1
217  s( i ) = max( s( i ), cabs1( a( i, j ) ) )
218  s( j ) = max( s( j ), cabs1( a( i, j ) ) )
219  amax = max( amax, cabs1( a( i, j ) ) )
220  END DO
221  s( j ) = max( s( j ), cabs1( a( j, j ) ) )
222  amax = max( amax, cabs1( a( j, j ) ) )
223  END DO
224  ELSE
225  DO j = 1, n
226  s( j ) = max( s( j ), cabs1( a( j, j ) ) )
227  amax = max( amax, cabs1( a( j, j ) ) )
228  DO i = j+1, n
229  s( i ) = max( s( i ), cabs1( a( i, j ) ) )
230  s( j ) = max( s( j ), cabs1( a( i, j ) ) )
231  amax = max( amax, cabs1( a( i, j ) ) )
232  END DO
233  END DO
234  END IF
235  DO j = 1, n
236  s( j ) = 1.0d0 / s( j )
237  END DO
238
239  tol = one / sqrt( 2.0d0 * n )
240
241  DO iter = 1, max_iter
242  scale = 0.0d0
243  sumsq = 0.0d0
244 * beta = |A|s
245  DO i = 1, n
246  work( i ) = zero
247  END DO
248  IF ( up ) THEN
249  DO j = 1, n
250  DO i = 1, j-1
251  work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
252  work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
253  END DO
254  work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
255  END DO
256  ELSE
257  DO j = 1, n
258  work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
259  DO i = j+1, n
260  work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
261  work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
262  END DO
263  END DO
264  END IF
265
266 * avg = s^T beta / n
267  avg = 0.0d0
268  DO i = 1, n
269  avg = avg + s( i )*work( i )
270  END DO
271  avg = avg / n
272
273  std = 0.0d0
274  DO i = n+1, 2*n
275  work( i ) = s( i-n ) * work( i-n ) - avg
276  END DO
277  CALL zlassq( n, work( n+1 ), 1, scale, sumsq )
278  std = scale * sqrt( sumsq / n )
279
280  IF ( std .LT. tol * avg ) GOTO 999
281
282  DO i = 1, n
283  t = cabs1( a( i, i ) )
284  si = s( i )
285  c2 = ( n-1 ) * t
286  c1 = ( n-2 ) * ( work( i ) - t*si )
287  c0 = -(t*si)*si + 2*work( i )*si - n*avg
288  d = c1*c1 - 4*c0*c2
289
290  IF ( d .LE. 0 ) THEN
291  info = -1
292  RETURN
293  END IF
294  si = -2*c0 / ( c1 + sqrt( d ) )
295
296  d = si - s( i )
297  u = zero
298  IF ( up ) THEN
299  DO j = 1, i
300  t = cabs1( a( j, i ) )
301  u = u + s( j )*t
302  work( j ) = work( j ) + d*t
303  END DO
304  DO j = i+1,n
305  t = cabs1( a( i, j ) )
306  u = u + s( j )*t
307  work( j ) = work( j ) + d*t
308  END DO
309  ELSE
310  DO j = 1, i
311  t = cabs1( a( i, j ) )
312  u = u + s( j )*t
313  work( j ) = work( j ) + d*t
314  END DO
315  DO j = i+1,n
316  t = cabs1( a( j, i ) )
317  u = u + s( j )*t
318  work( j ) = work( j ) + d*t
319  END DO
320  END IF
321
322  avg = avg + ( u + work( i ) ) * d / n
323  s( i ) = si
324  END DO
325  END DO
326
327  999 CONTINUE
328
329  smlnum = dlamch( 'SAFEMIN' )
330  bignum = one / smlnum
331  smin = bignum
332  smax = zero
333  t = one / sqrt( avg )
334  base = dlamch( 'B' )
335  u = one / log( base )
336  DO i = 1, n
337  s( i ) = base ** int( u * log( s( i ) * t ) )
338  smin = min( smin, s( i ) )
339  smax = max( smax, s( i ) )
340  END DO
341  scond = max( smin, smlnum ) / min( smax, bignum )
342 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
Here is the call graph for this function:
Here is the caller graph for this function: