00001 DOUBLE PRECISION FUNCTION ZLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF, 00002 $ LDAF, WORK ) 00003 * 00004 * -- LAPACK routine (version 3.2.2) -- 00005 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00006 * -- Jason Riedy of Univ. of California Berkeley. -- 00007 * -- June 2010 -- 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley and NAG Ltd. -- 00011 * 00012 IMPLICIT NONE 00013 * .. 00014 * .. Scalar Arguments .. 00015 CHARACTER*1 UPLO 00016 INTEGER NCOLS, LDA, LDAF 00017 * .. 00018 * .. Array Arguments .. 00019 COMPLEX*16 A( LDA, * ), AF( LDAF, * ) 00020 DOUBLE PRECISION WORK( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZLA_PORPVGRW computes the reciprocal pivot growth factor 00027 * norm(A)/norm(U). The "max absolute element" norm is used. If this is 00028 * much less than 1, the stability of the LU factorization of the 00029 * (equilibrated) matrix A could be poor. This also means that the 00030 * solution X, estimated condition numbers, and error bounds could be 00031 * unreliable. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * UPLO (input) CHARACTER*1 00037 * = 'U': Upper triangle of A is stored; 00038 * = 'L': Lower triangle of A is stored. 00039 * 00040 * NCOLS (input) INTEGER 00041 * The number of columns of the matrix A. NCOLS >= 0. 00042 * 00043 * A (input) COMPLEX*16 array, dimension (LDA,N) 00044 * On entry, the N-by-N matrix A. 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the array A. LDA >= max(1,N). 00048 * 00049 * AF (input) COMPLEX*16 array, dimension (LDAF,N) 00050 * The triangular factor U or L from the Cholesky factorization 00051 * A = U**T*U or A = L*L**T, as computed by ZPOTRF. 00052 * 00053 * LDAF (input) INTEGER 00054 * The leading dimension of the array AF. LDAF >= max(1,N). 00055 * 00056 * WORK (input) COMPLEX*16 array, dimension (2*N) 00057 * 00058 * ===================================================================== 00059 * 00060 * .. Local Scalars .. 00061 INTEGER I, J 00062 DOUBLE PRECISION AMAX, UMAX, RPVGRW 00063 LOGICAL UPPER 00064 COMPLEX*16 ZDUM 00065 * .. 00066 * .. External Functions .. 00067 EXTERNAL LSAME, ZLASET 00068 LOGICAL LSAME 00069 * .. 00070 * .. Intrinsic Functions .. 00071 INTRINSIC ABS, MAX, MIN, REAL, DIMAG 00072 * .. 00073 * .. Statement Functions .. 00074 DOUBLE PRECISION CABS1 00075 * .. 00076 * .. Statement Function Definitions .. 00077 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00078 * .. 00079 * .. Executable Statements .. 00080 UPPER = LSAME( 'Upper', UPLO ) 00081 * 00082 * DPOTRF will have factored only the NCOLSxNCOLS leading minor, so 00083 * we restrict the growth search to that minor and use only the first 00084 * 2*NCOLS workspace entries. 00085 * 00086 RPVGRW = 1.0D+0 00087 DO I = 1, 2*NCOLS 00088 WORK( I ) = 0.0D+0 00089 END DO 00090 * 00091 * Find the max magnitude entry of each column. 00092 * 00093 IF ( UPPER ) THEN 00094 DO J = 1, NCOLS 00095 DO I = 1, J 00096 WORK( NCOLS+J ) = 00097 $ MAX( CABS1( A( I, J ) ), WORK( NCOLS+J ) ) 00098 END DO 00099 END DO 00100 ELSE 00101 DO J = 1, NCOLS 00102 DO I = J, NCOLS 00103 WORK( NCOLS+J ) = 00104 $ MAX( CABS1( A( I, J ) ), WORK( NCOLS+J ) ) 00105 END DO 00106 END DO 00107 END IF 00108 * 00109 * Now find the max magnitude entry of each column of the factor in 00110 * AF. No pivoting, so no permutations. 00111 * 00112 IF ( LSAME( 'Upper', UPLO ) ) THEN 00113 DO J = 1, NCOLS 00114 DO I = 1, J 00115 WORK( J ) = MAX( CABS1( AF( I, J ) ), WORK( J ) ) 00116 END DO 00117 END DO 00118 ELSE 00119 DO J = 1, NCOLS 00120 DO I = J, NCOLS 00121 WORK( J ) = MAX( CABS1( AF( I, J ) ), WORK( J ) ) 00122 END DO 00123 END DO 00124 END IF 00125 * 00126 * Compute the *inverse* of the max element growth factor. Dividing 00127 * by zero would imply the largest entry of the factor's column is 00128 * zero. Than can happen when either the column of A is zero or 00129 * massive pivots made the factor underflow to zero. Neither counts 00130 * as growth in itself, so simply ignore terms with zero 00131 * denominators. 00132 * 00133 IF ( LSAME( 'Upper', UPLO ) ) THEN 00134 DO I = 1, NCOLS 00135 UMAX = WORK( I ) 00136 AMAX = WORK( NCOLS+I ) 00137 IF ( UMAX /= 0.0D+0 ) THEN 00138 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00139 END IF 00140 END DO 00141 ELSE 00142 DO I = 1, NCOLS 00143 UMAX = WORK( I ) 00144 AMAX = WORK( NCOLS+I ) 00145 IF ( UMAX /= 0.0D+0 ) THEN 00146 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00147 END IF 00148 END DO 00149 END IF 00150 00151 ZLA_PORPVGRW = RPVGRW 00152 END
1.7.2