SRC/dla_syrpvgrw.f

   1 *> \brief \b DLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLA_SYRPVGRW + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
  22 *                                               LDAF, IPIV, WORK )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER*1        UPLO
  26 *       INTEGER            N, INFO, LDA, LDAF
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IPIV( * )
  30 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *>
  40 *> DLA_SYRPVGRW computes the reciprocal pivot growth factor
  41 *> norm(A)/norm(U). The "max absolute element" norm is used. If this is
  42 *> much less than 1, the stability of the LU factorization of the
  43 *> (equilibrated) matrix A could be poor. This also means that the
  44 *> solution X, estimated condition numbers, and error bounds could be
  45 *> unreliable.
  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 number of linear equations, i.e., the order of the
  62 *>     matrix A.  N >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] INFO
  66 *> \verbatim
  67 *>          INFO is INTEGER
  68 *>     The value of INFO returned from DSYTRF, .i.e., the pivot in
  69 *>     column INFO is exactly 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] A
  73 *> \verbatim
  74 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  75 *>     On entry, the N-by-N matrix A.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LDA
  79 *> \verbatim
  80 *>          LDA is INTEGER
  81 *>     The leading dimension of the array A.  LDA >= max(1,N).
  82 *> \endverbatim
  83 *>
  84 *> \param[in] AF
  85 *> \verbatim
  86 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
  87 *>     The block diagonal matrix D and the multipliers used to
  88 *>     obtain the factor U or L as computed by DSYTRF.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] LDAF
  92 *> \verbatim
  93 *>          LDAF is INTEGER
  94 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
  95 *> \endverbatim
  96 *>
  97 *> \param[in] IPIV
  98 *> \verbatim
  99 *>          IPIV is INTEGER array, dimension (N)
 100 *>     Details of the interchanges and the block structure of D
 101 *>     as determined by DSYTRF.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] WORK
 105 *> \verbatim
 106 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
 107 *> \endverbatim
 108 *
 109 *  Authors:
 110 *  ========
 111 *
 112 *> \author Univ. of Tennessee
 113 *> \author Univ. of California Berkeley
 114 *> \author Univ. of Colorado Denver
 115 *> \author NAG Ltd.
 116 *
 117 *> \date September 2012
 118 *
 119 *> \ingroup doubleSYcomputational
 120 *
 121 *  =====================================================================
 122       DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
 123      $                                        LDAF, IPIV, WORK )
 124 *
 125 *  -- LAPACK computational routine (version 3.4.2) --
 126 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 127 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 128 *     September 2012
 129 *
 130 *     .. Scalar Arguments ..
 131       CHARACTER*1        UPLO
 132       INTEGER            N, INFO, LDA, LDAF
 133 *     ..
 134 *     .. Array Arguments ..
 135       INTEGER            IPIV( * )
 136       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
 137 *     ..
 138 *
 139 *  =====================================================================
 140 *
 141 *     .. Local Scalars ..
 142       INTEGER            NCOLS, I, J, K, KP
 143       DOUBLE PRECISION   AMAX, UMAX, RPVGRW, TMP
 144       LOGICAL            UPPER
 145 *     ..
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          ABS, MAX, MIN
 148 *     ..
 149 *     .. External Functions ..
 150       EXTERNAL           LSAME
 151       LOGICAL            LSAME
 152 *     ..
 153 *     .. Executable Statements ..
 154 *
 155       UPPER = LSAME( 'Upper', UPLO )
 156       IF ( INFO.EQ.0 ) THEN
 157          IF ( UPPER ) THEN
 158             NCOLS = 1
 159          ELSE
 160             NCOLS = N
 161          END IF
 162       ELSE
 163          NCOLS = INFO
 164       END IF
 165
 166       RPVGRW = 1.0D+0
 167       DO I = 1, 2*N
 168          WORK( I ) = 0.0D+0
 169       END DO
 170 *
 171 *     Find the max magnitude entry of each column of A.  Compute the max
 172 *     for all N columns so we can apply the pivot permutation while
 173 *     looping below.  Assume a full factorization is the common case.
 174 *
 175       IF ( UPPER ) THEN
 176          DO J = 1, N
 177             DO I = 1, J
 178                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
 179                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
 180             END DO
 181          END DO
 182       ELSE
 183          DO J = 1, N
 184             DO I = J, N
 185                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
 186                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
 187             END DO
 188          END DO
 189       END IF
 190 *
 191 *     Now find the max magnitude entry of each column of U or L.  Also
 192 *     permute the magnitudes of A above so they're in the same order as
 193 *     the factor.
 194 *
 195 *     The iteration orders and permutations were copied from dsytrs.
 196 *     Calls to SSWAP would be severe overkill.
 197 *
 198       IF ( UPPER ) THEN
 199          K = N
 200          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
 201             IF ( IPIV( K ).GT.0 ) THEN
 202 !              1x1 pivot
 203                KP = IPIV( K )
 204                IF ( KP .NE. K ) THEN
 205                   TMP = WORK( N+K )
 206                   WORK( N+K ) = WORK( N+KP )
 207                   WORK( N+KP ) = TMP
 208                END IF
 209                DO I = 1, K
 210                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
 211                END DO
 212                K = K - 1
 213             ELSE
 214 !              2x2 pivot
 215                KP = -IPIV( K )
 216                TMP = WORK( N+K-1 )
 217                WORK( N+K-1 ) = WORK( N+KP )
 218                WORK( N+KP ) = TMP
 219                DO I = 1, K-1
 220                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
 221                   WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
 222                END DO
 223                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
 224                K = K - 2
 225             END IF
 226          END DO
 227          K = NCOLS
 228          DO WHILE ( K .LE. N )
 229             IF ( IPIV( K ).GT.0 ) THEN
 230                KP = IPIV( K )
 231                IF ( KP .NE. K ) THEN
 232                   TMP = WORK( N+K )
 233                   WORK( N+K ) = WORK( N+KP )
 234                   WORK( N+KP ) = TMP
 235                END IF
 236                K = K + 1
 237             ELSE
 238                KP = -IPIV( K )
 239                TMP = WORK( N+K )
 240                WORK( N+K ) = WORK( N+KP )
 241                WORK( N+KP ) = TMP
 242                K = K + 2
 243             END IF
 244          END DO
 245       ELSE
 246          K = 1
 247          DO WHILE ( K .LE. NCOLS )
 248             IF ( IPIV( K ).GT.0 ) THEN
 249 !              1x1 pivot
 250                KP = IPIV( K )
 251                IF ( KP .NE. K ) THEN
 252                   TMP = WORK( N+K )
 253                   WORK( N+K ) = WORK( N+KP )
 254                   WORK( N+KP ) = TMP
 255                END IF
 256                DO I = K, N
 257                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
 258                END DO
 259                K = K + 1
 260             ELSE
 261 !              2x2 pivot
 262                KP = -IPIV( K )
 263                TMP = WORK( N+K+1 )
 264                WORK( N+K+1 ) = WORK( N+KP )
 265                WORK( N+KP ) = TMP
 266                DO I = K+1, N
 267                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
 268                   WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
 269                END DO
 270                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
 271                K = K + 2
 272             END IF
 273          END DO
 274          K = NCOLS
 275          DO WHILE ( K .GE. 1 )
 276             IF ( IPIV( K ).GT.0 ) THEN
 277                KP = IPIV( K )
 278                IF ( KP .NE. K ) THEN
 279                   TMP = WORK( N+K )
 280                   WORK( N+K ) = WORK( N+KP )
 281                   WORK( N+KP ) = TMP
 282                END IF
 283                K = K - 1
 284             ELSE
 285                KP = -IPIV( K )
 286                TMP = WORK( N+K )
 287                WORK( N+K ) = WORK( N+KP )
 288                WORK( N+KP ) = TMP
 289                K = K - 2
 290             ENDIF
 291          END DO
 292       END IF
 293 *
 294 *     Compute the *inverse* of the max element growth factor.  Dividing
 295 *     by zero would imply the largest entry of the factor's column is
 296 *     zero.  Than can happen when either the column of A is zero or
 297 *     massive pivots made the factor underflow to zero.  Neither counts
 298 *     as growth in itself, so simply ignore terms with zero
 299 *     denominators.
 300 *
 301       IF ( UPPER ) THEN
 302          DO I = NCOLS, N
 303             UMAX = WORK( I )
 304             AMAX = WORK( N+I )
 305             IF ( UMAX /= 0.0D+0 ) THEN
 306                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 307             END IF
 308          END DO
 309       ELSE
 310          DO I = 1, NCOLS
 311             UMAX = WORK( I )
 312             AMAX = WORK( N+I )
 313             IF ( UMAX /= 0.0D+0 ) THEN
 314                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 315             END IF
 316          END DO
 317       END IF
 318
 319       DLA_SYRPVGRW = RPVGRW
 320       END