LAPACK: SRC/zla_syrpvgrw.f Source File

Go to the documentation of this file.
00001       DOUBLE PRECISION FUNCTION ZLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
00002      $                                        LDAF, IPIV, 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       COMPLEX*16         A( LDA, * ), AF( LDAF, * )
00020       DOUBLE PRECISION   WORK( * )
00021       INTEGER            IPIV( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 * 
00027 *  ZLA_SYRPVGRW 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 ZSYTRF, .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*16 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*16 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 ZSYTRF.
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 ZSYTRF.
00068 *
00069 *     WORK    (input) COMPLEX*16 array, dimension (2*N)
00070 *
00071 *  =====================================================================
00072 *
00073 *     .. Local Scalars ..
00074       INTEGER            NCOLS, I, J, K, KP
00075       DOUBLE PRECISION   AMAX, UMAX, RPVGRW, TMP
00076       LOGICAL            UPPER
00077       COMPLEX*16         ZDUM
00078 *     ..
00079 *     .. Intrinsic Functions ..
00080       INTRINSIC          ABS, REAL, DIMAG, MAX, MIN
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           LSAME, ZLASET
00084       LOGICAL            LSAME
00085 *     ..
00086 *     .. Statement Functions ..
00087       DOUBLE PRECISION   CABS1
00088 *     ..
00089 *     .. Statement Function Definitions ..
00090       CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
00091 *     ..
00092 *     .. Executable Statements ..
00093 *
00094       UPPER = LSAME( 'Upper', UPLO )
00095       IF ( INFO.EQ.0 ) THEN
00096          IF ( UPPER ) THEN
00097             NCOLS = 1
00098          ELSE
00099             NCOLS = N
00100          END IF
00101       ELSE
00102          NCOLS = INFO
00103       END IF
00104 
00105       RPVGRW = 1.0D+0
00106       DO I = 1, 2*N
00107          WORK( I ) = 0.0D+0
00108       END DO
00109 *
00110 *     Find the max magnitude entry of each column of A.  Compute the max
00111 *     for all N columns so we can apply the pivot permutation while
00112 *     looping below.  Assume a full factorization is the common case.
00113 *
00114       IF ( UPPER ) THEN
00115          DO J = 1, N
00116             DO I = 1, J
00117                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
00118                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
00119             END DO
00120          END DO
00121       ELSE
00122          DO J = 1, N
00123             DO I = J, N
00124                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
00125                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
00126             END DO
00127          END DO
00128       END IF
00129 *
00130 *     Now find the max magnitude entry of each column of U or L.  Also
00131 *     permute the magnitudes of A above so they're in the same order as
00132 *     the factor.
00133 *
00134 *     The iteration orders and permutations were copied from zsytrs.
00135 *     Calls to SSWAP would be severe overkill.
00136 *
00137       IF ( UPPER ) THEN
00138          K = N
00139          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
00140             IF ( IPIV( K ).GT.0 ) THEN
00141 !              1x1 pivot
00142                KP = IPIV( K )
00143                IF ( KP .NE. K ) THEN
00144                   TMP = WORK( N+K )
00145                   WORK( N+K ) = WORK( N+KP )
00146                   WORK( N+KP ) = TMP
00147                END IF
00148                DO I = 1, K
00149                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00150                END DO
00151                K = K - 1
00152             ELSE
00153 !              2x2 pivot
00154                KP = -IPIV( K )
00155                TMP = WORK( N+K-1 )
00156                WORK( N+K-1 ) = WORK( N+KP )
00157                WORK( N+KP ) = TMP
00158                DO I = 1, K-1
00159                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00160                   WORK( K-1 ) =
00161      $                 MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
00162                END DO
00163                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00164                K = K - 2
00165             END IF
00166          END DO
00167          K = NCOLS
00168          DO WHILE ( K .LE. N )
00169             IF ( IPIV( K ).GT.0 ) THEN
00170                KP = IPIV( K )
00171                IF ( KP .NE. K ) THEN
00172                   TMP = WORK( N+K )
00173                   WORK( N+K ) = WORK( N+KP )
00174                   WORK( N+KP ) = TMP
00175                END IF
00176                K = K + 1
00177             ELSE
00178                KP = -IPIV( K )
00179                TMP = WORK( N+K )
00180                WORK( N+K ) = WORK( N+KP )
00181                WORK( N+KP ) = TMP
00182                K = K + 2
00183             END IF
00184          END DO
00185       ELSE
00186          K = 1
00187          DO WHILE ( K .LE. NCOLS )
00188             IF ( IPIV( K ).GT.0 ) THEN
00189 !              1x1 pivot
00190                KP = IPIV( K )
00191                IF ( KP .NE. K ) THEN
00192                   TMP = WORK( N+K )
00193                   WORK( N+K ) = WORK( N+KP )
00194                   WORK( N+KP ) = TMP
00195                END IF
00196                DO I = K, N
00197                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00198                END DO
00199                K = K + 1
00200             ELSE
00201 !              2x2 pivot
00202                KP = -IPIV( K )
00203                TMP = WORK( N+K+1 )
00204                WORK( N+K+1 ) = WORK( N+KP )
00205                WORK( N+KP ) = TMP
00206                DO I = K+1, N
00207                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00208                   WORK( K+1 ) =
00209      $                 MAX( CABS1( AF( I, K+1 ) ), WORK( K+1 ) )
00210                END DO
00211                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00212                K = K + 2
00213             END IF
00214          END DO
00215          K = NCOLS
00216          DO WHILE ( K .GE. 1 )
00217             IF ( IPIV( K ).GT.0 ) THEN
00218                KP = IPIV( K )
00219                IF ( KP .NE. K ) THEN
00220                   TMP = WORK( N+K )
00221                   WORK( N+K ) = WORK( N+KP )
00222                   WORK( N+KP ) = TMP
00223                END IF
00224                K = K - 1
00225             ELSE
00226                KP = -IPIV( K )
00227                TMP = WORK( N+K )
00228                WORK( N+K ) = WORK( N+KP )
00229                WORK( N+KP ) = TMP
00230                K = K - 2
00231             ENDIF
00232          END DO
00233       END IF
00234 *
00235 *     Compute the *inverse* of the max element growth factor.  Dividing
00236 *     by zero would imply the largest entry of the factor's column is
00237 *     zero.  Than can happen when either the column of A is zero or
00238 *     massive pivots made the factor underflow to zero.  Neither counts
00239 *     as growth in itself, so simply ignore terms with zero
00240 *     denominators.
00241 *
00242       IF ( UPPER ) THEN
00243          DO I = NCOLS, N
00244             UMAX = WORK( I )
00245             AMAX = WORK( N+I )
00246             IF ( UMAX /= 0.0D+0 ) THEN
00247                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00248             END IF
00249          END DO
00250       ELSE
00251          DO I = 1, NCOLS
00252             UMAX = WORK( I )
00253             AMAX = WORK( N+I )
00254             IF ( UMAX /= 0.0D+0 ) THEN
00255                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00256             END IF
00257          END DO
00258       END IF
00259 
00260       ZLA_SYRPVGRW = RPVGRW
00261       END