Return to dla_porpvgrw.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dla_porpvgrw.f, revision 1.4

1.1       bertrand    1:       DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF, 
                      2:      $                                        LDAF, WORK )
                      3: *
                      4: *     -- LAPACK routine (version 3.2.2)                                 --
                      5: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
                      6: *     -- Jason Riedy of Univ. of California Berkeley.                 --
                      7: *     -- June 2010                                                    --
                      8: *
                      9: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
                     10: *     -- Univ. of California Berkeley and NAG Ltd.                    --
                     11: *
                     12:       IMPLICIT NONE
                     13: *     ..
                     14: *     .. Scalar Arguments ..
                     15:       CHARACTER*1        UPLO
                     16:       INTEGER            NCOLS, LDA, LDAF
                     17: *     ..
                     18: *     .. Array Arguments ..
                     19:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
                     20: *     ..
                     21: *
                     22: *  Purpose
                     23: *  =======
                     24: * 
                     25: *  DLA_PORPVGRW computes the reciprocal pivot growth factor
                     26: *  norm(A)/norm(U). The "max absolute element" norm is used. If this is
                     27: *  much less than 1, the stability of the LU factorization of the
                     28: *  (equilibrated) matrix A could be poor. This also means that the
                     29: *  solution X, estimated condition numbers, and error bounds could be
                     30: *  unreliable.
                     31: *
                     32: *  Arguments
                     33: *  =========
                     34: *
                     35: *     UPLO    (input) CHARACTER*1
                     36: *       = 'U':  Upper triangle of A is stored;
                     37: *       = 'L':  Lower triangle of A is stored.
                     38: *
                     39: *     NCOLS   (input) INTEGER
                     40: *     The number of columns of the matrix A. NCOLS >= 0.
                     41: *
                     42: *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
                     43: *     On entry, the N-by-N matrix A.
                     44: *
                     45: *     LDA     (input) INTEGER
                     46: *     The leading dimension of the array A.  LDA >= max(1,N).
                     47: *
                     48: *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
                     49: *     The triangular factor U or L from the Cholesky factorization
                     50: *     A = U**T*U or A = L*L**T, as computed by DPOTRF.
                     51: *
                     52: *     LDAF    (input) INTEGER
                     53: *     The leading dimension of the array AF.  LDAF >= max(1,N).
                     54: *
                     55: *     WORK    (input) DOUBLE PRECISION array, dimension (2*N)
                     56: *
                     57: *  =====================================================================
                     58: *
                     59: *     .. Local Scalars ..
                     60:       INTEGER            I, J
                     61:       DOUBLE PRECISION   AMAX, UMAX, RPVGRW
                     62:       LOGICAL            UPPER
                     63: *     ..
                     64: *     .. Intrinsic Functions ..
                     65:       INTRINSIC          ABS, MAX, MIN
                     66: *     ..
                     67: *     .. External Functions ..
                     68:       EXTERNAL           LSAME, DLASET
                     69:       LOGICAL            LSAME
                     70: *     ..
                     71: *     .. Executable Statements ..
                     72: *
                     73:       UPPER = LSAME( 'Upper', UPLO )
                     74: *
                     75: *     DPOTRF will have factored only the NCOLSxNCOLS leading minor, so
                     76: *     we restrict the growth search to that minor and use only the first
                     77: *     2*NCOLS workspace entries.
                     78: *
                     79:       RPVGRW = 1.0D+0
                     80:       DO I = 1, 2*NCOLS
                     81:          WORK( I ) = 0.0D+0
                     82:       END DO
                     83: *
                     84: *     Find the max magnitude entry of each column.
                     85: *
                     86:       IF ( UPPER ) THEN
                     87:          DO J = 1, NCOLS
                     88:             DO I = 1, J
                     89:                WORK( NCOLS+J ) =
                     90:      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
                     91:             END DO
                     92:          END DO
                     93:       ELSE
                     94:          DO J = 1, NCOLS
                     95:             DO I = J, NCOLS
                     96:                WORK( NCOLS+J ) =
                     97:      $              MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
                     98:             END DO
                     99:          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