SRC/zpbcon.f

   1 *> \brief \b ZPBCON
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZPBCON + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpbcon.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpbcon.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpbcon.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
  22 *                          RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, KD, LDAB, N
  27 *       DOUBLE PRECISION   ANORM, RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   RWORK( * )
  31 *       COMPLEX*16         AB( LDAB, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZPBCON estimates the reciprocal of the condition number (in the
  41 *> 1-norm) of a complex Hermitian positive definite band matrix using
  42 *> the Cholesky factorization A = U**H*U or A = L*L**H computed by
  43 *> ZPBTRF.
  44 *>
  45 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
  46 *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] UPLO
  53 *> \verbatim
  54 *>          UPLO is CHARACTER*1
  55 *>          = 'U':  Upper triangular factor stored in AB;
  56 *>          = 'L':  Lower triangular factor stored in AB.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] N
  60 *> \verbatim
  61 *>          N is INTEGER
  62 *>          The order of the matrix A.  N >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] KD
  66 *> \verbatim
  67 *>          KD is INTEGER
  68 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  69 *>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] AB
  73 *> \verbatim
  74 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  75 *>          The triangular factor U or L from the Cholesky factorization
  76 *>          A = U**H*U or A = L*L**H of the band matrix A, stored in the
  77 *>          first KD+1 rows of the array.  The j-th column of U or L is
  78 *>          stored in the j-th column of the array AB as follows:
  79 *>          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
  80 *>          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).
  81 *> \endverbatim
  82 *>
  83 *> \param[in] LDAB
  84 *> \verbatim
  85 *>          LDAB is INTEGER
  86 *>          The leading dimension of the array AB.  LDAB >= KD+1.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] ANORM
  90 *> \verbatim
  91 *>          ANORM is DOUBLE PRECISION
  92 *>          The 1-norm (or infinity-norm) of the Hermitian band matrix A.
  93 *> \endverbatim
  94 *>
  95 *> \param[out] RCOND
  96 *> \verbatim
  97 *>          RCOND is DOUBLE PRECISION
  98 *>          The reciprocal of the condition number of the matrix A,
  99 *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
 100 *>          estimate of the 1-norm of inv(A) computed in this routine.
 101 *> \endverbatim
 102 *>
 103 *> \param[out] WORK
 104 *> \verbatim
 105 *>          WORK is COMPLEX*16 array, dimension (2*N)
 106 *> \endverbatim
 107 *>
 108 *> \param[out] RWORK
 109 *> \verbatim
 110 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 111 *> \endverbatim
 112 *>
 113 *> \param[out] INFO
 114 *> \verbatim
 115 *>          INFO is INTEGER
 116 *>          = 0:  successful exit
 117 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 118 *> \endverbatim
 119 *
 120 *  Authors:
 121 *  ========
 122 *
 123 *> \author Univ. of Tennessee
 124 *> \author Univ. of California Berkeley
 125 *> \author Univ. of Colorado Denver
 126 *> \author NAG Ltd.
 127 *
 128 *> \date November 2011
 129 *
 130 *> \ingroup complex16OTHERcomputational
 131 *
 132 *  =====================================================================
 133       SUBROUTINE ZPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
 134      $                   RWORK, INFO )
 135 *
 136 *  -- LAPACK computational routine (version 3.4.0) --
 137 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 138 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 139 *     November 2011
 140 *
 141 *     .. Scalar Arguments ..
 142       CHARACTER          UPLO
 143       INTEGER            INFO, KD, LDAB, N
 144       DOUBLE PRECISION   ANORM, RCOND
 145 *     ..
 146 *     .. Array Arguments ..
 147       DOUBLE PRECISION   RWORK( * )
 148       COMPLEX*16         AB( LDAB, * ), WORK( * )
 149 *     ..
 150 *
 151 *  =====================================================================
 152 *
 153 *     .. Parameters ..
 154       DOUBLE PRECISION   ONE, ZERO
 155       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 156 *     ..
 157 *     .. Local Scalars ..
 158       LOGICAL            UPPER
 159       CHARACTER          NORMIN
 160       INTEGER            IX, KASE
 161       DOUBLE PRECISION   AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
 162       COMPLEX*16         ZDUM
 163 *     ..
 164 *     .. Local Arrays ..
 165       INTEGER            ISAVE( 3 )
 166 *     ..
 167 *     .. External Functions ..
 168       LOGICAL            LSAME
 169       INTEGER            IZAMAX
 170       DOUBLE PRECISION   DLAMCH
 171       EXTERNAL           LSAME, IZAMAX, DLAMCH
 172 *     ..
 173 *     .. External Subroutines ..
 174       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLATBS
 175 *     ..
 176 *     .. Intrinsic Functions ..
 177       INTRINSIC          ABS, DBLE, DIMAG
 178 *     ..
 179 *     .. Statement Functions ..
 180       DOUBLE PRECISION   CABS1
 181 *     ..
 182 *     .. Statement Function definitions ..
 183       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187 *     Test the input parameters.
 188 *
 189       INFO = 0
 190       UPPER = LSAME( UPLO, 'U' )
 191       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 192          INFO = -1
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -2
 195       ELSE IF( KD.LT.0 ) THEN
 196          INFO = -3
 197       ELSE IF( LDAB.LT.KD+1 ) THEN
 198          INFO = -5
 199       ELSE IF( ANORM.LT.ZERO ) THEN
 200          INFO = -6
 201       END IF
 202       IF( INFO.NE.0 ) THEN
 203          CALL XERBLA( 'ZPBCON', -INFO )
 204          RETURN
 205       END IF
 206 *
 207 *     Quick return if possible
 208 *
 209       RCOND = ZERO
 210       IF( N.EQ.0 ) THEN
 211          RCOND = ONE
 212          RETURN
 213       ELSE IF( ANORM.EQ.ZERO ) THEN
 214          RETURN
 215       END IF
 216 *
 217       SMLNUM = DLAMCH( 'Safe minimum' )
 218 *
 219 *     Estimate the 1-norm of the inverse.
 220 *
 221       KASE = 0
 222       NORMIN = 'N'
 223    10 CONTINUE
 224       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
 225       IF( KASE.NE.0 ) THEN
 226          IF( UPPER ) THEN
 227 *
 228 *           Multiply by inv(U**H).
 229 *
 230             CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
 231      $                   NORMIN, N, KD, AB, LDAB, WORK, SCALEL, RWORK,
 232      $                   INFO )
 233             NORMIN = 'Y'
 234 *
 235 *           Multiply by inv(U).
 236 *
 237             CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
 238      $                   KD, AB, LDAB, WORK, SCALEU, RWORK, INFO )
 239          ELSE
 240 *
 241 *           Multiply by inv(L).
 242 *
 243             CALL ZLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
 244      $                   KD, AB, LDAB, WORK, SCALEL, RWORK, INFO )
 245             NORMIN = 'Y'
 246 *
 247 *           Multiply by inv(L**H).
 248 *
 249             CALL ZLATBS( 'Lower', 'Conjugate transpose', 'Non-unit',
 250      $                   NORMIN, N, KD, AB, LDAB, WORK, SCALEU, RWORK,
 251      $                   INFO )
 252          END IF
 253 *
 254 *        Multiply by 1/SCALE if doing so will not cause overflow.
 255 *
 256          SCALE = SCALEL*SCALEU
 257          IF( SCALE.NE.ONE ) THEN
 258             IX = IZAMAX( N, WORK, 1 )
 259             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
 260      $         GO TO 20
 261             CALL ZDRSCL( N, SCALE, WORK, 1 )
 262          END IF
 263          GO TO 10
 264       END IF
 265 *
 266 *     Compute the estimate of the reciprocal condition number.
 267 *
 268       IF( AINVNM.NE.ZERO )
 269      $   RCOND = ( ONE / AINVNM ) / ANORM
 270 *
 271    20 CONTINUE
 272 *
 273       RETURN
 274 *
 275 *     End of ZPBCON
 276 *
 277       END