LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sla_gbrcond.f
Go to the documentation of this file.
00001 *> \brief \b SLA_GBRCOND
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLA_GBRCOND + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gbrcond.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gbrcond.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gbrcond.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       REAL FUNCTION SLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB,
00022 *                                  IPIV, CMODE, C, INFO, WORK, IWORK )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          TRANS
00026 *       INTEGER            N, LDAB, LDAFB, INFO, KL, KU, CMODE
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * ), IPIV( * )
00030 *       REAL               AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
00031 *      $                   C( * )
00032 *      ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *>    SLA_GBRCOND Estimates the Skeel condition number of  op(A) * op2(C)
00041 *>    where op2 is determined by CMODE as follows
00042 *>    CMODE =  1    op2(C) = C
00043 *>    CMODE =  0    op2(C) = I
00044 *>    CMODE = -1    op2(C) = inv(C)
00045 *>    The Skeel condition number  cond(A) = norminf( |inv(A)||A| )
00046 *>    is computed by computing scaling factors R such that
00047 *>    diag(R)*A*op2(C) is row equilibrated and computing the standard
00048 *>    infinity-norm condition number.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] TRANS
00055 *> \verbatim
00056 *>          TRANS is CHARACTER*1
00057 *>     Specifies the form of the system of equations:
00058 *>       = 'N':  A * X = B     (No transpose)
00059 *>       = 'T':  A**T * X = B  (Transpose)
00060 *>       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
00061 *> \endverbatim
00062 *>
00063 *> \param[in] N
00064 *> \verbatim
00065 *>          N is INTEGER
00066 *>     The number of linear equations, i.e., the order of the
00067 *>     matrix A.  N >= 0.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] KL
00071 *> \verbatim
00072 *>          KL is INTEGER
00073 *>     The number of subdiagonals within the band of A.  KL >= 0.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] KU
00077 *> \verbatim
00078 *>          KU is INTEGER
00079 *>     The number of superdiagonals within the band of A.  KU >= 0.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] AB
00083 *> \verbatim
00084 *>          AB is REAL array, dimension (LDAB,N)
00085 *>     On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
00086 *>     The j-th column of A is stored in the j-th column of the
00087 *>     array AB as follows:
00088 *>     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDAB
00092 *> \verbatim
00093 *>          LDAB is INTEGER
00094 *>     The leading dimension of the array AB.  LDAB >= KL+KU+1.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] AFB
00098 *> \verbatim
00099 *>          AFB is REAL array, dimension (LDAFB,N)
00100 *>     Details of the LU factorization of the band matrix A, as
00101 *>     computed by SGBTRF.  U is stored as an upper triangular
00102 *>     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
00103 *>     and the multipliers used during the factorization are stored
00104 *>     in rows KL+KU+2 to 2*KL+KU+1.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] LDAFB
00108 *> \verbatim
00109 *>          LDAFB is INTEGER
00110 *>     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1.
00111 *> \endverbatim
00112 *>
00113 *> \param[in] IPIV
00114 *> \verbatim
00115 *>          IPIV is INTEGER array, dimension (N)
00116 *>     The pivot indices from the factorization A = P*L*U
00117 *>     as computed by SGBTRF; row i of the matrix was interchanged
00118 *>     with row IPIV(i).
00119 *> \endverbatim
00120 *>
00121 *> \param[in] CMODE
00122 *> \verbatim
00123 *>          CMODE is INTEGER
00124 *>     Determines op2(C) in the formula op(A) * op2(C) as follows:
00125 *>     CMODE =  1    op2(C) = C
00126 *>     CMODE =  0    op2(C) = I
00127 *>     CMODE = -1    op2(C) = inv(C)
00128 *> \endverbatim
00129 *>
00130 *> \param[in] C
00131 *> \verbatim
00132 *>          C is REAL array, dimension (N)
00133 *>     The vector C in the formula op(A) * op2(C).
00134 *> \endverbatim
00135 *>
00136 *> \param[out] INFO
00137 *> \verbatim
00138 *>          INFO is INTEGER
00139 *>       = 0:  Successful exit.
00140 *>     i > 0:  The ith argument is invalid.
00141 *> \endverbatim
00142 *>
00143 *> \param[in] WORK
00144 *> \verbatim
00145 *>          WORK is REAL array, dimension (5*N).
00146 *>     Workspace.
00147 *> \endverbatim
00148 *>
00149 *> \param[in] IWORK
00150 *> \verbatim
00151 *>          IWORK is INTEGER array, dimension (N).
00152 *>     Workspace.
00153 *> \endverbatim
00154 *
00155 *  Authors:
00156 *  ========
00157 *
00158 *> \author Univ. of Tennessee 
00159 *> \author Univ. of California Berkeley 
00160 *> \author Univ. of Colorado Denver 
00161 *> \author NAG Ltd. 
00162 *
00163 *> \date November 2011
00164 *
00165 *> \ingroup realGBcomputational
00166 *
00167 *  =====================================================================
00168       REAL FUNCTION SLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB,
00169      $                           IPIV, CMODE, C, INFO, WORK, IWORK )
00170 *
00171 *  -- LAPACK computational routine (version 3.4.0) --
00172 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00173 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00174 *     November 2011
00175 *
00176 *     .. Scalar Arguments ..
00177       CHARACTER          TRANS
00178       INTEGER            N, LDAB, LDAFB, INFO, KL, KU, CMODE
00179 *     ..
00180 *     .. Array Arguments ..
00181       INTEGER            IWORK( * ), IPIV( * )
00182       REAL               AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
00183      $                   C( * )
00184 *    ..
00185 *
00186 *  =====================================================================
00187 *
00188 *     .. Local Scalars ..
00189       LOGICAL            NOTRANS
00190       INTEGER            KASE, I, J, KD, KE
00191       REAL               AINVNM, TMP
00192 *     ..
00193 *     .. Local Arrays ..
00194       INTEGER            ISAVE( 3 )
00195 *     ..
00196 *     .. External Functions ..
00197       LOGICAL            LSAME
00198       EXTERNAL           LSAME
00199 *     ..
00200 *     .. External Subroutines ..
00201       EXTERNAL           SLACN2, SGBTRS, XERBLA
00202 *     ..
00203 *     .. Intrinsic Functions ..
00204       INTRINSIC          ABS, MAX
00205 *     ..
00206 *     .. Executable Statements ..
00207 *
00208       SLA_GBRCOND = 0.0
00209 *
00210       INFO = 0
00211       NOTRANS = LSAME( TRANS, 'N' )
00212       IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
00213      $     .AND. .NOT. LSAME(TRANS, 'C') ) THEN
00214          INFO = -1
00215       ELSE IF( N.LT.0 ) THEN
00216          INFO = -2
00217       ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
00218          INFO = -3
00219       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
00220          INFO = -4
00221       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
00222          INFO = -6
00223       ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
00224          INFO = -8
00225       END IF
00226       IF( INFO.NE.0 ) THEN
00227          CALL XERBLA( 'SLA_GBRCOND', -INFO )
00228          RETURN
00229       END IF
00230       IF( N.EQ.0 ) THEN
00231          SLA_GBRCOND = 1.0
00232          RETURN
00233       END IF
00234 *
00235 *     Compute the equilibration matrix R such that
00236 *     inv(R)*A*C has unit 1-norm.
00237 *
00238       KD = KU + 1
00239       KE = KL + 1
00240       IF ( NOTRANS ) THEN
00241          DO I = 1, N
00242             TMP = 0.0
00243                IF ( CMODE .EQ. 1 ) THEN
00244                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00245                   TMP = TMP + ABS( AB( KD+I-J, J ) * C( J ) )
00246                END DO
00247                ELSE IF ( CMODE .EQ. 0 ) THEN
00248                   DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00249                      TMP = TMP + ABS( AB( KD+I-J, J ) )
00250                   END DO
00251                ELSE
00252                   DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00253                      TMP = TMP + ABS( AB( KD+I-J, J ) / C( J ) )
00254                   END DO
00255                END IF
00256             WORK( 2*N+I ) = TMP
00257          END DO
00258       ELSE
00259          DO I = 1, N
00260             TMP = 0.0
00261             IF ( CMODE .EQ. 1 ) THEN
00262                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00263                   TMP = TMP + ABS( AB( KE-I+J, I ) * C( J ) )
00264                END DO
00265             ELSE IF ( CMODE .EQ. 0 ) THEN
00266                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00267                   TMP = TMP + ABS( AB( KE-I+J, I ) )
00268                END DO
00269             ELSE
00270                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
00271                   TMP = TMP + ABS( AB( KE-I+J, I ) / C( J ) )
00272                END DO
00273             END IF
00274             WORK( 2*N+I ) = TMP
00275          END DO
00276       END IF
00277 *
00278 *     Estimate the norm of inv(op(A)).
00279 *
00280       AINVNM = 0.0
00281 
00282       KASE = 0
00283    10 CONTINUE
00284       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00285       IF( KASE.NE.0 ) THEN
00286          IF( KASE.EQ.2 ) THEN
00287 *
00288 *           Multiply by R.
00289 *
00290             DO I = 1, N
00291                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00292             END DO
00293 
00294             IF ( NOTRANS ) THEN
00295                CALL SGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
00296      $              IPIV, WORK, N, INFO )
00297             ELSE
00298                CALL SGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV,
00299      $              WORK, N, INFO )
00300             END IF
00301 *
00302 *           Multiply by inv(C).
00303 *
00304             IF ( CMODE .EQ. 1 ) THEN
00305                DO I = 1, N
00306                   WORK( I ) = WORK( I ) / C( I )
00307                END DO
00308             ELSE IF ( CMODE .EQ. -1 ) THEN
00309                DO I = 1, N
00310                   WORK( I ) = WORK( I ) * C( I )
00311                END DO
00312             END IF
00313          ELSE
00314 *
00315 *           Multiply by inv(C**T).
00316 *
00317             IF ( CMODE .EQ. 1 ) THEN
00318                DO I = 1, N
00319                   WORK( I ) = WORK( I ) / C( I )
00320                END DO
00321             ELSE IF ( CMODE .EQ. -1 ) THEN
00322                DO I = 1, N
00323                   WORK( I ) = WORK( I ) * C( I )
00324                END DO
00325             END IF
00326 
00327             IF ( NOTRANS ) THEN
00328                CALL SGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV,
00329      $              WORK, N, INFO )
00330             ELSE
00331                CALL SGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
00332      $              IPIV, WORK, N, INFO )
00333             END IF
00334 *
00335 *           Multiply by R.
00336 *
00337             DO I = 1, N
00338                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00339             END DO
00340          END IF
00341          GO TO 10
00342       END IF
00343 *
00344 *     Compute the estimate of the reciprocal condition number.
00345 *
00346       IF( AINVNM .NE. 0.0 )
00347      $   SLA_GBRCOND = ( 1.0 / AINVNM )
00348 *
00349       RETURN
00350 *
00351       END
 All Files Functions