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