1*> \brief \b CGBT01 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 CGBT01( 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* REAL RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> CGBT01 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 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 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*> CGBTRF. 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 CGBTRF 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 CGBTRF. 100*> \endverbatim 101*> 102*> \param[out] WORK 103*> \verbatim 104*> WORK is COMPLEX array, dimension (2*KL+KU+1) 105*> \endverbatim 106*> 107*> \param[out] RESID 108*> \verbatim 109*> RESID is REAL 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*> \ingroup complex_lin 122* 123* ===================================================================== 124 SUBROUTINE CGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 125 $ RESID ) 126* 127* -- LAPACK test routine -- 128* -- LAPACK is a software package provided by Univ. of Tennessee, -- 129* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 130* 131* .. Scalar Arguments .. 132 INTEGER KL, KU, LDA, LDAFAC, M, N 133 REAL RESID 134* .. 135* .. Array Arguments .. 136 INTEGER IPIV( * ) 137 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 138* .. 139* 140* ===================================================================== 141* 142* .. Parameters .. 143 REAL ZERO, ONE 144 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 145* .. 146* .. Local Scalars .. 147 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ 148 REAL ANORM, EPS 149 COMPLEX T 150* .. 151* .. External Functions .. 152 REAL SCASUM, SLAMCH 153 EXTERNAL SCASUM, SLAMCH 154* .. 155* .. External Subroutines .. 156 EXTERNAL CAXPY, CCOPY 157* .. 158* .. Intrinsic Functions .. 159 INTRINSIC CMPLX, MAX, MIN, REAL 160* .. 161* .. Executable Statements .. 162* 163* Quick exit if M = 0 or N = 0. 164* 165 RESID = ZERO 166 IF( M.LE.0 .OR. N.LE.0 ) 167 $ RETURN 168* 169* Determine EPS and the norm of A. 170* 171 EPS = SLAMCH( 'Epsilon' ) 172 KD = KU + 1 173 ANORM = ZERO 174 DO 10 J = 1, N 175 I1 = MAX( KD+1-J, 1 ) 176 I2 = MIN( KD+M-J, KL+KD ) 177 IF( I2.GE.I1 ) 178 $ ANORM = MAX( ANORM, SCASUM( I2-I1+1, A( I1, J ), 1 ) ) 179 10 CONTINUE 180* 181* Compute one column at a time of L*U - A. 182* 183 KD = KL + KU + 1 184 DO 40 J = 1, N 185* 186* Copy the J-th column of U to WORK. 187* 188 JU = MIN( KL+KU, J-1 ) 189 JL = MIN( KL, M-J ) 190 LENJ = MIN( M, J ) - J + JU + 1 191 IF( LENJ.GT.0 ) THEN 192 CALL CCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 193 DO 20 I = LENJ + 1, JU + JL + 1 194 WORK( I ) = ZERO 195 20 CONTINUE 196* 197* Multiply by the unit lower triangular matrix L. Note that L 198* is stored as a product of transformations and permutations. 199* 200 DO 30 I = MIN( M-1, J ), J - JU, -1 201 IL = MIN( KL, M-I ) 202 IF( IL.GT.0 ) THEN 203 IW = I - J + JU + 1 204 T = WORK( IW ) 205 CALL CAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 206 $ 1 ) 207 IP = IPIV( I ) 208 IF( I.NE.IP ) THEN 209 IP = IP - J + JU + 1 210 WORK( IW ) = WORK( IP ) 211 WORK( IP ) = T 212 END IF 213 END IF 214 30 CONTINUE 215* 216* Subtract the corresponding column of A. 217* 218 JUA = MIN( JU, KU ) 219 IF( JUA+JL+1.GT.0 ) 220 $ CALL CAXPY( JUA+JL+1, -CMPLX( ONE ), A( KU+1-JUA, J ), 1, 221 $ WORK( JU+1-JUA ), 1 ) 222* 223* Compute the 1-norm of the column. 224* 225 RESID = MAX( RESID, SCASUM( JU+JL+1, WORK, 1 ) ) 226 END IF 227 40 CONTINUE 228* 229* Compute norm(L*U - A) / ( N * norm(A) * EPS ) 230* 231 IF( ANORM.LE.ZERO ) THEN 232 IF( RESID.NE.ZERO ) 233 $ RESID = ONE / EPS 234 ELSE 235 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 236 END IF 237* 238 RETURN 239* 240* End of CGBT01 241* 242 END 243