LAPACK 3.3.0

dla_syrpvgrw.f

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