1*> \brief \b SGTT01 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 SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 12* LDWORK, RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* INTEGER LDWORK, N 16* REAL RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* REAL D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 21* $ DU2( * ), DUF( * ), RWORK( * ), 22* $ WORK( LDWORK, * ) 23* .. 24* 25* 26*> \par Purpose: 27* ============= 28*> 29*> \verbatim 30*> 31*> SGTT01 reconstructs a tridiagonal matrix A from its LU factorization 32*> and computes the residual 33*> norm(L*U - A) / ( norm(A) * EPS ), 34*> where EPS is the machine epsilon. 35*> \endverbatim 36* 37* Arguments: 38* ========== 39* 40*> \param[in] N 41*> \verbatim 42*> N is INTEGTER 43*> The order of the matrix A. N >= 0. 44*> \endverbatim 45*> 46*> \param[in] DL 47*> \verbatim 48*> DL is REAL array, dimension (N-1) 49*> The (n-1) sub-diagonal elements of A. 50*> \endverbatim 51*> 52*> \param[in] D 53*> \verbatim 54*> D is REAL array, dimension (N) 55*> The diagonal elements of A. 56*> \endverbatim 57*> 58*> \param[in] DU 59*> \verbatim 60*> DU is REAL array, dimension (N-1) 61*> The (n-1) super-diagonal elements of A. 62*> \endverbatim 63*> 64*> \param[in] DLF 65*> \verbatim 66*> DLF is REAL array, dimension (N-1) 67*> The (n-1) multipliers that define the matrix L from the 68*> LU factorization of A. 69*> \endverbatim 70*> 71*> \param[in] DF 72*> \verbatim 73*> DF is REAL array, dimension (N) 74*> The n diagonal elements of the upper triangular matrix U from 75*> the LU factorization of A. 76*> \endverbatim 77*> 78*> \param[in] DUF 79*> \verbatim 80*> DUF is REAL array, dimension (N-1) 81*> The (n-1) elements of the first super-diagonal of U. 82*> \endverbatim 83*> 84*> \param[in] DU2 85*> \verbatim 86*> DU2 is REAL array, dimension (N-2) 87*> The (n-2) elements of the second super-diagonal of U. 88*> \endverbatim 89*> 90*> \param[in] IPIV 91*> \verbatim 92*> IPIV is INTEGER array, dimension (N) 93*> The pivot indices; for 1 <= i <= n, row i of the matrix was 94*> interchanged with row IPIV(i). IPIV(i) will always be either 95*> i or i+1; IPIV(i) = i indicates a row interchange was not 96*> required. 97*> \endverbatim 98*> 99*> \param[out] WORK 100*> \verbatim 101*> WORK is REAL array, dimension (LDWORK,N) 102*> \endverbatim 103*> 104*> \param[in] LDWORK 105*> \verbatim 106*> LDWORK is INTEGER 107*> The leading dimension of the array WORK. LDWORK >= max(1,N). 108*> \endverbatim 109*> 110*> \param[out] RWORK 111*> \verbatim 112*> RWORK is REAL array, dimension (N) 113*> \endverbatim 114*> 115*> \param[out] RESID 116*> \verbatim 117*> RESID is REAL 118*> The scaled residual: norm(L*U - A) / (norm(A) * EPS) 119*> \endverbatim 120* 121* Authors: 122* ======== 123* 124*> \author Univ. of Tennessee 125*> \author Univ. of California Berkeley 126*> \author Univ. of Colorado Denver 127*> \author NAG Ltd. 128* 129*> \date November 2011 130* 131*> \ingroup single_lin 132* 133* ===================================================================== 134 SUBROUTINE SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 135 $ LDWORK, RWORK, RESID ) 136* 137* -- LAPACK test routine (version 3.4.0) -- 138* -- LAPACK is a software package provided by Univ. of Tennessee, -- 139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 140* November 2011 141* 142* .. Scalar Arguments .. 143 INTEGER LDWORK, N 144 REAL RESID 145* .. 146* .. Array Arguments .. 147 INTEGER IPIV( * ) 148 REAL D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 149 $ DU2( * ), DUF( * ), RWORK( * ), 150 $ WORK( LDWORK, * ) 151* .. 152* 153* ===================================================================== 154* 155* .. Parameters .. 156 REAL ONE, ZERO 157 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 158* .. 159* .. Local Scalars .. 160 INTEGER I, IP, J, LASTJ 161 REAL ANORM, EPS, LI 162* .. 163* .. External Functions .. 164 REAL SLAMCH, SLANGT, SLANHS 165 EXTERNAL SLAMCH, SLANGT, SLANHS 166* .. 167* .. Intrinsic Functions .. 168 INTRINSIC MIN 169* .. 170* .. External Subroutines .. 171 EXTERNAL SAXPY, SSWAP 172* .. 173* .. Executable Statements .. 174* 175* Quick return if possible 176* 177 IF( N.LE.0 ) THEN 178 RESID = ZERO 179 RETURN 180 END IF 181* 182 EPS = SLAMCH( 'Epsilon' ) 183* 184* Copy the matrix U to WORK. 185* 186 DO 20 J = 1, N 187 DO 10 I = 1, N 188 WORK( I, J ) = ZERO 189 10 CONTINUE 190 20 CONTINUE 191 DO 30 I = 1, N 192 IF( I.EQ.1 ) THEN 193 WORK( I, I ) = DF( I ) 194 IF( N.GE.2 ) 195 $ WORK( I, I+1 ) = DUF( I ) 196 IF( N.GE.3 ) 197 $ WORK( I, I+2 ) = DU2( I ) 198 ELSE IF( I.EQ.N ) THEN 199 WORK( I, I ) = DF( I ) 200 ELSE 201 WORK( I, I ) = DF( I ) 202 WORK( I, I+1 ) = DUF( I ) 203 IF( I.LT.N-1 ) 204 $ WORK( I, I+2 ) = DU2( I ) 205 END IF 206 30 CONTINUE 207* 208* Multiply on the left by L. 209* 210 LASTJ = N 211 DO 40 I = N - 1, 1, -1 212 LI = DLF( I ) 213 CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK, 214 $ WORK( I+1, I ), LDWORK ) 215 IP = IPIV( I ) 216 IF( IP.EQ.I ) THEN 217 LASTJ = MIN( I+2, N ) 218 ELSE 219 CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ), 220 $ LDWORK ) 221 END IF 222 40 CONTINUE 223* 224* Subtract the matrix A. 225* 226 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 ) 227 IF( N.GT.1 ) THEN 228 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 ) 229 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 ) 230 WORK( N, N ) = WORK( N, N ) - D( N ) 231 DO 50 I = 2, N - 1 232 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 ) 233 WORK( I, I ) = WORK( I, I ) - D( I ) 234 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I ) 235 50 CONTINUE 236 END IF 237* 238* Compute the 1-norm of the tridiagonal matrix A. 239* 240 ANORM = SLANGT( '1', N, DL, D, DU ) 241* 242* Compute the 1-norm of WORK, which is only guaranteed to be 243* upper Hessenberg. 244* 245 RESID = SLANHS( '1', N, WORK, LDWORK, RWORK ) 246* 247* Compute norm(L*U - A) / (norm(A) * EPS) 248* 249 IF( ANORM.LE.ZERO ) THEN 250 IF( RESID.NE.ZERO ) 251 $ RESID = ONE / EPS 252 ELSE 253 RESID = ( RESID / ANORM ) / EPS 254 END IF 255* 256 RETURN 257* 258* End of SGTT01 259* 260 END 261