SRC/dla_porpvgrw.f

   1 *> \brief \b DLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite 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_PORPVGRW + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_porpvgrw.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_porpvgrw.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_porpvgrw.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF,
  22 *                                               LDAF, WORK )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER*1        UPLO
  26 *       INTEGER            NCOLS, LDA, LDAF
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *>
  39 *> DLA_PORPVGRW computes the reciprocal pivot growth factor
  40 *> norm(A)/norm(U). The "max absolute element" norm is used. If this is
  41 *> much less than 1, the stability of the LU factorization of the
  42 *> (equilibrated) matrix A could be poor. This also means that the
  43 *> solution X, estimated condition numbers, and error bounds could be
  44 *> unreliable.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] UPLO
  51 *> \verbatim
  52 *>          UPLO is CHARACTER*1
  53 *>       = 'U':  Upper triangle of A is stored;
  54 *>       = 'L':  Lower triangle of A is stored.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] NCOLS
  58 *> \verbatim
  59 *>          NCOLS is INTEGER
  60 *>     The number of columns of the matrix A. NCOLS >= 0.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] A
  64 *> \verbatim
  65 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  66 *>     On entry, the N-by-N matrix A.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] LDA
  70 *> \verbatim
  71 *>          LDA is INTEGER
  72 *>     The leading dimension of the array A.  LDA >= max(1,N).
  73 *> \endverbatim
  74 *>
  75 *> \param[in] AF
  76 *> \verbatim
  77 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
  78 *>     The triangular factor U or L from the Cholesky factorization
  79 *>     A = U**T*U or A = L*L**T, as computed by DPOTRF.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LDAF
  83 *> \verbatim
  84 *>          LDAF is INTEGER
  85 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
  86 *> \endverbatim
  87 *>
  88 *> \param[in] WORK
  89 *> \verbatim
  90 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
  91 *> \endverbatim
  92 *
  93 *  Authors:
  94 *  ========
  95 *
  96 *> \author Univ. of Tennessee
  97 *> \author Univ. of California Berkeley
  98 *> \author Univ. of Colorado Denver
  99 *> \author NAG Ltd.
 100 *
 101 *> \date September 2012
 102 *
 103 *> \ingroup doublePOcomputational
 104 *
 105 *  =====================================================================
 106       DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF,
 107      $                                        LDAF, WORK )
 108 *
 109 *  -- LAPACK computational routine (version 3.4.2) --
 110 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 111 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 112 *     September 2012
 113 *
 114 *     .. Scalar Arguments ..
 115       CHARACTER*1        UPLO
 116       INTEGER            NCOLS, LDA, LDAF
 117 *     ..
 118 *     .. Array Arguments ..
 119       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
 120 *     ..
 121 *
 122 *  =====================================================================
 123 *
 124 *     .. Local Scalars ..
 125       INTEGER            I, J
 126       DOUBLE PRECISION   AMAX, UMAX, RPVGRW
 127       LOGICAL            UPPER
 128 *     ..
 129 *     .. Intrinsic Functions ..
 130       INTRINSIC          ABS, MAX, MIN
 131 *     ..
 132 *     .. External Functions ..
 133       EXTERNAL           LSAME
 134       LOGICAL            LSAME
 135 *     ..
 136 *     .. Executable Statements ..
 137 *
 138       UPPER = LSAME( 'Upper', UPLO )
 139 *
 140 *     DPOTRF will have factored only the NCOLSxNCOLS leading minor, so
 141 *     we restrict the growth search to that minor and use only the first
 142 *     2*NCOLS workspace entries.
 143 *
 144       RPVGRW = 1.0D+0
 145       DO I = 1, 2*NCOLS
 146          WORK( I ) = 0.0D+0
 147       END DO
 148 *
 149 *     Find the max magnitude entry of each column.
 150 *
 151       IF ( UPPER ) THEN
 152          DO J = 1, NCOLS
 153             DO I = 1, J
 154                WORK( NCOLS+J ) =
 155      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
 156             END DO
 157          END DO
 158       ELSE
 159          DO J = 1, NCOLS
 160             DO I = J, NCOLS
 161                WORK( NCOLS+J ) =
 162      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
 163             END DO
 164          END DO
 165       END IF
 166 *
 167 *     Now find the max magnitude entry of each column of the factor in
 168 *     AF.  No pivoting, so no permutations.
 169 *
 170       IF ( LSAME( 'Upper', UPLO ) ) THEN
 171          DO J = 1, NCOLS
 172             DO I = 1, J
 173                WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
 174             END DO
 175          END DO
 176       ELSE
 177          DO J = 1, NCOLS
 178             DO I = J, NCOLS
 179                WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
 180             END DO
 181          END DO
 182       END IF
 183 *
 184 *     Compute the *inverse* of the max element growth factor.  Dividing
 185 *     by zero would imply the largest entry of the factor's column is
 186 *     zero.  Than can happen when either the column of A is zero or
 187 *     massive pivots made the factor underflow to zero.  Neither counts
 188 *     as growth in itself, so simply ignore terms with zero
 189 *     denominators.
 190 *
 191       IF ( LSAME( 'Upper', UPLO ) ) THEN
 192          DO I = 1, NCOLS
 193             UMAX = WORK( I )
 194             AMAX = WORK( NCOLS+I )
 195             IF ( UMAX /= 0.0D+0 ) THEN
 196                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 197             END IF
 198          END DO
 199       ELSE
 200          DO I = 1, NCOLS
 201             UMAX = WORK( I )
 202             AMAX = WORK( NCOLS+I )
 203             IF ( UMAX /= 0.0D+0 ) THEN
 204                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
 205             END IF
 206          END DO
 207       END IF
 208
 209       DLA_PORPVGRW = RPVGRW
 210       END