1*> \brief \b SGBT01 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 SGBT01( 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* REAL A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> SGBT01 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 REAL 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 REAL 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*> SGBTRF. 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 SGBTRF 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 SGBTRF. 100*> \endverbatim 101*> 102*> \param[out] WORK 103*> \verbatim 104*> WORK is REAL 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 single_lin 122* 123* ===================================================================== 124 SUBROUTINE SGBT01( 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 REAL 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, T 149* .. 150* .. External Functions .. 151 REAL SASUM, SLAMCH 152 EXTERNAL SASUM, SLAMCH 153* .. 154* .. External Subroutines .. 155 EXTERNAL SAXPY, SCOPY 156* .. 157* .. Intrinsic Functions .. 158 INTRINSIC MAX, MIN, REAL 159* .. 160* .. Executable Statements .. 161* 162* Quick exit if M = 0 or N = 0. 163* 164 RESID = ZERO 165 IF( M.LE.0 .OR. N.LE.0 ) 166 $ RETURN 167* 168* Determine EPS and the norm of A. 169* 170 EPS = SLAMCH( 'Epsilon' ) 171 KD = KU + 1 172 ANORM = ZERO 173 DO 10 J = 1, N 174 I1 = MAX( KD+1-J, 1 ) 175 I2 = MIN( KD+M-J, KL+KD ) 176 IF( I2.GE.I1 ) 177 $ ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) ) 178 10 CONTINUE 179* 180* Compute one column at a time of L*U - A. 181* 182 KD = KL + KU + 1 183 DO 40 J = 1, N 184* 185* Copy the J-th column of U to WORK. 186* 187 JU = MIN( KL+KU, J-1 ) 188 JL = MIN( KL, M-J ) 189 LENJ = MIN( M, J ) - J + JU + 1 190 IF( LENJ.GT.0 ) THEN 191 CALL SCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 192 DO 20 I = LENJ + 1, JU + JL + 1 193 WORK( I ) = ZERO 194 20 CONTINUE 195* 196* Multiply by the unit lower triangular matrix L. Note that L 197* is stored as a product of transformations and permutations. 198* 199 DO 30 I = MIN( M-1, J ), J - JU, -1 200 IL = MIN( KL, M-I ) 201 IF( IL.GT.0 ) THEN 202 IW = I - J + JU + 1 203 T = WORK( IW ) 204 CALL SAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 205 $ 1 ) 206 IP = IPIV( I ) 207 IF( I.NE.IP ) THEN 208 IP = IP - J + JU + 1 209 WORK( IW ) = WORK( IP ) 210 WORK( IP ) = T 211 END IF 212 END IF 213 30 CONTINUE 214* 215* Subtract the corresponding column of A. 216* 217 JUA = MIN( JU, KU ) 218 IF( JUA+JL+1.GT.0 ) 219 $ CALL SAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1, 220 $ WORK( JU+1-JUA ), 1 ) 221* 222* Compute the 1-norm of the column. 223* 224 RESID = MAX( RESID, SASUM( JU+JL+1, WORK, 1 ) ) 225 END IF 226 40 CONTINUE 227* 228* Compute norm(L*U - A) / ( N * norm(A) * EPS ) 229* 230 IF( ANORM.LE.ZERO ) THEN 231 IF( RESID.NE.ZERO ) 232 $ RESID = ONE / EPS 233 ELSE 234 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 235 END IF 236* 237 RETURN 238* 239* End of SGBT01 240* 241 END 242