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*> \ingroup complex16OTHERcomputational 129* 130* ===================================================================== 131 SUBROUTINE ZPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, 132 $ RWORK, INFO ) 133* 134* -- LAPACK computational routine -- 135* -- LAPACK is a software package provided by Univ. of Tennessee, -- 136* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 137* 138* .. Scalar Arguments .. 139 CHARACTER UPLO 140 INTEGER INFO, KD, LDAB, N 141 DOUBLE PRECISION ANORM, RCOND 142* .. 143* .. Array Arguments .. 144 DOUBLE PRECISION RWORK( * ) 145 COMPLEX*16 AB( LDAB, * ), WORK( * ) 146* .. 147* 148* ===================================================================== 149* 150* .. Parameters .. 151 DOUBLE PRECISION ONE, ZERO 152 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 153* .. 154* .. Local Scalars .. 155 LOGICAL UPPER 156 CHARACTER NORMIN 157 INTEGER IX, KASE 158 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 159 COMPLEX*16 ZDUM 160* .. 161* .. Local Arrays .. 162 INTEGER ISAVE( 3 ) 163* .. 164* .. External Functions .. 165 LOGICAL LSAME 166 INTEGER IZAMAX 167 DOUBLE PRECISION DLAMCH 168 EXTERNAL LSAME, IZAMAX, DLAMCH 169* .. 170* .. External Subroutines .. 171 EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATBS 172* .. 173* .. Intrinsic Functions .. 174 INTRINSIC ABS, DBLE, DIMAG 175* .. 176* .. Statement Functions .. 177 DOUBLE PRECISION CABS1 178* .. 179* .. Statement Function definitions .. 180 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 181* .. 182* .. Executable Statements .. 183* 184* Test the input parameters. 185* 186 INFO = 0 187 UPPER = LSAME( UPLO, 'U' ) 188 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 189 INFO = -1 190 ELSE IF( N.LT.0 ) THEN 191 INFO = -2 192 ELSE IF( KD.LT.0 ) THEN 193 INFO = -3 194 ELSE IF( LDAB.LT.KD+1 ) THEN 195 INFO = -5 196 ELSE IF( ANORM.LT.ZERO ) THEN 197 INFO = -6 198 END IF 199 IF( INFO.NE.0 ) THEN 200 CALL XERBLA( 'ZPBCON', -INFO ) 201 RETURN 202 END IF 203* 204* Quick return if possible 205* 206 RCOND = ZERO 207 IF( N.EQ.0 ) THEN 208 RCOND = ONE 209 RETURN 210 ELSE IF( ANORM.EQ.ZERO ) THEN 211 RETURN 212 END IF 213* 214 SMLNUM = DLAMCH( 'Safe minimum' ) 215* 216* Estimate the 1-norm of the inverse. 217* 218 KASE = 0 219 NORMIN = 'N' 220 10 CONTINUE 221 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 222 IF( KASE.NE.0 ) THEN 223 IF( UPPER ) THEN 224* 225* Multiply by inv(U**H). 226* 227 CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit', 228 $ NORMIN, N, KD, AB, LDAB, WORK, SCALEL, RWORK, 229 $ INFO ) 230 NORMIN = 'Y' 231* 232* Multiply by inv(U). 233* 234 CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 235 $ KD, AB, LDAB, WORK, SCALEU, RWORK, INFO ) 236 ELSE 237* 238* Multiply by inv(L). 239* 240 CALL ZLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 241 $ KD, AB, LDAB, WORK, SCALEL, RWORK, INFO ) 242 NORMIN = 'Y' 243* 244* Multiply by inv(L**H). 245* 246 CALL ZLATBS( 'Lower', 'Conjugate transpose', 'Non-unit', 247 $ NORMIN, N, KD, AB, LDAB, WORK, SCALEU, RWORK, 248 $ INFO ) 249 END IF 250* 251* Multiply by 1/SCALE if doing so will not cause overflow. 252* 253 SCALE = SCALEL*SCALEU 254 IF( SCALE.NE.ONE ) THEN 255 IX = IZAMAX( N, WORK, 1 ) 256 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 257 $ GO TO 20 258 CALL ZDRSCL( N, SCALE, WORK, 1 ) 259 END IF 260 GO TO 10 261 END IF 262* 263* Compute the estimate of the reciprocal condition number. 264* 265 IF( AINVNM.NE.ZERO ) 266 $ RCOND = ( ONE / AINVNM ) / ANORM 267* 268 20 CONTINUE 269* 270 RETURN 271* 272* End of ZPBCON 273* 274 END 275