SRC/dsyequb.f

   1 *> \brief \b DSYEQUB
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYEQUB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyequb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyequb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyequb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYEQUB( 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 *       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSYEQUB computes row and column scalings intended to equilibrate a
  39 *> symmetric matrix A and reduce its condition number
  40 *> (with respect to the two-norm).  S contains the scale factors,
  41 *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
  42 *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
  43 *> choice of S puts the condition number of B within a factor N of the
  44 *> 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 *>          Specifies whether the details of the factorization are stored
  55 *>          as an upper or lower triangular matrix.
  56 *>          = 'U':  Upper triangular, form is A = U*D*U**T;
  57 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] N
  61 *> \verbatim
  62 *>          N is INTEGER
  63 *>          The order of the matrix A.  N >= 0.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] A
  67 *> \verbatim
  68 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  69 *>          The N-by-N symmetric matrix whose scaling
  70 *>          factors are to be computed.  Only the diagonal elements of A
  71 *>          are referenced.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] LDA
  75 *> \verbatim
  76 *>          LDA is INTEGER
  77 *>          The leading dimension of the array A.  LDA >= max(1,N).
  78 *> \endverbatim
  79 *>
  80 *> \param[out] S
  81 *> \verbatim
  82 *>          S is DOUBLE PRECISION array, dimension (N)
  83 *>          If INFO = 0, S contains the scale factors for A.
  84 *> \endverbatim
  85 *>
  86 *> \param[out] SCOND
  87 *> \verbatim
  88 *>          SCOND is DOUBLE PRECISION
  89 *>          If INFO = 0, S contains the ratio of the smallest S(i) to
  90 *>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
  91 *>          large nor too small, it is not worth scaling by S.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] AMAX
  95 *> \verbatim
  96 *>          AMAX is DOUBLE PRECISION
  97 *>          Absolute value of largest matrix element.  If AMAX is very
  98 *>          close to overflow or very close to underflow, the matrix
  99 *>          should be scaled.
 100 *> \endverbatim
 101 *>
 102 *> \param[out] WORK
 103 *> \verbatim
 104 *>          WORK is DOUBLE PRECISION array, dimension (3*N)
 105 *> \endverbatim
 106 *>
 107 *> \param[out] INFO
 108 *> \verbatim
 109 *>          INFO is INTEGER
 110 *>          = 0:  successful exit
 111 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 112 *>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
 113 *> \endverbatim
 114 *
 115 *  Authors:
 116 *  ========
 117 *
 118 *> \author Univ. of Tennessee
 119 *> \author Univ. of California Berkeley
 120 *> \author Univ. of Colorado Denver
 121 *> \author NAG Ltd.
 122 *
 123 *> \date November 2011
 124 *
 125 *> \ingroup doubleSYcomputational
 126 *
 127 *> \par References:
 128 *  ================
 129 *>
 130 *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
 131 *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
 132 *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
 133 *>  Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf
 134 *>
 135 *  =====================================================================
 136       SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 137 *
 138 *  -- LAPACK computational routine (version 3.4.0) --
 139 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 140 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 141 *     November 2011
 142 *
 143 *     .. Scalar Arguments ..
 144       INTEGER            INFO, LDA, N
 145       DOUBLE PRECISION   AMAX, SCOND
 146       CHARACTER          UPLO
 147 *     ..
 148 *     .. Array Arguments ..
 149       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )
 150 *     ..
 151 *
 152 *  =====================================================================
 153 *
 154 *     .. Parameters ..
 155       DOUBLE PRECISION   ONE, ZERO
 156       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 157       INTEGER            MAX_ITER
 158       PARAMETER          ( MAX_ITER = 100 )
 159 *     ..
 160 *     .. Local Scalars ..
 161       INTEGER            I, J, ITER
 162       DOUBLE PRECISION   AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
 163      $                   SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
 164       LOGICAL            UP
 165 *     ..
 166 *     .. External Functions ..
 167       DOUBLE PRECISION   DLAMCH
 168       LOGICAL            LSAME
 169       EXTERNAL           DLAMCH, LSAME
 170 *     ..
 171 *     .. External Subroutines ..
 172       EXTERNAL           DLASSQ
 173 *     ..
 174 *     .. Intrinsic Functions ..
 175       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
 176 *     ..
 177 *     .. Executable Statements ..
 178 *
 179 *     Test input parameters.
 180 *
 181       INFO = 0
 182       IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
 183         INFO = -1
 184       ELSE IF ( N .LT. 0 ) THEN
 185         INFO = -2
 186       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
 187         INFO = -4
 188       END IF
 189       IF ( INFO .NE. 0 ) THEN
 190         CALL XERBLA( 'DSYEQUB', -INFO )
 191         RETURN
 192       END IF
 193
 194       UP = LSAME( UPLO, 'U' )
 195       AMAX = ZERO
 196 *
 197 *     Quick return if possible.
 198 *
 199       IF ( N .EQ. 0 ) THEN
 200         SCOND = ONE
 201         RETURN
 202       END IF
 203
 204       DO I = 1, N
 205         S( I ) = ZERO
 206       END DO
 207
 208       AMAX = ZERO
 209       IF ( UP ) THEN
 210          DO J = 1, N
 211             DO I = 1, J-1
 212                S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
 213                S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
 214                AMAX = MAX( AMAX, ABS( A(I, J) ) )
 215             END DO
 216             S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
 217             AMAX = MAX( AMAX, ABS( A( J, J ) ) )
 218          END DO
 219       ELSE
 220          DO J = 1, N
 221             S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
 222             AMAX = MAX( AMAX, ABS( A( J, J ) ) )
 223             DO I = J+1, N
 224                S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
 225                S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
 226                AMAX = MAX( AMAX, ABS( A( I, J ) ) )
 227             END DO
 228          END DO
 229       END IF
 230       DO J = 1, N
 231          S( J ) = 1.0D+0 / S( J )
 232       END DO
 233
 234       TOL = ONE / SQRT(2.0D0 * N)
 235
 236       DO ITER = 1, MAX_ITER
 237          SCALE = 0.0D+0
 238          SUMSQ = 0.0D+0
 239 *       BETA = |A|S
 240         DO I = 1, N
 241            WORK(I) = ZERO
 242         END DO
 243         IF ( UP ) THEN
 244            DO J = 1, N
 245               DO I = 1, J-1
 246                  T = ABS( A( I, J ) )
 247                  WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
 248                  WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
 249               END DO
 250               WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
 251            END DO
 252         ELSE
 253            DO J = 1, N
 254               WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
 255               DO I = J+1, N
 256                  T = ABS( A( I, J ) )
 257                  WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
 258                  WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
 259               END DO
 260            END DO
 261         END IF
 262
 263 *       avg = s^T beta / n
 264         AVG = 0.0D+0
 265         DO I = 1, N
 266           AVG = AVG + S( I )*WORK( I )
 267         END DO
 268         AVG = AVG / N
 269
 270         STD = 0.0D+0
 271         DO I = 2*N+1, 3*N
 272            WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG
 273         END DO
 274         CALL DLASSQ( N, WORK( 2*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 = ABS( A( I, I ) )
 281           SI = S( I )
 282           C2 = ( N-1 ) * T
 283           C1 = ( N-2 ) * ( WORK( I ) - T*SI )
 284           C0 = -(T*SI)*SI + 2*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 = ABS( 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 = ABS( 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 = ABS( 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 = ABS( 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 + WORK( I ) ) * D / N
 320           S( I ) = SI
 321
 322         END DO
 323
 324       END DO
 325
 326  999  CONTINUE
 327
 328       SMLNUM = DLAMCH( 'SAFEMIN' )
 329       BIGNUM = ONE / SMLNUM
 330       SMIN = BIGNUM
 331       SMAX = ZERO
 332       T = ONE / SQRT(AVG)
 333       BASE = DLAMCH( 'B' )
 334       U = ONE / LOG( BASE )
 335       DO I = 1, N
 336         S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
 337         SMIN = MIN( SMIN, S( I ) )
 338         SMAX = MAX( SMAX, S( I ) )
 339       END DO
 340       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
 341 *
 342       END