LAPACK: SRC/cla_herpvgrw.f Source File

Go to the documentation of this file.
00001       REAL FUNCTION CLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, IPIV,
00002      $                            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            N, INFO, LDA, LDAF
00017 *     ..
00018 *     .. Array Arguments ..
00019       INTEGER            IPIV( * )
00020       COMPLEX            A( LDA, * ), AF( LDAF, * )
00021       REAL               WORK( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 * 
00027 *  CLA_HERPVGRW computes the reciprocal pivot growth factor
00028 *  norm(A)/norm(U). The "max absolute element" norm is used. If this is
00029 *  much less than 1, the stability of the LU factorization of the
00030 *  (equilibrated) matrix A could be poor. This also means that the
00031 *  solution X, estimated condition numbers, and error bounds could be
00032 *  unreliable.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *     UPLO    (input) CHARACTER*1
00038 *       = 'U':  Upper triangle of A is stored;
00039 *       = 'L':  Lower triangle of A is stored.
00040 *
00041 *     N       (input) INTEGER
00042 *     The number of linear equations, i.e., the order of the
00043 *     matrix A.  N >= 0.
00044 *
00045 *     INFO    (input) INTEGER
00046 *     The value of INFO returned from SSYTRF, .i.e., the pivot in
00047 *     column INFO is exactly 0.
00048 *
00049 *     NCOLS   (input) INTEGER
00050 *     The number of columns of the matrix A. NCOLS >= 0.
00051 *
00052 *     A       (input) COMPLEX array, dimension (LDA,N)
00053 *     On entry, the N-by-N matrix A.
00054 *
00055 *     LDA     (input) INTEGER
00056 *     The leading dimension of the array A.  LDA >= max(1,N).
00057 *
00058 *     AF      (input) COMPLEX array, dimension (LDAF,N)
00059 *     The block diagonal matrix D and the multipliers used to
00060 *     obtain the factor U or L as computed by CHETRF.
00061 *
00062 *     LDAF    (input) INTEGER
00063 *     The leading dimension of the array AF.  LDAF >= max(1,N).
00064 *
00065 *     IPIV    (input) INTEGER array, dimension (N)
00066 *     Details of the interchanges and the block structure of D
00067 *     as determined by CHETRF.
00068 *
00069 *     WORK    (input) COMPLEX array, dimension (2*N)
00070 *
00071 *  =====================================================================
00072 *
00073 *     .. Local Scalars ..
00074       INTEGER            NCOLS, I, J, K, KP
00075       REAL               AMAX, UMAX, RPVGRW, TMP
00076       LOGICAL            UPPER, LSAME
00077       COMPLEX            ZDUM
00078 *     ..
00079 *     .. External Functions ..
00080       EXTERNAL           LSAME, CLASET
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          ABS, REAL, AIMAG, MAX, MIN
00084 *     ..
00085 *     .. Statement Functions ..
00086       REAL               CABS1
00087 *     ..
00088 *     .. Statement Function Definitions ..
00089       CABS1( ZDUM ) = ABS( REAL ( ZDUM ) ) + ABS( AIMAG ( ZDUM ) )
00090 *     ..
00091 *     .. Executable Statements ..
00092 *
00093       UPPER = LSAME( 'Upper', UPLO )
00094       IF ( INFO.EQ.0 ) THEN
00095          IF (UPPER) THEN
00096             NCOLS = 1
00097          ELSE
00098             NCOLS = N
00099          END IF
00100       ELSE
00101          NCOLS = INFO
00102       END IF
00103 
00104       RPVGRW = 1.0
00105       DO I = 1, 2*N
00106          WORK( I ) = 0.0
00107       END DO
00108 *
00109 *     Find the max magnitude entry of each column of A.  Compute the max
00110 *     for all N columns so we can apply the pivot permutation while
00111 *     looping below.  Assume a full factorization is the common case.
00112 *
00113       IF ( UPPER ) THEN
00114          DO J = 1, N
00115             DO I = 1, J
00116                WORK( N+I ) = MAX( CABS1( A( I,J ) ), WORK( N+I ) )
00117                WORK( N+J ) = MAX( CABS1( A( I,J ) ), WORK( N+J ) )
00118             END DO
00119          END DO
00120       ELSE
00121          DO J = 1, N
00122             DO I = J, N
00123                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
00124                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
00125             END DO
00126          END DO
00127       END IF
00128 *
00129 *     Now find the max magnitude entry of each column of U or L.  Also
00130 *     permute the magnitudes of A above so they're in the same order as
00131 *     the factor.
00132 *
00133 *     The iteration orders and permutations were copied from csytrs.
00134 *     Calls to SSWAP would be severe overkill.
00135 *
00136       IF ( UPPER ) THEN
00137          K = N
00138          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
00139             IF ( IPIV( K ).GT.0 ) THEN
00140 !              1x1 pivot
00141                KP = IPIV( K )
00142                IF ( KP .NE. K ) THEN
00143                   TMP = WORK( N+K )
00144                   WORK( N+K ) = WORK( N+KP )
00145                   WORK( N+KP ) = TMP
00146                END IF
00147                DO I = 1, K
00148                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00149                END DO
00150                K = K - 1
00151             ELSE
00152 !              2x2 pivot
00153                KP = -IPIV( K )
00154                TMP = WORK( N+K-1 )
00155                WORK( N+K-1 ) = WORK( N+KP )
00156                WORK( N+KP ) = TMP
00157                DO I = 1, K-1
00158                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00159                   WORK( K-1 ) =
00160      $                 MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
00161                END DO
00162                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00163                K = K - 2
00164             END IF
00165          END DO
00166          K = NCOLS
00167          DO WHILE ( K .LE. N )
00168             IF ( IPIV( K ).GT.0 ) THEN
00169                KP = IPIV( K )
00170                IF ( KP .NE. K ) THEN
00171                   TMP = WORK( N+K )
00172                   WORK( N+K ) = WORK( N+KP )
00173                   WORK( N+KP ) = TMP
00174                END IF
00175                K = K + 1
00176             ELSE
00177                KP = -IPIV( K )
00178                TMP = WORK( N+K )
00179                WORK( N+K ) = WORK( N+KP )
00180                WORK( N+KP ) = TMP
00181                K = K + 2
00182             END IF
00183          END DO
00184       ELSE
00185          K = 1
00186          DO WHILE ( K .LE. NCOLS )
00187             IF ( IPIV( K ).GT.0 ) THEN
00188 !              1x1 pivot
00189                KP = IPIV( K )
00190                IF ( KP .NE. K ) THEN
00191                   TMP = WORK( N+K )
00192                   WORK( N+K ) = WORK( N+KP )
00193                   WORK( N+KP ) = TMP
00194                END IF
00195                DO I = K, N
00196                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00197                END DO
00198                K = K + 1
00199             ELSE
00200 !              2x2 pivot
00201                KP = -IPIV( K )
00202                TMP = WORK( N+K+1 )
00203                WORK( N+K+1 ) = WORK( N+KP )
00204                WORK( N+KP ) = TMP
00205                DO I = K+1, N
00206                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00207                   WORK( K+1 ) =
00208      $                 MAX( CABS1( AF( I, K+1 ) ) , WORK( K+1 ) )
00209                END DO
00210                WORK(K) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00211                K = K + 2
00212             END IF
00213          END DO
00214          K = NCOLS
00215          DO WHILE ( K .GE. 1 )
00216             IF ( IPIV( K ).GT.0 ) THEN
00217                KP = IPIV( K )
00218                IF ( KP .NE. K ) THEN
00219                   TMP = WORK( N+K )
00220                   WORK( N+K ) = WORK( N+KP )
00221                   WORK( N+KP ) = TMP
00222                END IF
00223                K = K - 1
00224             ELSE
00225                KP = -IPIV( K )
00226                TMP = WORK( N+K )
00227                WORK( N+K ) = WORK( N+KP )
00228                WORK( N+KP ) = TMP
00229                K = K - 2
00230             ENDIF
00231          END DO
00232       END IF
00233 *
00234 *     Compute the *inverse* of the max element growth factor.  Dividing
00235 *     by zero would imply the largest entry of the factor's column is
00236 *     zero.  Than can happen when either the column of A is zero or
00237 *     massive pivots made the factor underflow to zero.  Neither counts
00238 *     as growth in itself, so simply ignore terms with zero
00239 *     denominators.
00240 *
00241       IF ( UPPER ) THEN
00242          DO I = NCOLS, N
00243             UMAX = WORK( I )
00244             AMAX = WORK( N+I )
00245             IF ( UMAX /= 0.0 ) THEN
00246                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00247             END IF
00248          END DO
00249       ELSE
00250          DO I = 1, NCOLS
00251             UMAX = WORK( I )
00252             AMAX = WORK( N+I )
00253             IF ( UMAX /= 0.0 ) THEN
00254                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00255             END IF
00256          END DO
00257       END IF
00258 
00259       CLA_HERPVGRW = RPVGRW
00260       END