SRC/zgbcon.f

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