1*> \brief \b ZGBT01 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 12* RESID ) 13* 14* .. Scalar Arguments .. 15* INTEGER KL, KU, LDA, LDAFAC, M, N 16* DOUBLE PRECISION RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZGBT01 reconstructs a band matrix A from its L*U factorization and 30*> computes the residual: 31*> norm(L*U - A) / ( N * norm(A) * EPS ), 32*> where EPS is the machine epsilon. 33*> 34*> The expression L*U - A is computed one column at a time, so A and 35*> AFAC are not modified. 36*> \endverbatim 37* 38* Arguments: 39* ========== 40* 41*> \param[in] M 42*> \verbatim 43*> M is INTEGER 44*> The number of rows of the matrix A. M >= 0. 45*> \endverbatim 46*> 47*> \param[in] N 48*> \verbatim 49*> N is INTEGER 50*> The number of columns of the matrix A. N >= 0. 51*> \endverbatim 52*> 53*> \param[in] KL 54*> \verbatim 55*> KL is INTEGER 56*> The number of subdiagonals within the band of A. KL >= 0. 57*> \endverbatim 58*> 59*> \param[in] KU 60*> \verbatim 61*> KU is INTEGER 62*> The number of superdiagonals within the band of A. KU >= 0. 63*> \endverbatim 64*> 65*> \param[in,out] A 66*> \verbatim 67*> A is COMPLEX*16 array, dimension (LDA,N) 68*> The original matrix A in band storage, stored in rows 1 to 69*> KL+KU+1. 70*> \endverbatim 71*> 72*> \param[in] LDA 73*> \verbatim 74*> LDA is INTEGER. 75*> The leading dimension of the array A. LDA >= max(1,KL+KU+1). 76*> \endverbatim 77*> 78*> \param[in] AFAC 79*> \verbatim 80*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) 81*> The factored form of the matrix A. AFAC contains the banded 82*> factors L and U from the L*U factorization, as computed by 83*> ZGBTRF. U is stored as an upper triangular band matrix with 84*> KL+KU superdiagonals in rows 1 to KL+KU+1, and the 85*> multipliers used during the factorization are stored in rows 86*> KL+KU+2 to 2*KL+KU+1. See ZGBTRF for further details. 87*> \endverbatim 88*> 89*> \param[in] LDAFAC 90*> \verbatim 91*> LDAFAC is INTEGER 92*> The leading dimension of the array AFAC. 93*> LDAFAC >= max(1,2*KL*KU+1). 94*> \endverbatim 95*> 96*> \param[in] IPIV 97*> \verbatim 98*> IPIV is INTEGER array, dimension (min(M,N)) 99*> The pivot indices from ZGBTRF. 100*> \endverbatim 101*> 102*> \param[out] WORK 103*> \verbatim 104*> WORK is COMPLEX*16 array, dimension (2*KL+KU+1) 105*> \endverbatim 106*> 107*> \param[out] RESID 108*> \verbatim 109*> RESID is DOUBLE PRECISION 110*> norm(L*U - A) / ( N * norm(A) * EPS ) 111*> \endverbatim 112* 113* Authors: 114* ======== 115* 116*> \author Univ. of Tennessee 117*> \author Univ. of California Berkeley 118*> \author Univ. of Colorado Denver 119*> \author NAG Ltd. 120* 121*> \date November 2011 122* 123*> \ingroup complex16_lin 124* 125* ===================================================================== 126 SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 127 $ RESID ) 128* 129* -- LAPACK test routine (version 3.4.0) -- 130* -- LAPACK is a software package provided by Univ. of Tennessee, -- 131* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 132* November 2011 133* 134* .. Scalar Arguments .. 135 INTEGER KL, KU, LDA, LDAFAC, M, N 136 DOUBLE PRECISION RESID 137* .. 138* .. Array Arguments .. 139 INTEGER IPIV( * ) 140 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 141* .. 142* 143* ===================================================================== 144* 145* .. Parameters .. 146 DOUBLE PRECISION ZERO, ONE 147 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 148* .. 149* .. Local Scalars .. 150 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ 151 DOUBLE PRECISION ANORM, EPS 152 COMPLEX*16 T 153* .. 154* .. External Functions .. 155 DOUBLE PRECISION DLAMCH, DZASUM 156 EXTERNAL DLAMCH, DZASUM 157* .. 158* .. External Subroutines .. 159 EXTERNAL ZAXPY, ZCOPY 160* .. 161* .. Intrinsic Functions .. 162 INTRINSIC DBLE, DCMPLX, MAX, MIN 163* .. 164* .. Executable Statements .. 165* 166* Quick exit if M = 0 or N = 0. 167* 168 RESID = ZERO 169 IF( M.LE.0 .OR. N.LE.0 ) 170 $ RETURN 171* 172* Determine EPS and the norm of A. 173* 174 EPS = DLAMCH( 'Epsilon' ) 175 KD = KU + 1 176 ANORM = ZERO 177 DO 10 J = 1, N 178 I1 = MAX( KD+1-J, 1 ) 179 I2 = MIN( KD+M-J, KL+KD ) 180 IF( I2.GE.I1 ) 181 $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) ) 182 10 CONTINUE 183* 184* Compute one column at a time of L*U - A. 185* 186 KD = KL + KU + 1 187 DO 40 J = 1, N 188* 189* Copy the J-th column of U to WORK. 190* 191 JU = MIN( KL+KU, J-1 ) 192 JL = MIN( KL, M-J ) 193 LENJ = MIN( M, J ) - J + JU + 1 194 IF( LENJ.GT.0 ) THEN 195 CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 196 DO 20 I = LENJ + 1, JU + JL + 1 197 WORK( I ) = ZERO 198 20 CONTINUE 199* 200* Multiply by the unit lower triangular matrix L. Note that L 201* is stored as a product of transformations and permutations. 202* 203 DO 30 I = MIN( M-1, J ), J - JU, -1 204 IL = MIN( KL, M-I ) 205 IF( IL.GT.0 ) THEN 206 IW = I - J + JU + 1 207 T = WORK( IW ) 208 CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 209 $ 1 ) 210 IP = IPIV( I ) 211 IF( I.NE.IP ) THEN 212 IP = IP - J + JU + 1 213 WORK( IW ) = WORK( IP ) 214 WORK( IP ) = T 215 END IF 216 END IF 217 30 CONTINUE 218* 219* Subtract the corresponding column of A. 220* 221 JUA = MIN( JU, KU ) 222 IF( JUA+JL+1.GT.0 ) 223 $ CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ), 224 $ 1, WORK( JU+1-JUA ), 1 ) 225* 226* Compute the 1-norm of the column. 227* 228 RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) ) 229 END IF 230 40 CONTINUE 231* 232* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 233* 234 IF( ANORM.LE.ZERO ) THEN 235 IF( RESID.NE.ZERO ) 236 $ RESID = ONE / EPS 237 ELSE 238 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 239 END IF 240* 241 RETURN 242* 243* End of ZGBT01 244* 245 END 246