SRC/zsyequb.f

   1 *> \brief \b ZSYEQUB
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZSYEQUB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsyequb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsyequb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsyequb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZSYEQUB( 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 *> ZSYEQUB 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*16 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 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 *> \date November 2011
 121 *
 122 *> \ingroup complex16SYcomputational
 123 *
 124 *> \par References:
 125 *  ================
 126 *>
 127 *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
 128 *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
 129 *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
 130 *>  Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
 131 *>
 132 *  =====================================================================
 133       SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 134 *
 135 *  -- LAPACK computational routine (version 3.4.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 2011
 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
 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 *
 343       END