SRC/zla_syrpvgrw.f

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