dgbcon_8f_source.html

*> \brief \b DGBCON

*

*  =========== DOCUMENTATION ===========

*

* Online html documentation available at

*            http://www.netlib.org/lapack/explore-html/

*

*> \htmlonly

*> Download DGBCON + dependencies

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbcon.f">

*> [TGZ]</a>

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbcon.f">

*> [ZIP]</a>

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbcon.f">

*> [TXT]</a>

*> \endhtmlonly

*

*  Definition:

*  ===========

*

*       SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,

*                          WORK, IWORK, INFO )

*

*       .. Scalar Arguments ..

*       CHARACTER          NORM

*       INTEGER            INFO, KL, KU, LDAB, N

*       DOUBLE PRECISION   ANORM, RCOND

*       ..

*       .. Array Arguments ..

*       INTEGER            IPIV( * ), IWORK( * )

*       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )

*       ..

*

*

*> \par Purpose:

*  =============

*>

*> \verbatim

*>

*> DGBCON estimates the reciprocal of the condition number of a real

*> general band matrix A, in either the 1-norm or the infinity-norm,

*> using the LU factorization computed by DGBTRF.

*>

*> An estimate is obtained for norm(inv(A)), and the reciprocal of the

*> condition number is computed as

*>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).

*> \endverbatim

*

*  Arguments:

*  ==========

*

*> \param[in] NORM

*> \verbatim

*>          NORM is CHARACTER*1

*>          Specifies whether the 1-norm condition number or the

*>          infinity-norm condition number is required:

*>          = '1' or 'O':  1-norm;

*>          = 'I':         Infinity-norm.

*> \endverbatim

*>

*> \param[in] N

*> \verbatim

*>          N is INTEGER

*>          The order of the matrix A.  N >= 0.

*> \endverbatim

*>

*> \param[in] KL

*> \verbatim

*>          KL is INTEGER

*>          The number of subdiagonals within the band of A.  KL >= 0.

*> \endverbatim

*>

*> \param[in] KU

*> \verbatim

*>          KU is INTEGER

*>          The number of superdiagonals within the band of A.  KU >= 0.

*> \endverbatim

*>

*> \param[in] AB

*> \verbatim

*>          AB is DOUBLE PRECISION array, dimension (LDAB,N)

*>          Details of the LU factorization of the band matrix A, as

*>          computed by DGBTRF.  U is stored as an upper triangular band

*>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and

*>          the multipliers used during the factorization are stored in

*>          rows KL+KU+2 to 2*KL+KU+1.

*> \endverbatim

*>

*> \param[in] LDAB

*> \verbatim

*>          LDAB is INTEGER

*>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.

*> \endverbatim

*>

*> \param[in] IPIV

*> \verbatim

*>          IPIV is INTEGER array, dimension (N)

*>          The pivot indices; for 1 <= i <= N, row i of the matrix was

*>          interchanged with row IPIV(i).

*> \endverbatim

*>

*> \param[in] ANORM

*> \verbatim

*>          ANORM is DOUBLE PRECISION

*>          If NORM = '1' or 'O', the 1-norm of the original matrix A.

*>          If NORM = 'I', the infinity-norm of the original matrix A.

*> \endverbatim

*>

*> \param[out] RCOND

*> \verbatim

*>          RCOND is DOUBLE PRECISION

*>          The reciprocal of the condition number of the matrix A,

*>          computed as RCOND = 1/(norm(A) * norm(inv(A))).

*> \endverbatim

*>

*> \param[out] WORK

*> \verbatim

*>          WORK is DOUBLE PRECISION array, dimension (3*N)

*> \endverbatim

*>

*> \param[out] IWORK

*> \verbatim

*>          IWORK is INTEGER array, dimension (N)

*> \endverbatim

*>

*> \param[out] INFO

*> \verbatim

*>          INFO is INTEGER

*>          = 0:  successful exit

*>          < 0: if INFO = -i, the i-th argument had an illegal value

*> \endverbatim

*

*  Authors:

*  ========

*

*> \author Univ. of Tennessee

*> \author Univ. of California Berkeley

*> \author Univ. of Colorado Denver

*> \author NAG Ltd.

*

*> \ingroup doubleGBcomputational

*

*  =====================================================================


      SUBROUTINE dgbcon( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,

     $                   WORK, IWORK, INFO )

*

*  -- LAPACK computational routine --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*

*     .. Scalar Arguments ..

      CHARACTER          NORM

      INTEGER            INFO, KL, KU, LDAB, N

      DOUBLE PRECISION   ANORM, RCOND

*     ..

*     .. Array Arguments ..

      INTEGER            IPIV( * ), IWORK( * )

      DOUBLE PRECISION   AB( LDAB, * ), WORK( * )

*     ..

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ONE, ZERO

      parameter( one = 1.0d+0, zero = 0.0d+0 )

*     ..

*     .. Local Scalars ..

      LOGICAL            LNOTI, ONENRM

      CHARACTER          NORMIN

      INTEGER            IX, J, JP, KASE, KASE1, KD, LM

      DOUBLE PRECISION   AINVNM, SCALE, SMLNUM, T

*     ..

*     .. Local Arrays ..

      INTEGER            ISAVE( 3 )

*     ..

*     .. External Functions ..

      LOGICAL            LSAME

      INTEGER            IDAMAX

      DOUBLE PRECISION   DDOT, DLAMCH

      EXTERNAL           lsame, idamax, ddot, dlamch

*     ..

*     .. External Subroutines ..

      EXTERNAL           daxpy, dlacn2, dlatbs, drscl, xerbla

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC          abs, min

*     ..

*     .. Executable Statements ..

*

*     Test the input parameters.

*

      info = 0

      onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )

      IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN

         info = -1

      ELSE IF( n.LT.0 ) THEN

         info = -2

      ELSE IF( kl.LT.0 ) THEN

         info = -3

      ELSE IF( ku.LT.0 ) THEN

         info = -4

      ELSE IF( ldab.LT.2*kl+ku+1 ) THEN

         info = -6

      ELSE IF( anorm.LT.zero ) THEN

         info = -8

      END IF

      IF( info.NE.0 ) THEN

         CALL xerbla( 'DGBCON', -info )

         RETURN

      END IF

*

*     Quick return if possible

*

      rcond = zero

      IF( n.EQ.0 ) THEN

         rcond = one

         RETURN

      ELSE IF( anorm.EQ.zero ) THEN

         RETURN

      END IF

*

      smlnum = dlamch( 'Safe minimum' )

*

*     Estimate the norm of inv(A).

*

      ainvnm = zero

      normin = 'N'

      IF( onenrm ) THEN

         kase1 = 1

      ELSE

         kase1 = 2

      END IF

      kd = kl + ku + 1

      lnoti = kl.GT.0

      kase = 0

   10 CONTINUE

      CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )

      IF( kase.NE.0 ) THEN

         IF( kase.EQ.kase1 ) THEN

*

*           Multiply by inv(L).

*

            IF( lnoti ) THEN

               DO 20 j = 1, n - 1

                  lm = min( kl, n-j )

                  jp = ipiv( j )

                  t = work( jp )

                  IF( jp.NE.j ) THEN

                     work( jp ) = work( j )

                     work( j ) = t

                  END IF

                  CALL daxpy( lm, -t, ab( kd+1, j ), 1, work( j+1 ), 1 )

   20          CONTINUE

            END IF

*

*           Multiply by inv(U).

*

            CALL dlatbs( 'Upper', 'No transpose', 'Non-unit', normin, n,

     $                   kl+ku, ab, ldab, work, scale, work( 2*n+1 ),

     $                   info )

         ELSE

*

*           Multiply by inv(U**T).

*

            CALL dlatbs( 'Upper', 'Transpose', 'Non-unit', normin, n,

     $                   kl+ku, ab, ldab, work, scale, work( 2*n+1 ),

     $                   info )

*

*           Multiply by inv(L**T).

*

            IF( lnoti ) THEN

               DO 30 j = n - 1, 1, -1

                  lm = min( kl, n-j )

                  work( j ) = work( j ) - ddot( lm, ab( kd+1, j ), 1,

     $                        work( j+1 ), 1 )

                  jp = ipiv( j )

                  IF( jp.NE.j ) THEN

                     t = work( jp )

                     work( jp ) = work( j )

                     work( j ) = t

                  END IF

   30          CONTINUE

            END IF

         END IF

*

*        Divide X by 1/SCALE if doing so will not cause overflow.

*

         normin = 'Y'

         IF( scale.NE.one ) THEN

            ix = idamax( n, work, 1 )

            IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )

     $         GO TO 40

            CALL drscl( n, scale, work, 1 )

         END IF

         GO TO 10

      END IF

*

*     Compute the estimate of the reciprocal condition number.

*

      IF( ainvnm.NE.zero )

     $   rcond = ( one / ainvnm ) / anorm

*

   40 CONTINUE

      RETURN

*

*     End of DGBCON

*


      END

xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60

dgbcon
subroutine dgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
DGBCON
Definition dgbcon.f:146

drscl
subroutine drscl(n, sa, sx, incx)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition drscl.f:84

dlatbs
subroutine dlatbs(uplo, trans, diag, normin, n, kd, ab, ldab, x, scale, cnorm, info)
DLATBS solves a triangular banded system of equations.
Definition dlatbs.f:242

dlacn2
subroutine dlacn2(n, v, x, isgn, est, kase, isave)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition dlacn2.f:136

daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89

min
#define min(a, b)
Definition macros.h:20