001:       DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF, 
002:      $                                        LDAF, WORK )
003: *
004: *     -- LAPACK routine (version 3.2.1)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- April 2009                                                   --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       CHARACTER*1        UPLO
016:       INTEGER            NCOLS, LDA, LDAF
017: *     ..
018: *     .. Array Arguments ..
019:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: * 
025: *  DLA_PORPVGRW computes the reciprocal pivot growth factor
026: *  norm(A)/norm(U). The "max absolute element" norm is used. If this is
027: *  much less than 1, the stability of the LU factorization of the
028: *  (equilibrated) matrix A could be poor. This also means that the
029: *  solution X, estimated condition numbers, and error bounds could be
030: *  unreliable.
031: *
032: *  Arguments
033: *  =========
034: *
035: *     UPLO    (input) CHARACTER*1
036: *       = 'U':  Upper triangle of A is stored;
037: *       = 'L':  Lower triangle of A is stored.
038: *
039: *     NCOLS   (input) INTEGER
040: *     The number of columns of the matrix A. NCOLS >= 0.
041: *
042: *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
043: *     On entry, the N-by-N matrix A.
044: *
045: *     LDA     (input) INTEGER
046: *     The leading dimension of the array A.  LDA >= max(1,N).
047: *
048: *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
049: *     The triangular factor U or L from the Cholesky factorization
050: *     A = U**T*U or A = L*L**T, as computed by DPOTRF.
051: *
052: *     LDAF    (input) INTEGER
053: *     The leading dimension of the array AF.  LDAF >= max(1,N).
054: *
055: *     WORK    (input) DOUBLE PRECISION array, dimension (2*N)
056: *
057: *  =====================================================================
058: *
059: *     .. Local Scalars ..
060:       INTEGER            I, J
061:       DOUBLE PRECISION   AMAX, UMAX, RPVGRW
062:       LOGICAL            UPPER
063: *     ..
064: *     .. Intrinsic Functions ..
065:       INTRINSIC          ABS, MAX, MIN
066: *     ..
067: *     .. External Functions ..
068:       EXTERNAL           LSAME, DLASET
069:       LOGICAL            LSAME
070: *     ..
071: *     .. Executable Statements ..
072: *
073:       UPPER = LSAME( 'Upper', UPLO )
074: *
075: *     DPOTRF will have factored only the NCOLSxNCOLS leading minor, so
076: *     we restrict the growth search to that minor and use only the first
077: *     2*NCOLS workspace entries.
078: *
079:       RPVGRW = 1.0D+0
080:       DO I = 1, 2*NCOLS
081:          WORK( I ) = 0.0D+0
082:       END DO
083: *
084: *     Find the max magnitude entry of each column.
085: *
086:       IF ( UPPER ) THEN
087:          DO J = 1, NCOLS
088:             DO I = 1, J
089:                WORK( NCOLS+J ) =
090:      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
091:             END DO
092:          END DO
093:       ELSE
094:          DO J = 1, NCOLS
095:             DO I = J, NCOLS
096:                WORK( NCOLS+J ) =
097:      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
098:             END DO
099:          END DO
100:       END IF
101: *
102: *     Now find the max magnitude entry of each column of the factor in
103: *     AF.  No pivoting, so no permutations.
104: *
105:       IF ( LSAME( 'Upper', UPLO ) ) THEN
106:          DO J = 1, NCOLS
107:             DO I = 1, J
108:                WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
109:             END DO
110:          END DO
111:       ELSE
112:          DO J = 1, NCOLS
113:             DO I = J, NCOLS
114:                WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
115:             END DO
116:          END DO
117:       END IF
118: *
119: *     Compute the *inverse* of the max element growth factor.  Dividing
120: *     by zero would imply the largest entry of the factor's column is
121: *     zero.  Than can happen when either the column of A is zero or
122: *     massive pivots made the factor underflow to zero.  Neither counts
123: *     as growth in itself, so simply ignore terms with zero
124: *     denominators.
125: *
126:       IF ( LSAME( 'Upper', UPLO ) ) THEN
127:          DO I = 1, NCOLS
128:             UMAX = WORK( I )
129:             AMAX = WORK( NCOLS+I )
130:             IF ( UMAX /= 0.0D+0 ) THEN
131:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
132:             END IF
133:          END DO
134:       ELSE
135:          DO I = 1, NCOLS
136:             UMAX = WORK( I )
137:             AMAX = WORK( NCOLS+I )
138:             IF ( UMAX /= 0.0D+0 ) THEN
139:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
140:             END IF
141:          END DO
142:       END IF
143: 
144:       DLA_PORPVGRW = RPVGRW
145:       END FUNCTION
146: