LAPACK: SRC/sgbcon.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00002      $                   WORK, IWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH.
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          NORM
00013       INTEGER            INFO, KL, KU, LDAB, N
00014       REAL               ANORM, RCOND
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IPIV( * ), IWORK( * )
00018       REAL               AB( LDAB, * ), WORK( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  SGBCON estimates the reciprocal of the condition number of a real
00025 *  general band matrix A, in either the 1-norm or the infinity-norm,
00026 *  using the LU factorization computed by SGBTRF.
00027 *
00028 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
00029 *  condition number is computed as
00030 *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  NORM    (input) CHARACTER*1
00036 *          Specifies whether the 1-norm condition number or the
00037 *          infinity-norm condition number is required:
00038 *          = '1' or 'O':  1-norm;
00039 *          = 'I':         Infinity-norm.
00040 *
00041 *  N       (input) INTEGER
00042 *          The order of the matrix A.  N >= 0.
00043 *
00044 *  KL      (input) INTEGER
00045 *          The number of subdiagonals within the band of A.  KL >= 0.
00046 *
00047 *  KU      (input) INTEGER
00048 *          The number of superdiagonals within the band of A.  KU >= 0.
00049 *
00050 *  AB      (input) REAL array, dimension (LDAB,N)
00051 *          Details of the LU factorization of the band matrix A, as
00052 *          computed by SGBTRF.  U is stored as an upper triangular band
00053 *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
00054 *          the multipliers used during the factorization are stored in
00055 *          rows KL+KU+2 to 2*KL+KU+1.
00056 *
00057 *  LDAB    (input) INTEGER
00058 *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00059 *
00060 *  IPIV    (input) INTEGER array, dimension (N)
00061 *          The pivot indices; for 1 <= i <= N, row i of the matrix was
00062 *          interchanged with row IPIV(i).
00063 *
00064 *  ANORM   (input) REAL
00065 *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
00066 *          If NORM = 'I', the infinity-norm of the original matrix A.
00067 *
00068 *  RCOND   (output) REAL
00069 *          The reciprocal of the condition number of the matrix A,
00070 *          computed as RCOND = 1/(norm(A) * norm(inv(A))).
00071 *
00072 *  WORK    (workspace) REAL array, dimension (3*N)
00073 *
00074 *  IWORK   (workspace) INTEGER array, dimension (N)
00075 *
00076 *  INFO    (output) INTEGER
00077 *          = 0:  successful exit
00078 *          < 0: if INFO = -i, the i-th argument had an illegal value
00079 *
00080 *  =====================================================================
00081 *
00082 *     .. Parameters ..
00083       REAL               ONE, ZERO
00084       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00085 *     ..
00086 *     .. Local Scalars ..
00087       LOGICAL            LNOTI, ONENRM
00088       CHARACTER          NORMIN
00089       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
00090       REAL               AINVNM, SCALE, SMLNUM, T
00091 *     ..
00092 *     .. Local Arrays ..
00093       INTEGER            ISAVE( 3 )
00094 *     ..
00095 *     .. External Functions ..
00096       LOGICAL            LSAME
00097       INTEGER            ISAMAX
00098       REAL               SDOT, SLAMCH
00099       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           SAXPY, SLACN2, SLATBS, SRSCL, XERBLA
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          ABS, MIN
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109 *     Test the input parameters.
00110 *
00111       INFO = 0
00112       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
00113       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
00114          INFO = -1
00115       ELSE IF( N.LT.0 ) THEN
00116          INFO = -2
00117       ELSE IF( KL.LT.0 ) THEN
00118          INFO = -3
00119       ELSE IF( KU.LT.0 ) THEN
00120          INFO = -4
00121       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
00122          INFO = -6
00123       ELSE IF( ANORM.LT.ZERO ) THEN
00124          INFO = -8
00125       END IF
00126       IF( INFO.NE.0 ) THEN
00127          CALL XERBLA( 'SGBCON', -INFO )
00128          RETURN
00129       END IF
00130 *
00131 *     Quick return if possible
00132 *
00133       RCOND = ZERO
00134       IF( N.EQ.0 ) THEN
00135          RCOND = ONE
00136          RETURN
00137       ELSE IF( ANORM.EQ.ZERO ) THEN
00138          RETURN
00139       END IF
00140 *
00141       SMLNUM = SLAMCH( 'Safe minimum' )
00142 *
00143 *     Estimate the norm of inv(A).
00144 *
00145       AINVNM = ZERO
00146       NORMIN = 'N'
00147       IF( ONENRM ) THEN
00148          KASE1 = 1
00149       ELSE
00150          KASE1 = 2
00151       END IF
00152       KD = KL + KU + 1
00153       LNOTI = KL.GT.0
00154       KASE = 0
00155    10 CONTINUE
00156       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00157       IF( KASE.NE.0 ) THEN
00158          IF( KASE.EQ.KASE1 ) THEN
00159 *
00160 *           Multiply by inv(L).
00161 *
00162             IF( LNOTI ) THEN
00163                DO 20 J = 1, N - 1
00164                   LM = MIN( KL, N-J )
00165                   JP = IPIV( J )
00166                   T = WORK( JP )
00167                   IF( JP.NE.J ) THEN
00168                      WORK( JP ) = WORK( J )
00169                      WORK( J ) = T
00170                   END IF
00171                   CALL SAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
00172    20          CONTINUE
00173             END IF
00174 *
00175 *           Multiply by inv(U).
00176 *
00177             CALL SLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00178      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
00179      $                   INFO )
00180          ELSE
00181 *
00182 *           Multiply by inv(U').
00183 *
00184             CALL SLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
00185      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
00186      $                   INFO )
00187 *
00188 *           Multiply by inv(L').
00189 *
00190             IF( LNOTI ) THEN
00191                DO 30 J = N - 1, 1, -1
00192                   LM = MIN( KL, N-J )
00193                   WORK( J ) = WORK( J ) - SDOT( LM, AB( KD+1, J ), 1,
00194      $                        WORK( J+1 ), 1 )
00195                   JP = IPIV( J )
00196                   IF( JP.NE.J ) THEN
00197                      T = WORK( JP )
00198                      WORK( JP ) = WORK( J )
00199                      WORK( J ) = T
00200                   END IF
00201    30          CONTINUE
00202             END IF
00203          END IF
00204 *
00205 *        Divide X by 1/SCALE if doing so will not cause overflow.
00206 *
00207          NORMIN = 'Y'
00208          IF( SCALE.NE.ONE ) THEN
00209             IX = ISAMAX( N, WORK, 1 )
00210             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00211      $         GO TO 40
00212             CALL SRSCL( N, SCALE, WORK, 1 )
00213          END IF
00214          GO TO 10
00215       END IF
00216 *
00217 *     Compute the estimate of the reciprocal condition number.
00218 *
00219       IF( AINVNM.NE.ZERO )
00220      $   RCOND = ( ONE / AINVNM ) / ANORM
00221 *
00222    40 CONTINUE
00223       RETURN
00224 *
00225 *     End of SGBCON
00226 *
00227       END