LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dla_syrpvgrw.f
Go to the documentation of this file.
00001 *> \brief \b DLA_SYRPVGRW
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLA_SYRPVGRW + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
00022 *                                               LDAF, IPIV, WORK )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER*1        UPLO
00026 *       INTEGER            N, INFO, LDA, LDAF
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IPIV( * )
00030 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> 
00040 *> DLA_SYRPVGRW computes the reciprocal pivot growth factor
00041 *> norm(A)/norm(U). The "max absolute element" norm is used. If this is
00042 *> much less than 1, the stability of the LU factorization of the
00043 *> (equilibrated) matrix A could be poor. This also means that the
00044 *> solution X, estimated condition numbers, and error bounds could be
00045 *> unreliable.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] UPLO
00052 *> \verbatim
00053 *>          UPLO is CHARACTER*1
00054 *>       = 'U':  Upper triangle of A is stored;
00055 *>       = 'L':  Lower triangle of A is stored.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] N
00059 *> \verbatim
00060 *>          N is INTEGER
00061 *>     The number of linear equations, i.e., the order of the
00062 *>     matrix A.  N >= 0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] INFO
00066 *> \verbatim
00067 *>          INFO is INTEGER
00068 *>     The value of INFO returned from DSYTRF, .i.e., the pivot in
00069 *>     column INFO is exactly 0.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] A
00073 *> \verbatim
00074 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00075 *>     On entry, the N-by-N matrix A.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDA
00079 *> \verbatim
00080 *>          LDA is INTEGER
00081 *>     The leading dimension of the array A.  LDA >= max(1,N).
00082 *> \endverbatim
00083 *>
00084 *> \param[in] AF
00085 *> \verbatim
00086 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
00087 *>     The block diagonal matrix D and the multipliers used to
00088 *>     obtain the factor U or L as computed by DSYTRF.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDAF
00092 *> \verbatim
00093 *>          LDAF is INTEGER
00094 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
00095 *> \endverbatim
00096 *>
00097 *> \param[in] IPIV
00098 *> \verbatim
00099 *>          IPIV is INTEGER array, dimension (N)
00100 *>     Details of the interchanges and the block structure of D
00101 *>     as determined by DSYTRF.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] WORK
00105 *> \verbatim
00106 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
00107 *> \endverbatim
00108 *
00109 *  Authors:
00110 *  ========
00111 *
00112 *> \author Univ. of Tennessee 
00113 *> \author Univ. of California Berkeley 
00114 *> \author Univ. of Colorado Denver 
00115 *> \author NAG Ltd. 
00116 *
00117 *> \date November 2011
00118 *
00119 *> \ingroup doubleSYcomputational
00120 *
00121 *  =====================================================================
00122       DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
00123      $                                        LDAF, IPIV, WORK )
00124 *
00125 *  -- LAPACK computational routine (version 3.4.0) --
00126 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00127 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00128 *     November 2011
00129 *
00130 *     .. Scalar Arguments ..
00131       CHARACTER*1        UPLO
00132       INTEGER            N, INFO, LDA, LDAF
00133 *     ..
00134 *     .. Array Arguments ..
00135       INTEGER            IPIV( * )
00136       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
00137 *     ..
00138 *
00139 *  =====================================================================
00140 *
00141 *     .. Local Scalars ..
00142       INTEGER            NCOLS, I, J, K, KP
00143       DOUBLE PRECISION   AMAX, UMAX, RPVGRW, TMP
00144       LOGICAL            UPPER
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          ABS, MAX, MIN
00148 *     ..
00149 *     .. External Functions ..
00150       EXTERNAL           LSAME, DLASET
00151       LOGICAL            LSAME
00152 *     ..
00153 *     .. Executable Statements ..
00154 *
00155       UPPER = LSAME( 'Upper', UPLO )
00156       IF ( INFO.EQ.0 ) THEN
00157          IF ( UPPER ) THEN
00158             NCOLS = 1
00159          ELSE
00160             NCOLS = N
00161          END IF
00162       ELSE
00163          NCOLS = INFO
00164       END IF
00165 
00166       RPVGRW = 1.0D+0
00167       DO I = 1, 2*N
00168          WORK( I ) = 0.0D+0
00169       END DO
00170 *
00171 *     Find the max magnitude entry of each column of A.  Compute the max
00172 *     for all N columns so we can apply the pivot permutation while
00173 *     looping below.  Assume a full factorization is the common case.
00174 *
00175       IF ( UPPER ) THEN
00176          DO J = 1, N
00177             DO I = 1, J
00178                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
00179                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
00180             END DO
00181          END DO
00182       ELSE
00183          DO J = 1, N
00184             DO I = J, N
00185                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
00186                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
00187             END DO
00188          END DO
00189       END IF
00190 *
00191 *     Now find the max magnitude entry of each column of U or L.  Also
00192 *     permute the magnitudes of A above so they're in the same order as
00193 *     the factor.
00194 *
00195 *     The iteration orders and permutations were copied from dsytrs.
00196 *     Calls to SSWAP would be severe overkill.
00197 *
00198       IF ( UPPER ) THEN
00199          K = N
00200          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
00201             IF ( IPIV( K ).GT.0 ) THEN
00202 !              1x1 pivot
00203                KP = IPIV( K )
00204                IF ( KP .NE. K ) THEN
00205                   TMP = WORK( N+K )
00206                   WORK( N+K ) = WORK( N+KP )
00207                   WORK( N+KP ) = TMP
00208                END IF
00209                DO I = 1, K
00210                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
00211                END DO
00212                K = K - 1
00213             ELSE
00214 !              2x2 pivot
00215                KP = -IPIV( K )
00216                TMP = WORK( N+K-1 )
00217                WORK( N+K-1 ) = WORK( N+KP )
00218                WORK( N+KP ) = TMP
00219                DO I = 1, K-1
00220                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
00221                   WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
00222                END DO
00223                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
00224                K = K - 2
00225             END IF
00226          END DO
00227          K = NCOLS
00228          DO WHILE ( K .LE. N )
00229             IF ( IPIV( K ).GT.0 ) THEN
00230                KP = IPIV( K )
00231                IF ( KP .NE. K ) THEN
00232                   TMP = WORK( N+K )
00233                   WORK( N+K ) = WORK( N+KP )
00234                   WORK( N+KP ) = TMP
00235                END IF
00236                K = K + 1
00237             ELSE
00238                KP = -IPIV( K )
00239                TMP = WORK( N+K )
00240                WORK( N+K ) = WORK( N+KP )
00241                WORK( N+KP ) = TMP
00242                K = K + 2
00243             END IF
00244          END DO
00245       ELSE
00246          K = 1
00247          DO WHILE ( K .LE. NCOLS )
00248             IF ( IPIV( K ).GT.0 ) THEN
00249 !              1x1 pivot
00250                KP = IPIV( K )
00251                IF ( KP .NE. K ) THEN
00252                   TMP = WORK( N+K )
00253                   WORK( N+K ) = WORK( N+KP )
00254                   WORK( N+KP ) = TMP
00255                END IF
00256                DO I = K, N
00257                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
00258                END DO
00259                K = K + 1
00260             ELSE
00261 !              2x2 pivot
00262                KP = -IPIV( K )
00263                TMP = WORK( N+K+1 )
00264                WORK( N+K+1 ) = WORK( N+KP )
00265                WORK( N+KP ) = TMP
00266                DO I = K+1, N
00267                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
00268                   WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
00269                END DO
00270                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
00271                K = K + 2
00272             END IF
00273          END DO
00274          K = NCOLS
00275          DO WHILE ( K .GE. 1 )
00276             IF ( IPIV( K ).GT.0 ) THEN
00277                KP = IPIV( K )
00278                IF ( KP .NE. K ) THEN
00279                   TMP = WORK( N+K )
00280                   WORK( N+K ) = WORK( N+KP )
00281                   WORK( N+KP ) = TMP
00282                END IF
00283                K = K - 1
00284             ELSE
00285                KP = -IPIV( K )
00286                TMP = WORK( N+K )
00287                WORK( N+K ) = WORK( N+KP )
00288                WORK( N+KP ) = TMP
00289                K = K - 2
00290             ENDIF
00291          END DO
00292       END IF
00293 *
00294 *     Compute the *inverse* of the max element growth factor.  Dividing
00295 *     by zero would imply the largest entry of the factor's column is
00296 *     zero.  Than can happen when either the column of A is zero or
00297 *     massive pivots made the factor underflow to zero.  Neither counts
00298 *     as growth in itself, so simply ignore terms with zero
00299 *     denominators.
00300 *
00301       IF ( UPPER ) THEN
00302          DO I = NCOLS, N
00303             UMAX = WORK( I )
00304             AMAX = WORK( N+I )
00305             IF ( UMAX /= 0.0D+0 ) THEN
00306                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00307             END IF
00308          END DO
00309       ELSE
00310          DO I = 1, NCOLS
00311             UMAX = WORK( I )
00312             AMAX = WORK( N+I )
00313             IF ( UMAX /= 0.0D+0 ) THEN
00314                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00315             END IF
00316          END DO
00317       END IF
00318 
00319       DLA_SYRPVGRW = RPVGRW
00320       END
 All Files Functions