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*> \date December 2016 122* 123*> \ingroup single_lin 124* 125* ===================================================================== 126 SUBROUTINE SGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 127 $ RESID ) 128* 129* -- LAPACK test routine (version 3.7.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* December 2016 133* 134* .. Scalar Arguments .. 135 INTEGER KL, KU, LDA, LDAFAC, M, N 136 REAL RESID 137* .. 138* .. Array Arguments .. 139 INTEGER IPIV( * ) 140 REAL A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 141* .. 142* 143* ===================================================================== 144* 145* .. Parameters .. 146 REAL ZERO, ONE 147 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 148* .. 149* .. Local Scalars .. 150 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ 151 REAL ANORM, EPS, T 152* .. 153* .. External Functions .. 154 REAL SASUM, SLAMCH 155 EXTERNAL SASUM, SLAMCH 156* .. 157* .. External Subroutines .. 158 EXTERNAL SAXPY, SCOPY 159* .. 160* .. Intrinsic Functions .. 161 INTRINSIC MAX, MIN, REAL 162* .. 163* .. Executable Statements .. 164* 165* Quick exit if M = 0 or N = 0. 166* 167 RESID = ZERO 168 IF( M.LE.0 .OR. N.LE.0 ) 169 $ RETURN 170* 171* Determine EPS and the norm of A. 172* 173 EPS = SLAMCH( 'Epsilon' ) 174 KD = KU + 1 175 ANORM = ZERO 176 DO 10 J = 1, N 177 I1 = MAX( KD+1-J, 1 ) 178 I2 = MIN( KD+M-J, KL+KD ) 179 IF( I2.GE.I1 ) 180 $ ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) ) 181 10 CONTINUE 182* 183* Compute one column at a time of L*U - A. 184* 185 KD = KL + KU + 1 186 DO 40 J = 1, N 187* 188* Copy the J-th column of U to WORK. 189* 190 JU = MIN( KL+KU, J-1 ) 191 JL = MIN( KL, M-J ) 192 LENJ = MIN( M, J ) - J + JU + 1 193 IF( LENJ.GT.0 ) THEN 194 CALL SCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 195 DO 20 I = LENJ + 1, JU + JL + 1 196 WORK( I ) = ZERO 197 20 CONTINUE 198* 199* Multiply by the unit lower triangular matrix L. Note that L 200* is stored as a product of transformations and permutations. 201* 202 DO 30 I = MIN( M-1, J ), J - JU, -1 203 IL = MIN( KL, M-I ) 204 IF( IL.GT.0 ) THEN 205 IW = I - J + JU + 1 206 T = WORK( IW ) 207 CALL SAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 208 $ 1 ) 209 IP = IPIV( I ) 210 IF( I.NE.IP ) THEN 211 IP = IP - J + JU + 1 212 WORK( IW ) = WORK( IP ) 213 WORK( IP ) = T 214 END IF 215 END IF 216 30 CONTINUE 217* 218* Subtract the corresponding column of A. 219* 220 JUA = MIN( JU, KU ) 221 IF( JUA+JL+1.GT.0 ) 222 $ CALL SAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1, 223 $ WORK( JU+1-JUA ), 1 ) 224* 225* Compute the 1-norm of the column. 226* 227 RESID = MAX( RESID, SASUM( JU+JL+1, WORK, 1 ) ) 228 END IF 229 40 CONTINUE 230* 231* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 232* 233 IF( ANORM.LE.ZERO ) THEN 234 IF( RESID.NE.ZERO ) 235 $ RESID = ONE / EPS 236 ELSE 237 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 238 END IF 239* 240 RETURN 241* 242* End of SGBT01 243* 244 END 245