SRC/cheequb.f

   1 *> \brief \b CHEEQUB
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CHEEQUB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheequb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheequb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheequb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHEEQUB( 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 *> CHEEQUB computes row and column scalings intended to equilibrate a
  40 *> Hermitian matrix A and reduce its condition number
  41 *> (with respect to the two-norm).  S contains the scale factors,
  42 *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
  43 *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
  44 *> choice of S puts the condition number of B within a factor N of the
  45 *> smallest possible condition number over all possible diagonal
  46 *> scalings.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] UPLO
  53 *> \verbatim
  54 *>          UPLO is CHARACTER*1
  55 *>          = 'U':  Upper triangles of A and B are stored;
  56 *>          = 'L':  Lower triangles of A and B are stored.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] N
  60 *> \verbatim
  61 *>          N is INTEGER
  62 *>          The order of the matrix A.  N >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] A
  66 *> \verbatim
  67 *>          A is COMPLEX array, dimension (LDA,N)
  68 *>          The N-by-N Hermitian matrix whose scaling
  69 *>          factors are to be computed.  Only the diagonal elements of A
  70 *>          are referenced.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] LDA
  74 *> \verbatim
  75 *>          LDA is INTEGER
  76 *>          The leading dimension of the array A.  LDA >= max(1,N).
  77 *> \endverbatim
  78 *>
  79 *> \param[out] S
  80 *> \verbatim
  81 *>          S is REAL array, dimension (N)
  82 *>          If INFO = 0, S contains the scale factors for A.
  83 *> \endverbatim
  84 *>
  85 *> \param[out] SCOND
  86 *> \verbatim
  87 *>          SCOND is REAL
  88 *>          If INFO = 0, S contains the ratio of the smallest S(i) to
  89 *>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
  90 *>          large nor too small, it is not worth scaling by S.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] AMAX
  94 *> \verbatim
  95 *>          AMAX is REAL
  96 *>          Absolute value of largest matrix element.  If AMAX is very
  97 *>          close to overflow or very close to underflow, the matrix
  98 *>          should be scaled.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] WORK
 102 *> \verbatim
 103 *>          WORK is COMPLEX array, dimension (3*N)
 104 *> \endverbatim
 105 *>
 106 *> \param[out] INFO
 107 *> \verbatim
 108 *>          INFO is INTEGER
 109 *>          = 0:  successful exit
 110 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 111 *>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
 112 *> \endverbatim
 113 *
 114 *  Authors:
 115 *  ========
 116 *
 117 *> \author Univ. of Tennessee
 118 *> \author Univ. of California Berkeley
 119 *> \author Univ. of Colorado Denver
 120 *> \author NAG Ltd.
 121 *
 122 *> \date April 2012
 123 *
 124 *> \ingroup complexHEcomputational
 125 *
 126 *  =====================================================================
 127       SUBROUTINE CHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 128 *
 129 *  -- LAPACK computational routine (version 3.4.1) --
 130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 132 *     April 2012
 133 *
 134 *     .. Scalar Arguments ..
 135       INTEGER            INFO, LDA, N
 136       REAL               AMAX, SCOND
 137       CHARACTER          UPLO
 138 *     ..
 139 *     .. Array Arguments ..
 140       COMPLEX            A( LDA, * ), WORK( * )
 141       REAL               S( * )
 142 *     ..
 143 *
 144 *  =====================================================================
 145 *
 146 *     .. Parameters ..
 147       REAL               ONE, ZERO
 148       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 149       INTEGER            MAX_ITER
 150       PARAMETER          ( MAX_ITER = 100 )
 151 *     ..
 152 *     .. Local Scalars ..
 153       INTEGER            I, J, ITER
 154       REAL               AVG, STD, TOL, C0, C1, C2, T, U, SI, D,
 155      $                   BASE, SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
 156       LOGICAL            UP
 157       COMPLEX            ZDUM
 158 *     ..
 159 *     .. External Functions ..
 160       REAL               SLAMCH
 161       LOGICAL            LSAME
 162       EXTERNAL           LSAME, SLAMCH
 163 *     ..
 164 *     .. External Subroutines ..
 165       EXTERNAL           CLASSQ
 166 *     ..
 167 *     .. Intrinsic Functions ..
 168       INTRINSIC          ABS, AIMAG, INT, LOG, MAX, MIN, REAL, SQRT
 169 *     ..
 170 *     .. Statement Functions ..
 171       REAL               CABS1
 172 *     ..
 173 *     .. Statement Function Definitions ..
 174       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
 175 *
 176 *     Test input parameters.
 177 *
 178       INFO = 0
 179       IF (.NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
 180         INFO = -1
 181       ELSE IF ( N .LT. 0 ) THEN
 182         INFO = -2
 183       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
 184         INFO = -4
 185       END IF
 186       IF ( INFO .NE. 0 ) THEN
 187         CALL XERBLA( 'CHEEQUB', -INFO )
 188         RETURN
 189       END IF
 190
 191       UP = LSAME( UPLO, 'U' )
 192       AMAX = ZERO
 193 *
 194 *     Quick return if possible.
 195 *
 196       IF ( N .EQ. 0 ) THEN
 197         SCOND = ONE
 198         RETURN
 199       END IF
 200
 201       DO I = 1, N
 202         S( I ) = ZERO
 203       END DO
 204
 205       AMAX = ZERO
 206       IF ( UP ) THEN
 207          DO J = 1, N
 208             DO I = 1, J-1
 209                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
 210                S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )
 211                AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
 212             END DO
 213             S( J ) = MAX( S( J ), CABS1( A( J, J ) ) )
 214             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
 215          END DO
 216       ELSE
 217          DO J = 1, N
 218             S( J ) = MAX( S( J ), CABS1( A( J, J ) ) )
 219             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
 220             DO I = J+1, N
 221                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
 222                S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )
 223                AMAX = MAX( AMAX, CABS1( A(I, J ) ) )
 224             END DO
 225          END DO
 226       END IF
 227       DO J = 1, N
 228          S( J ) = 1.0 / S( J )
 229       END DO
 230
 231       TOL = ONE / SQRT( 2.0E0 * N )
 232
 233       DO ITER = 1, MAX_ITER
 234          SCALE = 0.0
 235          SUMSQ = 0.0
 236 *       beta = |A|s
 237         DO I = 1, N
 238            WORK( I ) = ZERO
 239         END DO
 240         IF ( UP ) THEN
 241            DO J = 1, N
 242               DO I = 1, J-1
 243                  T = CABS1( A( I, J ) )
 244                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
 245                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
 246               END DO
 247               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
 248            END DO
 249         ELSE
 250            DO J = 1, N
 251               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
 252               DO I = J+1, N
 253                  T = CABS1( A( I, J ) )
 254                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
 255                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
 256               END DO
 257            END DO
 258         END IF
 259
 260 *       avg = s^T beta / n
 261         AVG = 0.0
 262         DO I = 1, N
 263           AVG = AVG + S( I )*WORK( I )
 264         END DO
 265         AVG = AVG / N
 266
 267         STD = 0.0
 268         DO I = 2*N+1, 3*N
 269            WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG
 270         END DO
 271         CALL CLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ )
 272         STD = SCALE * SQRT( SUMSQ / N )
 273
 274         IF ( STD .LT. TOL * AVG ) GOTO 999
 275
 276         DO I = 1, N
 277           T = CABS1( A( I, I ) )
 278           SI = S( I )
 279           C2 = ( N-1 ) * T
 280           C1 = ( N-2 ) * ( WORK( I ) - T*SI )
 281           C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
 282
 283           D = C1*C1 - 4*C0*C2
 284           IF ( D .LE. 0 ) THEN
 285             INFO = -1
 286             RETURN
 287           END IF
 288           SI = -2*C0 / ( C1 + SQRT( D ) )
 289
 290           D = SI - S(I)
 291           U = ZERO
 292           IF ( UP ) THEN
 293             DO J = 1, I
 294               T = CABS1( A( J, I ) )
 295               U = U + S( J )*T
 296               WORK( J ) = WORK( J ) + D*T
 297             END DO
 298             DO J = I+1,N
 299               T = CABS1( A( I, J ) )
 300               U = U + S( J )*T
 301               WORK( J ) = WORK( J ) + D*T
 302             END DO
 303           ELSE
 304             DO J = 1, I
 305               T = CABS1( A( I, J ) )
 306               U = U + S( J )*T
 307               WORK( J ) = WORK( J ) + D*T
 308             END DO
 309             DO J = I+1,N
 310               T = CABS1( A( J, I ) )
 311               U = U + S( J )*T
 312               WORK( J ) = WORK( J ) + D*T
 313             END DO
 314           END IF
 315           AVG = AVG + ( U + WORK( I ) ) * D / N
 316           S( I ) = SI
 317         END DO
 318
 319       END DO
 320
 321  999  CONTINUE
 322
 323       SMLNUM = SLAMCH( 'SAFEMIN' )
 324       BIGNUM = ONE / SMLNUM
 325       SMIN = BIGNUM
 326       SMAX = ZERO
 327       T = ONE / SQRT( AVG )
 328       BASE = SLAMCH( 'B' )
 329       U = ONE / LOG( BASE )
 330       DO I = 1, N
 331         S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
 332         SMIN = MIN( SMIN, S( I ) )
 333         SMAX = MAX( SMAX, S( I ) )
 334       END DO
 335       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
 336
 337       END