1*> \brief \b ZHEEQUB
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHEEQUB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheequb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheequb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheequb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZHEEQUB( 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*       COMPLEX*16         A( LDA, * ), WORK( * )
30*       DOUBLE PRECISION   S( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> ZHEEQUB computes row and column scalings intended to equilibrate a
40*> Hermitian 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*16 array, dimension (LDA,N)
67*>          The N-by-N Hermitian 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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*16 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 complex16HEcomputational
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 ZHEEQUB( 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      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*
340      END
341