SRC/zla_herpvgrw.f

   1 *> \brief \b ZLA_HERPVGRW
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLA_HERPVGRW + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_herpvgrw.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_herpvgrw.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_herpvgrw.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       DOUBLE PRECISION FUNCTION ZLA_HERPVGRW( 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 *       COMPLEX*16         A( LDA, * ), AF( LDAF, * )
  31 *       DOUBLE PRECISION   WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *>
  41 *> ZLA_HERPVGRW 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 ZHETRF, .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 ZHETRF.
  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 ZHETRF.
 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 June 2016
 119 *
 120 *> \ingroup complex16HEcomputational
 121 *
 122 *  =====================================================================
 123       DOUBLE PRECISION FUNCTION ZLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF,
 124      $                                        LDAF, IPIV, WORK )
 125 *
 126 *  -- LAPACK computational routine (version 3.6.1) --
 127 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 128 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 129 *     June 2016
 130 *
 131 *     .. Scalar Arguments ..
 132       CHARACTER*1        UPLO
 133       INTEGER            N, INFO, LDA, LDAF
 134 *     ..
 135 *     .. Array Arguments ..
 136       INTEGER            IPIV( * )
 137       COMPLEX*16         A( LDA, * ), AF( LDAF, * )
 138       DOUBLE PRECISION   WORK( * )
 139 *     ..
 140 *
 141 *  =====================================================================
 142 *
 143 *     .. Local Scalars ..
 144       INTEGER            NCOLS, I, J, K, KP
 145       DOUBLE PRECISION   AMAX, UMAX, RPVGRW, TMP
 146       LOGICAL            UPPER, LSAME
 147       COMPLEX*16         ZDUM
 148 *     ..
 149 *     .. External Functions ..
 150       EXTERNAL           LSAME, ZLASET
 151 *     ..
 152 *     .. Intrinsic Functions ..
 153       INTRINSIC          ABS, REAL, DIMAG, MAX, MIN
 154 *     ..
 155 *     .. Statement Functions ..
 156       DOUBLE PRECISION   CABS1
 157 *     ..
 158 *     .. Statement Function Definitions ..
 159       CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
 160 *     ..
 161 *     .. Executable Statements ..
 162 *
 163       UPPER = LSAME( 'Upper', UPLO )
 164       IF ( INFO.EQ.0 ) THEN
 165          IF (UPPER) THEN
 166             NCOLS = 1
 167          ELSE
 168             NCOLS = N
 169          END IF
 170       ELSE
 171          NCOLS = INFO
 172       END IF
 173
 174       RPVGRW = 1.0D+0
 175       DO I = 1, 2*N
 176          WORK( I ) = 0.0D+0
 177       END DO
 178 *
 179 *     Find the max magnitude entry of each column of A.  Compute the max
 180 *     for all N columns so we can apply the pivot permutation while
 181 *     looping below.  Assume a full factorization is the common case.
 182 *
 183       IF ( UPPER ) THEN
 184          DO J = 1, N
 185             DO I = 1, J
 186                WORK( N+I ) = MAX( CABS1( A( I,J ) ), WORK( N+I ) )
 187                WORK( N+J ) = MAX( CABS1( A( I,J ) ), WORK( N+J ) )
 188             END DO
 189          END DO
 190       ELSE
 191          DO J = 1, N
 192             DO I = J, N
 193                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
 194                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
 195             END DO
 196          END DO
 197       END IF
 198 *
 199 *     Now find the max magnitude entry of each column of U or L.  Also
 200 *     permute the magnitudes of A above so they're in the same order as
 201 *     the factor.
 202 *
 203 *     The iteration orders and permutations were copied from zsytrs.
 204 *     Calls to SSWAP would be severe overkill.
 205 *
 206       IF ( UPPER ) THEN
 207          K = N
 208          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
 209             IF ( IPIV( K ).GT.0 ) THEN
 210 !              1x1 pivot
 211                KP = IPIV( K )
 212                IF ( KP .NE. K ) THEN
 213                   TMP = WORK( N+K )
 214                   WORK( N+K ) = WORK( N+KP )
 215                   WORK( N+KP ) = TMP
 216                END IF
 217                DO I = 1, K
 218                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
 219                END DO
 220                K = K - 1
 221             ELSE
 222 !              2x2 pivot
 223                KP = -IPIV( K )
 224                TMP = WORK( N+K-1 )
 225                WORK( N+K-1 ) = WORK( N+KP )
 226                WORK( N+KP ) = TMP
 227                DO I = 1, K-1
 228                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
 229                   WORK( K-1 ) =
 230      $                 MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
 231                END DO
 232                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
 233                K = K - 2
 234             END IF
 235          END DO
 236          K = NCOLS
 237          DO WHILE ( K .LE. N )
 238             IF ( IPIV( K ).GT.0 ) THEN
 239                KP = IPIV( K )
 240                IF ( KP .NE. K ) THEN
 241                   TMP = WORK( N+K )
 242                   WORK( N+K ) = WORK( N+KP )
 243                   WORK( N+KP ) = TMP
 244                END IF
 245                K = K + 1
 246             ELSE
 247                KP = -IPIV( K )
 248                TMP = WORK( N+K )
 249                WORK( N+K ) = WORK( N+KP )
 250                WORK( N+KP ) = TMP
 251                K = K + 2
 252             END IF
 253          END DO
 254       ELSE
 255          K = 1
 256          DO WHILE ( K .LE. NCOLS )
 257             IF ( IPIV( K ).GT.0 ) THEN
 258 !              1x1 pivot
 259                KP = IPIV( K )
 260                IF ( KP .NE. K ) THEN
 261                   TMP = WORK( N+K )
 262                   WORK( N+K ) = WORK( N+KP )
 263                   WORK( N+KP ) = TMP
 264                END IF
 265                DO I = K, N
 266                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
 267                END DO
 268                K = K + 1
 269             ELSE
 270 !              2x2 pivot
 271                KP = -IPIV( K )
 272                TMP = WORK( N+K+1 )
 273                WORK( N+K+1 ) = WORK( N+KP )
 274                WORK( N+KP ) = TMP
 275                DO I = K+1, N
 276                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
 277                   WORK( K+1 ) =
 278      $                 MAX( CABS1( AF( I, K+1 ) ) , WORK( K+1 ) )
 279                END DO
 280                WORK(K) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
 281                K = K + 2
 282             END IF
 283          END DO
 284          K = NCOLS
 285          DO WHILE ( K .GE. 1 )
 286             IF ( IPIV( K ).GT.0 ) THEN
 287                KP = IPIV( K )
 288                IF ( KP .NE. K ) THEN
 289                   TMP = WORK( N+K )
 290                   WORK( N+K ) = WORK( N+KP )
 291                   WORK( N+KP ) = TMP
 292                END IF
 293                K = K - 1
 294             ELSE
 295                KP = -IPIV( K )
 296                TMP = WORK( N+K )
 297                WORK( N+K ) = WORK( N+KP )
 298                WORK( N+KP ) = TMP
 299                K = K - 2
 300             ENDIF
 301          END DO
 302       END IF
 303 *
 304 *     Compute the *inverse* of the max element growth factor.  Dividing
 305 *     by zero would imply the largest entry of the factor's column is
 306 *     zero.  Than can happen when either the column of A is zero or
 307 *     massive pivots made the factor underflow to zero.  Neither counts
 308 *     as growth in itself, so simply ignore terms with zero
 309 *     denominators.
 310 *
 311       IF ( UPPER ) THEN
 312          DO I = NCOLS, N
 313             UMAX = WORK( I )
 314             AMAX = WORK( N+I )
 315             IF ( UMAX /= 0.0D+0 ) THEN
 316                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 317             END IF
 318          END DO
 319       ELSE
 320          DO I = 1, NCOLS
 321             UMAX = WORK( I )
 322             AMAX = WORK( N+I )
 323             IF ( UMAX /= 0.0D+0 ) THEN
 324                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 325             END IF
 326          END DO
 327       END IF
 328
 329       ZLA_HERPVGRW = RPVGRW
 330       END