SRC/dgbcon.f

   1 *> \brief \b DGBCON
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGBCON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbcon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbcon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbcon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
  22 *                          WORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          NORM
  26 *       INTEGER            INFO, KL, KU, LDAB, N
  27 *       DOUBLE PRECISION   ANORM, RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IPIV( * ), IWORK( * )
  31 *       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DGBCON estimates the reciprocal of the condition number of a real
  41 *> general band matrix A, in either the 1-norm or the infinity-norm,
  42 *> using the LU factorization computed by DGBTRF.
  43 *>
  44 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
  45 *> condition number is computed as
  46 *>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] NORM
  53 *> \verbatim
  54 *>          NORM is CHARACTER*1
  55 *>          Specifies whether the 1-norm condition number or the
  56 *>          infinity-norm condition number is required:
  57 *>          = '1' or 'O':  1-norm;
  58 *>          = 'I':         Infinity-norm.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>          The order of the matrix A.  N >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] KL
  68 *> \verbatim
  69 *>          KL is INTEGER
  70 *>          The number of subdiagonals within the band of A.  KL >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] KU
  74 *> \verbatim
  75 *>          KU is INTEGER
  76 *>          The number of superdiagonals within the band of A.  KU >= 0.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] AB
  80 *> \verbatim
  81 *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
  82 *>          Details of the LU factorization of the band matrix A, as
  83 *>          computed by DGBTRF.  U is stored as an upper triangular band
  84 *>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
  85 *>          the multipliers used during the factorization are stored in
  86 *>          rows KL+KU+2 to 2*KL+KU+1.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDAB
  90 *> \verbatim
  91 *>          LDAB is INTEGER
  92 *>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] IPIV
  96 *> \verbatim
  97 *>          IPIV is INTEGER array, dimension (N)
  98 *>          The pivot indices; for 1 <= i <= N, row i of the matrix was
  99 *>          interchanged with row IPIV(i).
 100 *> \endverbatim
 101 *>
 102 *> \param[in] ANORM
 103 *> \verbatim
 104 *>          ANORM is DOUBLE PRECISION
 105 *>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
 106 *>          If NORM = 'I', the infinity-norm of the original matrix A.
 107 *> \endverbatim
 108 *>
 109 *> \param[out] RCOND
 110 *> \verbatim
 111 *>          RCOND is DOUBLE PRECISION
 112 *>          The reciprocal of the condition number of the matrix A,
 113 *>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
 114 *> \endverbatim
 115 *>
 116 *> \param[out] WORK
 117 *> \verbatim
 118 *>          WORK is DOUBLE PRECISION array, dimension (3*N)
 119 *> \endverbatim
 120 *>
 121 *> \param[out] IWORK
 122 *> \verbatim
 123 *>          IWORK is INTEGER array, dimension (N)
 124 *> \endverbatim
 125 *>
 126 *> \param[out] INFO
 127 *> \verbatim
 128 *>          INFO is INTEGER
 129 *>          = 0:  successful exit
 130 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 131 *> \endverbatim
 132 *
 133 *  Authors:
 134 *  ========
 135 *
 136 *> \author Univ. of Tennessee
 137 *> \author Univ. of California Berkeley
 138 *> \author Univ. of Colorado Denver
 139 *> \author NAG Ltd.
 140 *
 141 *> \date November 2011
 142 *
 143 *> \ingroup doubleGBcomputational
 144 *
 145 *  =====================================================================
 146       SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
 147      $                   WORK, IWORK, INFO )
 148 *
 149 *  -- LAPACK computational routine (version 3.4.0) --
 150 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 151 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 152 *     November 2011
 153 *
 154 *     .. Scalar Arguments ..
 155       CHARACTER          NORM
 156       INTEGER            INFO, KL, KU, LDAB, N
 157       DOUBLE PRECISION   ANORM, RCOND
 158 *     ..
 159 *     .. Array Arguments ..
 160       INTEGER            IPIV( * ), IWORK( * )
 161       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
 162 *     ..
 163 *
 164 *  =====================================================================
 165 *
 166 *     .. Parameters ..
 167       DOUBLE PRECISION   ONE, ZERO
 168       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 169 *     ..
 170 *     .. Local Scalars ..
 171       LOGICAL            LNOTI, ONENRM
 172       CHARACTER          NORMIN
 173       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
 174       DOUBLE PRECISION   AINVNM, SCALE, SMLNUM, T
 175 *     ..
 176 *     .. Local Arrays ..
 177       INTEGER            ISAVE( 3 )
 178 *     ..
 179 *     .. External Functions ..
 180       LOGICAL            LSAME
 181       INTEGER            IDAMAX
 182       DOUBLE PRECISION   DDOT, DLAMCH
 183       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
 184 *     ..
 185 *     .. External Subroutines ..
 186       EXTERNAL           DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
 187 *     ..
 188 *     .. Intrinsic Functions ..
 189       INTRINSIC          ABS, MIN
 190 *     ..
 191 *     .. Executable Statements ..
 192 *
 193 *     Test the input parameters.
 194 *
 195       INFO = 0
 196       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
 197       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
 198          INFO = -1
 199       ELSE IF( N.LT.0 ) THEN
 200          INFO = -2
 201       ELSE IF( KL.LT.0 ) THEN
 202          INFO = -3
 203       ELSE IF( KU.LT.0 ) THEN
 204          INFO = -4
 205       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
 206          INFO = -6
 207       ELSE IF( ANORM.LT.ZERO ) THEN
 208          INFO = -8
 209       END IF
 210       IF( INFO.NE.0 ) THEN
 211          CALL XERBLA( 'DGBCON', -INFO )
 212          RETURN
 213       END IF
 214 *
 215 *     Quick return if possible
 216 *
 217       RCOND = ZERO
 218       IF( N.EQ.0 ) THEN
 219          RCOND = ONE
 220          RETURN
 221       ELSE IF( ANORM.EQ.ZERO ) THEN
 222          RETURN
 223       END IF
 224 *
 225       SMLNUM = DLAMCH( 'Safe minimum' )
 226 *
 227 *     Estimate the norm of inv(A).
 228 *
 229       AINVNM = ZERO
 230       NORMIN = 'N'
 231       IF( ONENRM ) THEN
 232          KASE1 = 1
 233       ELSE
 234          KASE1 = 2
 235       END IF
 236       KD = KL + KU + 1
 237       LNOTI = KL.GT.0
 238       KASE = 0
 239    10 CONTINUE
 240       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
 241       IF( KASE.NE.0 ) THEN
 242          IF( KASE.EQ.KASE1 ) THEN
 243 *
 244 *           Multiply by inv(L).
 245 *
 246             IF( LNOTI ) THEN
 247                DO 20 J = 1, N - 1
 248                   LM = MIN( KL, N-J )
 249                   JP = IPIV( J )
 250                   T = WORK( JP )
 251                   IF( JP.NE.J ) THEN
 252                      WORK( JP ) = WORK( J )
 253                      WORK( J ) = T
 254                   END IF
 255                   CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
 256    20          CONTINUE
 257             END IF
 258 *
 259 *           Multiply by inv(U).
 260 *
 261             CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
 262      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
 263      $                   INFO )
 264          ELSE
 265 *
 266 *           Multiply by inv(U**T).
 267 *
 268             CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
 269      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
 270      $                   INFO )
 271 *
 272 *           Multiply by inv(L**T).
 273 *
 274             IF( LNOTI ) THEN
 275                DO 30 J = N - 1, 1, -1
 276                   LM = MIN( KL, N-J )
 277                   WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
 278      $                        WORK( J+1 ), 1 )
 279                   JP = IPIV( J )
 280                   IF( JP.NE.J ) THEN
 281                      T = WORK( JP )
 282                      WORK( JP ) = WORK( J )
 283                      WORK( J ) = T
 284                   END IF
 285    30          CONTINUE
 286             END IF
 287          END IF
 288 *
 289 *        Divide X by 1/SCALE if doing so will not cause overflow.
 290 *
 291          NORMIN = 'Y'
 292          IF( SCALE.NE.ONE ) THEN
 293             IX = IDAMAX( N, WORK, 1 )
 294             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
 295      $         GO TO 40
 296             CALL DRSCL( N, SCALE, WORK, 1 )
 297          END IF
 298          GO TO 10
 299       END IF
 300 *
 301 *     Compute the estimate of the reciprocal condition number.
 302 *
 303       IF( AINVNM.NE.ZERO )
 304      $   RCOND = ( ONE / AINVNM ) / ANORM
 305 *
 306    40 CONTINUE
 307       RETURN
 308 *
 309 *     End of DGBCON
 310 *
 311       END