LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cla_syrpvgrw.f
Go to the documentation of this file.
00001 *> \brief \b CLA_SYRPVGRW
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLA_SYRPVGRW + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_syrpvgrw.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syrpvgrw.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syrpvgrw.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       REAL FUNCTION CLA_SYRPVGRW( 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 *       COMPLEX            A( LDA, * ), AF( LDAF, * )
00030 *       REAL               WORK( * )
00031 *       INTEGER            IPIV( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> 
00041 *> CLA_SYRPVGRW 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 CSYTRF, .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 CSYTRF.
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 CSYTRF.
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 complexSYcomputational
00121 *
00122 *  =====================================================================
00123       REAL FUNCTION CLA_SYRPVGRW( 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       COMPLEX            A( LDA, * ), AF( LDAF, * )
00137       REAL               WORK( * )
00138       INTEGER            IPIV( * )
00139 *     ..
00140 *
00141 *  =====================================================================
00142 *
00143 *     .. Local Scalars ..
00144       INTEGER            NCOLS, I, J, K, KP
00145       REAL               AMAX, UMAX, RPVGRW, TMP
00146       LOGICAL            UPPER
00147       COMPLEX            ZDUM
00148 *     ..
00149 *     .. Intrinsic Functions ..
00150       INTRINSIC          ABS, REAL, AIMAG, MAX, MIN
00151 *     ..
00152 *     .. External Subroutines ..
00153       EXTERNAL           LSAME, CLASET
00154       LOGICAL            LSAME
00155 *     ..
00156 *     .. Statement Functions ..
00157       REAL               CABS1
00158 *     ..
00159 *     .. Statement Function Definitions ..
00160       CABS1( ZDUM ) = ABS( REAL ( ZDUM ) ) + ABS( AIMAG ( ZDUM ) )
00161 *     ..
00162 *     .. Executable Statements ..
00163 *
00164       UPPER = LSAME( 'Upper', UPLO )
00165       IF ( INFO.EQ.0 ) THEN
00166          IF ( UPPER ) THEN
00167             NCOLS = 1
00168          ELSE
00169             NCOLS = N
00170          END IF
00171       ELSE
00172          NCOLS = INFO
00173       END IF
00174 
00175       RPVGRW = 1.0
00176       DO I = 1, 2*N
00177          WORK( I ) = 0.0
00178       END DO
00179 *
00180 *     Find the max magnitude entry of each column of A.  Compute the max
00181 *     for all N columns so we can apply the pivot permutation while
00182 *     looping below.  Assume a full factorization is the common case.
00183 *
00184       IF ( UPPER ) THEN
00185          DO J = 1, N
00186             DO I = 1, J
00187                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
00188                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
00189             END DO
00190          END DO
00191       ELSE
00192          DO J = 1, N
00193             DO I = J, N
00194                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
00195                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
00196             END DO
00197          END DO
00198       END IF
00199 *
00200 *     Now find the max magnitude entry of each column of U or L.  Also
00201 *     permute the magnitudes of A above so they're in the same order as
00202 *     the factor.
00203 *
00204 *     The iteration orders and permutations were copied from csytrs.
00205 *     Calls to SSWAP would be severe overkill.
00206 *
00207       IF ( UPPER ) THEN
00208          K = N
00209          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
00210             IF ( IPIV( K ).GT.0 ) THEN
00211 !              1x1 pivot
00212                KP = IPIV( K )
00213                IF ( KP .NE. K ) THEN
00214                   TMP = WORK( N+K )
00215                   WORK( N+K ) = WORK( N+KP )
00216                   WORK( N+KP ) = TMP
00217                END IF
00218                DO I = 1, K
00219                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00220                END DO
00221                K = K - 1
00222             ELSE
00223 !              2x2 pivot
00224                KP = -IPIV( K )
00225                TMP = WORK( N+K-1 )
00226                WORK( N+K-1 ) = WORK( N+KP )
00227                WORK( N+KP ) = TMP
00228                DO I = 1, K-1
00229                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00230                   WORK( K-1 ) =
00231      $                 MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
00232                END DO
00233                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00234                K = K - 2
00235             END IF
00236          END DO
00237          K = NCOLS
00238          DO WHILE ( K .LE. N )
00239             IF ( IPIV( K ).GT.0 ) THEN
00240                KP = IPIV( K )
00241                IF ( KP .NE. K ) THEN
00242                   TMP = WORK( N+K )
00243                   WORK( N+K ) = WORK( N+KP )
00244                   WORK( N+KP ) = TMP
00245                END IF
00246                K = K + 1
00247             ELSE
00248                KP = -IPIV( K )
00249                TMP = WORK( N+K )
00250                WORK( N+K ) = WORK( N+KP )
00251                WORK( N+KP ) = TMP
00252                K = K + 2
00253             END IF
00254          END DO
00255       ELSE
00256          K = 1
00257          DO WHILE ( K .LE. NCOLS )
00258             IF ( IPIV( K ).GT.0 ) THEN
00259 !              1x1 pivot
00260                KP = IPIV( K )
00261                IF ( KP .NE. K ) THEN
00262                   TMP = WORK( N+K )
00263                   WORK( N+K ) = WORK( N+KP )
00264                   WORK( N+KP ) = TMP
00265                END IF
00266                DO I = K, N
00267                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00268                END DO
00269                K = K + 1
00270             ELSE
00271 !              2x2 pivot
00272                KP = -IPIV( K )
00273                TMP = WORK( N+K+1 )
00274                WORK( N+K+1 ) = WORK( N+KP )
00275                WORK( N+KP ) = TMP
00276                DO I = K+1, N
00277                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
00278                   WORK( K+1 ) =
00279      $                 MAX( CABS1( AF( I, K+1 ) ), WORK( K+1 ) )
00280                END DO
00281                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
00282                K = K + 2
00283             END IF
00284          END DO
00285          K = NCOLS
00286          DO WHILE ( K .GE. 1 )
00287             IF ( IPIV( K ).GT.0 ) THEN
00288                KP = IPIV( K )
00289                IF ( KP .NE. K ) THEN
00290                   TMP = WORK( N+K )
00291                   WORK( N+K ) = WORK( N+KP )
00292                   WORK( N+KP ) = TMP
00293                END IF
00294                K = K - 1
00295             ELSE
00296                KP = -IPIV( K )
00297                TMP = WORK( N+K )
00298                WORK( N+K ) = WORK( N+KP )
00299                WORK( N+KP ) = TMP
00300                K = K - 2
00301             ENDIF
00302          END DO
00303       END IF
00304 *
00305 *     Compute the *inverse* of the max element growth factor.  Dividing
00306 *     by zero would imply the largest entry of the factor's column is
00307 *     zero.  Than can happen when either the column of A is zero or
00308 *     massive pivots made the factor underflow to zero.  Neither counts
00309 *     as growth in itself, so simply ignore terms with zero
00310 *     denominators.
00311 *
00312       IF ( UPPER ) THEN
00313          DO I = NCOLS, N
00314             UMAX = WORK( I )
00315             AMAX = WORK( N+I )
00316             IF ( UMAX /= 0.0 ) THEN
00317                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00318             END IF
00319          END DO
00320       ELSE
00321          DO I = 1, NCOLS
00322             UMAX = WORK( I )
00323             AMAX = WORK( N+I )
00324             IF ( UMAX /= 0.0 ) THEN
00325                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
00326             END IF
00327          END DO
00328       END IF
00329 
00330       CLA_SYRPVGRW = RPVGRW
00331       END
 All Files Functions