001:       DOUBLE PRECISION FUNCTION ZLA_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:       COMPLEX*16         A( LDA, * ), AF( LDAF, * )
020:       DOUBLE PRECISION   WORK( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: * 
026: *  ZLA_PORPVGRW computes the reciprocal pivot growth factor
027: *  norm(A)/norm(U). The "max absolute element" norm is used. If this is
028: *  much less than 1, the stability of the LU factorization of the
029: *  (equilibrated) matrix A could be poor. This also means that the
030: *  solution X, estimated condition numbers, and error bounds could be
031: *  unreliable.
032: *
033: *  Arguments
034: *  =========
035: *
036: *     UPLO    (input) CHARACTER*1
037: *       = 'U':  Upper triangle of A is stored;
038: *       = 'L':  Lower triangle of A is stored.
039: *
040: *     NCOLS   (input) INTEGER
041: *     The number of columns of the matrix A. NCOLS >= 0.
042: *
043: *     A       (input) COMPLEX*16 array, dimension (LDA,N)
044: *     On entry, the N-by-N matrix A.
045: *
046: *     LDA     (input) INTEGER
047: *     The leading dimension of the array A.  LDA >= max(1,N).
048: *
049: *     AF      (input) COMPLEX*16 array, dimension (LDAF,N)
050: *     The triangular factor U or L from the Cholesky factorization
051: *     A = U**T*U or A = L*L**T, as computed by ZPOTRF.
052: *
053: *     LDAF    (input) INTEGER
054: *     The leading dimension of the array AF.  LDAF >= max(1,N).
055: *
056: *     WORK    (input) COMPLEX*16 array, dimension (2*N)
057: *
058: *  =====================================================================
059: *
060: *     .. Local Scalars ..
061:       INTEGER            I, J
062:       DOUBLE PRECISION   AMAX, UMAX, RPVGRW
063:       LOGICAL            UPPER
064:       COMPLEX*16         ZDUM
065: *     ..
066: *     .. External Functions ..
067:       EXTERNAL           LSAME, ZLASET
068:       LOGICAL            LSAME
069: *     ..
070: *     .. Intrinsic Functions ..
071:       INTRINSIC          ABS, MAX, MIN, REAL, DIMAG
072: *     ..
073: *     .. Statement Functions ..
074:       DOUBLE PRECISION   CABS1
075: *     ..
076: *     .. Statement Function Definitions ..
077:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
078: *     ..
079: *     .. Executable Statements ..
080:       UPPER = LSAME( 'Upper', UPLO )
081: *
082: *     DPOTRF will have factored only the NCOLSxNCOLS leading minor, so
083: *     we restrict the growth search to that minor and use only the first
084: *     2*NCOLS workspace entries.
085: *
086:       RPVGRW = 1.0D+0
087:       DO I = 1, 2*NCOLS
088:          WORK( I ) = 0.0D+0
089:       END DO
090: *
091: *     Find the max magnitude entry of each column.
092: *
093:       IF ( UPPER ) THEN
094:          DO J = 1, NCOLS
095:             DO I = 1, J
096:                WORK( NCOLS+J ) =
097:      $              MAX( CABS1( A( I, J ) ), WORK( NCOLS+J ) )
098:             END DO
099:          END DO
100:       ELSE
101:          DO J = 1, NCOLS
102:             DO I = J, NCOLS
103:                WORK( NCOLS+J ) =
104:      $              MAX( CABS1( A( I, J ) ), WORK( NCOLS+J ) )
105:             END DO
106:          END DO
107:       END IF
108: *
109: *     Now find the max magnitude entry of each column of the factor in
110: *     AF.  No pivoting, so no permutations.
111: *
112:       IF ( LSAME( 'Upper', UPLO ) ) THEN
113:          DO J = 1, NCOLS
114:             DO I = 1, J
115:                WORK( J ) = MAX( CABS1( AF( I, J ) ), WORK( J ) )
116:             END DO
117:          END DO
118:       ELSE
119:          DO J = 1, NCOLS
120:             DO I = J, NCOLS
121:                WORK( J ) = MAX( CABS1( AF( I, J ) ), WORK( J ) )
122:             END DO
123:          END DO
124:       END IF
125: *
126: *     Compute the *inverse* of the max element growth factor.  Dividing
127: *     by zero would imply the largest entry of the factor's column is
128: *     zero.  Than can happen when either the column of A is zero or
129: *     massive pivots made the factor underflow to zero.  Neither counts
130: *     as growth in itself, so simply ignore terms with zero
131: *     denominators.
132: *
133:       IF ( LSAME( 'Upper', UPLO ) ) THEN
134:          DO I = 1, NCOLS
135:             UMAX = WORK( I )
136:             AMAX = WORK( NCOLS+I )
137:             IF ( UMAX /= 0.0D+0 ) THEN
138:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
139:             END IF
140:          END DO
141:       ELSE
142:          DO I = 1, NCOLS
143:             UMAX = WORK( I )
144:             AMAX = WORK( NCOLS+I )
145:             IF ( UMAX /= 0.0D+0 ) THEN
146:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
147:             END IF
148:          END DO
149:       END IF
150: 
151:       ZLA_PORPVGRW = RPVGRW
152:       END FUNCTION
153: