LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgbcon.f
Go to the documentation of this file.
00001 *> \brief \b SGBCON
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGBCON + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgbcon.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgbcon.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgbcon.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00022 *                          WORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          NORM
00026 *       INTEGER            INFO, KL, KU, LDAB, N
00027 *       REAL               ANORM, RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IPIV( * ), IWORK( * )
00031 *       REAL               AB( LDAB, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> SGBCON estimates the reciprocal of the condition number of a real
00041 *> general band matrix A, in either the 1-norm or the infinity-norm,
00042 *> using the LU factorization computed by SGBTRF.
00043 *>
00044 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
00045 *> condition number is computed as
00046 *>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] NORM
00053 *> \verbatim
00054 *>          NORM is CHARACTER*1
00055 *>          Specifies whether the 1-norm condition number or the
00056 *>          infinity-norm condition number is required:
00057 *>          = '1' or 'O':  1-norm;
00058 *>          = 'I':         Infinity-norm.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] N
00062 *> \verbatim
00063 *>          N is INTEGER
00064 *>          The order of the matrix A.  N >= 0.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] KL
00068 *> \verbatim
00069 *>          KL is INTEGER
00070 *>          The number of subdiagonals within the band of A.  KL >= 0.
00071 *> \endverbatim
00072 *>
00073 *> \param[in] KU
00074 *> \verbatim
00075 *>          KU is INTEGER
00076 *>          The number of superdiagonals within the band of A.  KU >= 0.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] AB
00080 *> \verbatim
00081 *>          AB is REAL array, dimension (LDAB,N)
00082 *>          Details of the LU factorization of the band matrix A, as
00083 *>          computed by SGBTRF.  U is stored as an upper triangular band
00084 *>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
00085 *>          the multipliers used during the factorization are stored in
00086 *>          rows KL+KU+2 to 2*KL+KU+1.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDAB
00090 *> \verbatim
00091 *>          LDAB is INTEGER
00092 *>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] IPIV
00096 *> \verbatim
00097 *>          IPIV is INTEGER array, dimension (N)
00098 *>          The pivot indices; for 1 <= i <= N, row i of the matrix was
00099 *>          interchanged with row IPIV(i).
00100 *> \endverbatim
00101 *>
00102 *> \param[in] ANORM
00103 *> \verbatim
00104 *>          ANORM is REAL
00105 *>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
00106 *>          If NORM = 'I', the infinity-norm of the original matrix A.
00107 *> \endverbatim
00108 *>
00109 *> \param[out] RCOND
00110 *> \verbatim
00111 *>          RCOND is REAL
00112 *>          The reciprocal of the condition number of the matrix A,
00113 *>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
00114 *> \endverbatim
00115 *>
00116 *> \param[out] WORK
00117 *> \verbatim
00118 *>          WORK is REAL array, dimension (3*N)
00119 *> \endverbatim
00120 *>
00121 *> \param[out] IWORK
00122 *> \verbatim
00123 *>          IWORK is INTEGER array, dimension (N)
00124 *> \endverbatim
00125 *>
00126 *> \param[out] INFO
00127 *> \verbatim
00128 *>          INFO is INTEGER
00129 *>          = 0:  successful exit
00130 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00131 *> \endverbatim
00132 *
00133 *  Authors:
00134 *  ========
00135 *
00136 *> \author Univ. of Tennessee 
00137 *> \author Univ. of California Berkeley 
00138 *> \author Univ. of Colorado Denver 
00139 *> \author NAG Ltd. 
00140 *
00141 *> \date November 2011
00142 *
00143 *> \ingroup realGBcomputational
00144 *
00145 *  =====================================================================
00146       SUBROUTINE SGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00147      $                   WORK, IWORK, INFO )
00148 *
00149 *  -- LAPACK computational routine (version 3.4.0) --
00150 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00151 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00152 *     November 2011
00153 *
00154 *     .. Scalar Arguments ..
00155       CHARACTER          NORM
00156       INTEGER            INFO, KL, KU, LDAB, N
00157       REAL               ANORM, RCOND
00158 *     ..
00159 *     .. Array Arguments ..
00160       INTEGER            IPIV( * ), IWORK( * )
00161       REAL               AB( LDAB, * ), WORK( * )
00162 *     ..
00163 *
00164 *  =====================================================================
00165 *
00166 *     .. Parameters ..
00167       REAL               ONE, ZERO
00168       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00169 *     ..
00170 *     .. Local Scalars ..
00171       LOGICAL            LNOTI, ONENRM
00172       CHARACTER          NORMIN
00173       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
00174       REAL               AINVNM, SCALE, SMLNUM, T
00175 *     ..
00176 *     .. Local Arrays ..
00177       INTEGER            ISAVE( 3 )
00178 *     ..
00179 *     .. External Functions ..
00180       LOGICAL            LSAME
00181       INTEGER            ISAMAX
00182       REAL               SDOT, SLAMCH
00183       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
00184 *     ..
00185 *     .. External Subroutines ..
00186       EXTERNAL           SAXPY, SLACN2, SLATBS, SRSCL, XERBLA
00187 *     ..
00188 *     .. Intrinsic Functions ..
00189       INTRINSIC          ABS, MIN
00190 *     ..
00191 *     .. Executable Statements ..
00192 *
00193 *     Test the input parameters.
00194 *
00195       INFO = 0
00196       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
00197       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
00198          INFO = -1
00199       ELSE IF( N.LT.0 ) THEN
00200          INFO = -2
00201       ELSE IF( KL.LT.0 ) THEN
00202          INFO = -3
00203       ELSE IF( KU.LT.0 ) THEN
00204          INFO = -4
00205       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
00206          INFO = -6
00207       ELSE IF( ANORM.LT.ZERO ) THEN
00208          INFO = -8
00209       END IF
00210       IF( INFO.NE.0 ) THEN
00211          CALL XERBLA( 'SGBCON', -INFO )
00212          RETURN
00213       END IF
00214 *
00215 *     Quick return if possible
00216 *
00217       RCOND = ZERO
00218       IF( N.EQ.0 ) THEN
00219          RCOND = ONE
00220          RETURN
00221       ELSE IF( ANORM.EQ.ZERO ) THEN
00222          RETURN
00223       END IF
00224 *
00225       SMLNUM = SLAMCH( 'Safe minimum' )
00226 *
00227 *     Estimate the norm of inv(A).
00228 *
00229       AINVNM = ZERO
00230       NORMIN = 'N'
00231       IF( ONENRM ) THEN
00232          KASE1 = 1
00233       ELSE
00234          KASE1 = 2
00235       END IF
00236       KD = KL + KU + 1
00237       LNOTI = KL.GT.0
00238       KASE = 0
00239    10 CONTINUE
00240       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00241       IF( KASE.NE.0 ) THEN
00242          IF( KASE.EQ.KASE1 ) THEN
00243 *
00244 *           Multiply by inv(L).
00245 *
00246             IF( LNOTI ) THEN
00247                DO 20 J = 1, N - 1
00248                   LM = MIN( KL, N-J )
00249                   JP = IPIV( J )
00250                   T = WORK( JP )
00251                   IF( JP.NE.J ) THEN
00252                      WORK( JP ) = WORK( J )
00253                      WORK( J ) = T
00254                   END IF
00255                   CALL SAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
00256    20          CONTINUE
00257             END IF
00258 *
00259 *           Multiply by inv(U).
00260 *
00261             CALL SLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00262      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
00263      $                   INFO )
00264          ELSE
00265 *
00266 *           Multiply by inv(U**T).
00267 *
00268             CALL SLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
00269      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
00270      $                   INFO )
00271 *
00272 *           Multiply by inv(L**T).
00273 *
00274             IF( LNOTI ) THEN
00275                DO 30 J = N - 1, 1, -1
00276                   LM = MIN( KL, N-J )
00277                   WORK( J ) = WORK( J ) - SDOT( LM, AB( KD+1, J ), 1,
00278      $                        WORK( J+1 ), 1 )
00279                   JP = IPIV( J )
00280                   IF( JP.NE.J ) THEN
00281                      T = WORK( JP )
00282                      WORK( JP ) = WORK( J )
00283                      WORK( J ) = T
00284                   END IF
00285    30          CONTINUE
00286             END IF
00287          END IF
00288 *
00289 *        Divide X by 1/SCALE if doing so will not cause overflow.
00290 *
00291          NORMIN = 'Y'
00292          IF( SCALE.NE.ONE ) THEN
00293             IX = ISAMAX( N, WORK, 1 )
00294             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00295      $         GO TO 40
00296             CALL SRSCL( N, SCALE, WORK, 1 )
00297          END IF
00298          GO TO 10
00299       END IF
00300 *
00301 *     Compute the estimate of the reciprocal condition number.
00302 *
00303       IF( AINVNM.NE.ZERO )
00304      $   RCOND = ( ONE / AINVNM ) / ANORM
00305 *
00306    40 CONTINUE
00307       RETURN
00308 *
00309 *     End of SGBCON
00310 *
00311       END
 All Files Functions