1*> \brief \b ZGTT01 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 ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 12* LDWORK, RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* INTEGER LDWORK, N 16* DOUBLE PRECISION RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* DOUBLE PRECISION RWORK( * ) 21* COMPLEX*16 D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 22* $ DU2( * ), DUF( * ), WORK( LDWORK, * ) 23* .. 24* 25* 26*> \par Purpose: 27* ============= 28*> 29*> \verbatim 30*> 31*> ZGTT01 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 COMPLEX*16 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 COMPLEX*16 array, dimension (N) 55*> The diagonal elements of A. 56*> \endverbatim 57*> 58*> \param[in] DU 59*> \verbatim 60*> DU is COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 DOUBLE PRECISION array, dimension (N) 113*> \endverbatim 114*> 115*> \param[out] RESID 116*> \verbatim 117*> RESID is DOUBLE PRECISION 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 December 2016 130* 131*> \ingroup complex16_lin 132* 133* ===================================================================== 134 SUBROUTINE ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 135 $ LDWORK, RWORK, RESID ) 136* 137* -- LAPACK test routine (version 3.7.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* December 2016 141* 142* .. Scalar Arguments .. 143 INTEGER LDWORK, N 144 DOUBLE PRECISION RESID 145* .. 146* .. Array Arguments .. 147 INTEGER IPIV( * ) 148 DOUBLE PRECISION RWORK( * ) 149 COMPLEX*16 D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 150 $ DU2( * ), DUF( * ), WORK( LDWORK, * ) 151* .. 152* 153* ===================================================================== 154* 155* .. Parameters .. 156 DOUBLE PRECISION ONE, ZERO 157 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 158* .. 159* .. Local Scalars .. 160 INTEGER I, IP, J, LASTJ 161 DOUBLE PRECISION ANORM, EPS 162 COMPLEX*16 LI 163* .. 164* .. External Functions .. 165 DOUBLE PRECISION DLAMCH, ZLANGT, ZLANHS 166 EXTERNAL DLAMCH, ZLANGT, ZLANHS 167* .. 168* .. Intrinsic Functions .. 169 INTRINSIC MIN 170* .. 171* .. External Subroutines .. 172 EXTERNAL ZAXPY, ZSWAP 173* .. 174* .. Executable Statements .. 175* 176* Quick return if possible 177* 178 IF( N.LE.0 ) THEN 179 RESID = ZERO 180 RETURN 181 END IF 182* 183 EPS = DLAMCH( 'Epsilon' ) 184* 185* Copy the matrix U to WORK. 186* 187 DO 20 J = 1, N 188 DO 10 I = 1, N 189 WORK( I, J ) = ZERO 190 10 CONTINUE 191 20 CONTINUE 192 DO 30 I = 1, N 193 IF( I.EQ.1 ) THEN 194 WORK( I, I ) = DF( I ) 195 IF( N.GE.2 ) 196 $ WORK( I, I+1 ) = DUF( I ) 197 IF( N.GE.3 ) 198 $ WORK( I, I+2 ) = DU2( I ) 199 ELSE IF( I.EQ.N ) THEN 200 WORK( I, I ) = DF( I ) 201 ELSE 202 WORK( I, I ) = DF( I ) 203 WORK( I, I+1 ) = DUF( I ) 204 IF( I.LT.N-1 ) 205 $ WORK( I, I+2 ) = DU2( I ) 206 END IF 207 30 CONTINUE 208* 209* Multiply on the left by L. 210* 211 LASTJ = N 212 DO 40 I = N - 1, 1, -1 213 LI = DLF( I ) 214 CALL ZAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK, 215 $ WORK( I+1, I ), LDWORK ) 216 IP = IPIV( I ) 217 IF( IP.EQ.I ) THEN 218 LASTJ = MIN( I+2, N ) 219 ELSE 220 CALL ZSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ), 221 $ LDWORK ) 222 END IF 223 40 CONTINUE 224* 225* Subtract the matrix A. 226* 227 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 ) 228 IF( N.GT.1 ) THEN 229 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 ) 230 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 ) 231 WORK( N, N ) = WORK( N, N ) - D( N ) 232 DO 50 I = 2, N - 1 233 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 ) 234 WORK( I, I ) = WORK( I, I ) - D( I ) 235 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I ) 236 50 CONTINUE 237 END IF 238* 239* Compute the 1-norm of the tridiagonal matrix A. 240* 241 ANORM = ZLANGT( '1', N, DL, D, DU ) 242* 243* Compute the 1-norm of WORK, which is only guaranteed to be 244* upper Hessenberg. 245* 246 RESID = ZLANHS( '1', N, WORK, LDWORK, RWORK ) 247* 248* Compute norm(L*U - A) / (norm(A) * EPS) 249* 250 IF( ANORM.LE.ZERO ) THEN 251 IF( RESID.NE.ZERO ) 252 $ RESID = ONE / EPS 253 ELSE 254 RESID = ( RESID / ANORM ) / EPS 255 END IF 256* 257 RETURN 258* 259* End of ZGTT01 260* 261 END 262