LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ zheequb()

 subroutine zheequb ( 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 )

ZHEEQUB

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

Purpose:
``` ZHEEQUB computes row and column scalings intended to equilibrate a
Hermitian 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 Hermitian 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.```
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 131 of file zheequb.f.

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  DOUBLE PRECISION AMAX, SCOND
140  CHARACTER UPLO
141 * ..
142 * .. Array Arguments ..
143  COMPLEX*16 A( LDA, * ), WORK( * )
144  DOUBLE PRECISION S( * )
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  DOUBLE PRECISION ONE, ZERO
151  parameter( one = 1.0d0, zero = 0.0d0 )
152  INTEGER MAX_ITER
153  parameter( max_iter = 100 )
154 * ..
155 * .. Local Scalars ..
156  INTEGER I, J, ITER
157  DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
158  \$ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
159  LOGICAL UP
160  COMPLEX*16 ZDUM
161 * ..
162 * .. External Functions ..
163  DOUBLE PRECISION DLAMCH
164  LOGICAL LSAME
165  EXTERNAL dlamch, lsame
166 * ..
167 * .. External Subroutines ..
168  EXTERNAL zlassq, xerbla
169 * ..
170 * .. Intrinsic Functions ..
171  INTRINSIC abs, dble, dimag, int, log, max, min, sqrt
172 * ..
173 * .. Statement Functions ..
174  DOUBLE PRECISION CABS1
175 * ..
176 * .. Statement Function Definitions ..
177  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZHEEQUB', -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.0d0 / s( j )
234  END DO
235
236  tol = one / sqrt( 2.0d0 * n )
237
238  DO iter = 1, max_iter
239  scale = 0.0d0
240  sumsq = 0.0d0
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.0d0
265  DO i = 1, n
266  avg = avg + dble( s( i )*work( i ) )
267  END DO
268  avg = avg / n
269
270  std = 0.0d0
271  DO i = n+1, 2*n
272  work( i ) = s( i-n ) * work( i-n ) - avg
273  END DO
274  CALL zlassq( 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 = ( n-2 ) * ( dble( work( i ) ) - t*si )
284  c0 = -(t*si)*si + 2 * dble( 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 + dble( work( i ) ) ) * d / n
320  s( i ) = si
321  END DO
322  END DO
323
324  999 CONTINUE
325
326  smlnum = dlamch( 'SAFEMIN' )
327  bignum = one / smlnum
328  smin = bignum
329  smax = zero
330  t = one / sqrt( avg )
331  base = dlamch( '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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f90:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
Here is the call graph for this function:
Here is the caller graph for this function: