1*> \brief \b ZLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded matrices. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLA_GBRCOND_C + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_gbrcond_c.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_gbrcond_c.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_gbrcond_c.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, 22* LDAB, AFB, LDAFB, IPIV, 23* C, CAPPLY, INFO, WORK, 24* RWORK ) 25* 26* .. Scalar Arguments .. 27* CHARACTER TRANS 28* LOGICAL CAPPLY 29* INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO 30* .. 31* .. Array Arguments .. 32* INTEGER IPIV( * ) 33* COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ) 34* DOUBLE PRECISION C( * ), RWORK( * ) 35* 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> ZLA_GBRCOND_C Computes the infinity norm condition number of 44*> op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector. 45*> \endverbatim 46* 47* Arguments: 48* ========== 49* 50*> \param[in] TRANS 51*> \verbatim 52*> TRANS is CHARACTER*1 53*> Specifies the form of the system of equations: 54*> = 'N': A * X = B (No transpose) 55*> = 'T': A**T * X = B (Transpose) 56*> = 'C': A**H * X = B (Conjugate Transpose = Transpose) 57*> \endverbatim 58*> 59*> \param[in] N 60*> \verbatim 61*> N is INTEGER 62*> The number of linear equations, i.e., the order of the 63*> matrix A. N >= 0. 64*> \endverbatim 65*> 66*> \param[in] KL 67*> \verbatim 68*> KL is INTEGER 69*> The number of subdiagonals within the band of A. KL >= 0. 70*> \endverbatim 71*> 72*> \param[in] KU 73*> \verbatim 74*> KU is INTEGER 75*> The number of superdiagonals within the band of A. KU >= 0. 76*> \endverbatim 77*> 78*> \param[in] AB 79*> \verbatim 80*> AB is COMPLEX*16 array, dimension (LDAB,N) 81*> On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 82*> The j-th column of A is stored in the j-th column of the 83*> array AB as follows: 84*> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 85*> \endverbatim 86*> 87*> \param[in] LDAB 88*> \verbatim 89*> LDAB is INTEGER 90*> The leading dimension of the array AB. LDAB >= KL+KU+1. 91*> \endverbatim 92*> 93*> \param[in] AFB 94*> \verbatim 95*> AFB is COMPLEX*16 array, dimension (LDAFB,N) 96*> Details of the LU factorization of the band matrix A, as 97*> computed by ZGBTRF. U is stored as an upper triangular 98*> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 99*> and the multipliers used during the factorization are stored 100*> in rows KL+KU+2 to 2*KL+KU+1. 101*> \endverbatim 102*> 103*> \param[in] LDAFB 104*> \verbatim 105*> LDAFB is INTEGER 106*> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 107*> \endverbatim 108*> 109*> \param[in] IPIV 110*> \verbatim 111*> IPIV is INTEGER array, dimension (N) 112*> The pivot indices from the factorization A = P*L*U 113*> as computed by ZGBTRF; row i of the matrix was interchanged 114*> with row IPIV(i). 115*> \endverbatim 116*> 117*> \param[in] C 118*> \verbatim 119*> C is DOUBLE PRECISION array, dimension (N) 120*> The vector C in the formula op(A) * inv(diag(C)). 121*> \endverbatim 122*> 123*> \param[in] CAPPLY 124*> \verbatim 125*> CAPPLY is LOGICAL 126*> If .TRUE. then access the vector C in the formula above. 127*> \endverbatim 128*> 129*> \param[out] INFO 130*> \verbatim 131*> INFO is INTEGER 132*> = 0: Successful exit. 133*> i > 0: The ith argument is invalid. 134*> \endverbatim 135*> 136*> \param[out] WORK 137*> \verbatim 138*> WORK is COMPLEX*16 array, dimension (2*N). 139*> Workspace. 140*> \endverbatim 141*> 142*> \param[out] RWORK 143*> \verbatim 144*> RWORK is DOUBLE PRECISION array, dimension (N). 145*> Workspace. 146*> \endverbatim 147* 148* Authors: 149* ======== 150* 151*> \author Univ. of Tennessee 152*> \author Univ. of California Berkeley 153*> \author Univ. of Colorado Denver 154*> \author NAG Ltd. 155* 156*> \ingroup complex16GBcomputational 157* 158* ===================================================================== 159 DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, 160 $ LDAB, AFB, LDAFB, IPIV, 161 $ C, CAPPLY, INFO, WORK, 162 $ RWORK ) 163* 164* -- LAPACK computational routine -- 165* -- LAPACK is a software package provided by Univ. of Tennessee, -- 166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 167* 168* .. Scalar Arguments .. 169 CHARACTER TRANS 170 LOGICAL CAPPLY 171 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO 172* .. 173* .. Array Arguments .. 174 INTEGER IPIV( * ) 175 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ) 176 DOUBLE PRECISION C( * ), RWORK( * ) 177* 178* 179* ===================================================================== 180* 181* .. Local Scalars .. 182 LOGICAL NOTRANS 183 INTEGER KASE, I, J 184 DOUBLE PRECISION AINVNM, ANORM, TMP 185 COMPLEX*16 ZDUM 186* .. 187* .. Local Arrays .. 188 INTEGER ISAVE( 3 ) 189* .. 190* .. External Functions .. 191 LOGICAL LSAME 192 EXTERNAL LSAME 193* .. 194* .. External Subroutines .. 195 EXTERNAL ZLACN2, ZGBTRS, XERBLA 196* .. 197* .. Intrinsic Functions .. 198 INTRINSIC ABS, MAX 199* .. 200* .. Statement Functions .. 201 DOUBLE PRECISION CABS1 202* .. 203* .. Statement Function Definitions .. 204 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 205* .. 206* .. Executable Statements .. 207 ZLA_GBRCOND_C = 0.0D+0 208* 209 INFO = 0 210 NOTRANS = LSAME( TRANS, 'N' ) 211 IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT. 212 $ LSAME( TRANS, 'C' ) ) THEN 213 INFO = -1 214 ELSE IF( N.LT.0 ) THEN 215 INFO = -2 216 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 217 INFO = -3 218 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 219 INFO = -4 220 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 221 INFO = -6 222 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 223 INFO = -8 224 END IF 225 IF( INFO.NE.0 ) THEN 226 CALL XERBLA( 'ZLA_GBRCOND_C', -INFO ) 227 RETURN 228 END IF 229* 230* Compute norm of op(A)*op2(C). 231* 232 ANORM = 0.0D+0 233 KD = KU + 1 234 KE = KL + 1 235 IF ( NOTRANS ) THEN 236 DO I = 1, N 237 TMP = 0.0D+0 238 IF ( CAPPLY ) THEN 239 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 240 TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J ) 241 END DO 242 ELSE 243 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 244 TMP = TMP + CABS1( AB( KD+I-J, J ) ) 245 END DO 246 END IF 247 RWORK( I ) = TMP 248 ANORM = MAX( ANORM, TMP ) 249 END DO 250 ELSE 251 DO I = 1, N 252 TMP = 0.0D+0 253 IF ( CAPPLY ) THEN 254 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 255 TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J ) 256 END DO 257 ELSE 258 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 259 TMP = TMP + CABS1( AB( KE-I+J, I ) ) 260 END DO 261 END IF 262 RWORK( I ) = TMP 263 ANORM = MAX( ANORM, TMP ) 264 END DO 265 END IF 266* 267* Quick return if possible. 268* 269 IF( N.EQ.0 ) THEN 270 ZLA_GBRCOND_C = 1.0D+0 271 RETURN 272 ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 273 RETURN 274 END IF 275* 276* Estimate the norm of inv(op(A)). 277* 278 AINVNM = 0.0D+0 279* 280 KASE = 0 281 10 CONTINUE 282 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 283 IF( KASE.NE.0 ) THEN 284 IF( KASE.EQ.2 ) THEN 285* 286* Multiply by R. 287* 288 DO I = 1, N 289 WORK( I ) = WORK( I ) * RWORK( I ) 290 END DO 291* 292 IF ( NOTRANS ) THEN 293 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 294 $ IPIV, WORK, N, INFO ) 295 ELSE 296 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 297 $ LDAFB, IPIV, WORK, N, INFO ) 298 ENDIF 299* 300* Multiply by inv(C). 301* 302 IF ( CAPPLY ) THEN 303 DO I = 1, N 304 WORK( I ) = WORK( I ) * C( I ) 305 END DO 306 END IF 307 ELSE 308* 309* Multiply by inv(C**H). 310* 311 IF ( CAPPLY ) THEN 312 DO I = 1, N 313 WORK( I ) = WORK( I ) * C( I ) 314 END DO 315 END IF 316* 317 IF ( NOTRANS ) THEN 318 CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 319 $ LDAFB, IPIV, WORK, N, INFO ) 320 ELSE 321 CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 322 $ IPIV, WORK, N, INFO ) 323 END IF 324* 325* Multiply by R. 326* 327 DO I = 1, N 328 WORK( I ) = WORK( I ) * RWORK( I ) 329 END DO 330 END IF 331 GO TO 10 332 END IF 333* 334* Compute the estimate of the reciprocal condition number. 335* 336 IF( AINVNM .NE. 0.0D+0 ) 337 $ ZLA_GBRCOND_C = 1.0D+0 / AINVNM 338* 339 RETURN 340* 341* End of ZLA_GBRCOND_C 342* 343 END 344