![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGTT01 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 00012 * LDWORK, RWORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDWORK, N 00016 * DOUBLE PRECISION RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * DOUBLE PRECISION D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 00021 * $ DU2( * ), DUF( * ), RWORK( * ), 00022 * $ WORK( LDWORK, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> DGTT01 reconstructs a tridiagonal matrix A from its LU factorization 00032 *> and computes the residual 00033 *> norm(L*U - A) / ( norm(A) * EPS ), 00034 *> where EPS is the machine epsilon. 00035 *> \endverbatim 00036 * 00037 * Arguments: 00038 * ========== 00039 * 00040 *> \param[in] N 00041 *> \verbatim 00042 *> N is INTEGTER 00043 *> The order of the matrix A. N >= 0. 00044 *> \endverbatim 00045 *> 00046 *> \param[in] DL 00047 *> \verbatim 00048 *> DL is DOUBLE PRECISION array, dimension (N-1) 00049 *> The (n-1) sub-diagonal elements of A. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] D 00053 *> \verbatim 00054 *> D is DOUBLE PRECISION array, dimension (N) 00055 *> The diagonal elements of A. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] DU 00059 *> \verbatim 00060 *> DU is DOUBLE PRECISION array, dimension (N-1) 00061 *> The (n-1) super-diagonal elements of A. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] DLF 00065 *> \verbatim 00066 *> DLF is DOUBLE PRECISION array, dimension (N-1) 00067 *> The (n-1) multipliers that define the matrix L from the 00068 *> LU factorization of A. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] DF 00072 *> \verbatim 00073 *> DF is DOUBLE PRECISION array, dimension (N) 00074 *> The n diagonal elements of the upper triangular matrix U from 00075 *> the LU factorization of A. 00076 *> \endverbatim 00077 *> 00078 *> \param[in] DUF 00079 *> \verbatim 00080 *> DUF is DOUBLE PRECISION array, dimension (N-1) 00081 *> The (n-1) elements of the first super-diagonal of U. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] DU2 00085 *> \verbatim 00086 *> DU2 is DOUBLE PRECISION array, dimension (N-2) 00087 *> The (n-2) elements of the second super-diagonal of U. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] IPIV 00091 *> \verbatim 00092 *> IPIV is INTEGER array, dimension (N) 00093 *> The pivot indices; for 1 <= i <= n, row i of the matrix was 00094 *> interchanged with row IPIV(i). IPIV(i) will always be either 00095 *> i or i+1; IPIV(i) = i indicates a row interchange was not 00096 *> required. 00097 *> \endverbatim 00098 *> 00099 *> \param[out] WORK 00100 *> \verbatim 00101 *> WORK is DOUBLE PRECISION array, dimension (LDWORK,N) 00102 *> \endverbatim 00103 *> 00104 *> \param[in] LDWORK 00105 *> \verbatim 00106 *> LDWORK is INTEGER 00107 *> The leading dimension of the array WORK. LDWORK >= max(1,N). 00108 *> \endverbatim 00109 *> 00110 *> \param[out] RWORK 00111 *> \verbatim 00112 *> RWORK is DOUBLE PRECISION array, dimension (N) 00113 *> \endverbatim 00114 *> 00115 *> \param[out] RESID 00116 *> \verbatim 00117 *> RESID is DOUBLE PRECISION 00118 *> The scaled residual: norm(L*U - A) / (norm(A) * EPS) 00119 *> \endverbatim 00120 * 00121 * Authors: 00122 * ======== 00123 * 00124 *> \author Univ. of Tennessee 00125 *> \author Univ. of California Berkeley 00126 *> \author Univ. of Colorado Denver 00127 *> \author NAG Ltd. 00128 * 00129 *> \date November 2011 00130 * 00131 *> \ingroup double_lin 00132 * 00133 * ===================================================================== 00134 SUBROUTINE DGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 00135 $ LDWORK, RWORK, RESID ) 00136 * 00137 * -- LAPACK test routine (version 3.4.0) -- 00138 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00140 * November 2011 00141 * 00142 * .. Scalar Arguments .. 00143 INTEGER LDWORK, N 00144 DOUBLE PRECISION RESID 00145 * .. 00146 * .. Array Arguments .. 00147 INTEGER IPIV( * ) 00148 DOUBLE PRECISION D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 00149 $ DU2( * ), DUF( * ), RWORK( * ), 00150 $ WORK( LDWORK, * ) 00151 * .. 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 DOUBLE PRECISION ONE, ZERO 00157 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 INTEGER I, IP, J, LASTJ 00161 DOUBLE PRECISION ANORM, EPS, LI 00162 * .. 00163 * .. External Functions .. 00164 DOUBLE PRECISION DLAMCH, DLANGT, DLANHS 00165 EXTERNAL DLAMCH, DLANGT, DLANHS 00166 * .. 00167 * .. Intrinsic Functions .. 00168 INTRINSIC MIN 00169 * .. 00170 * .. External Subroutines .. 00171 EXTERNAL DAXPY, DSWAP 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Quick return if possible 00176 * 00177 IF( N.LE.0 ) THEN 00178 RESID = ZERO 00179 RETURN 00180 END IF 00181 * 00182 EPS = DLAMCH( 'Epsilon' ) 00183 * 00184 * Copy the matrix U to WORK. 00185 * 00186 DO 20 J = 1, N 00187 DO 10 I = 1, N 00188 WORK( I, J ) = ZERO 00189 10 CONTINUE 00190 20 CONTINUE 00191 DO 30 I = 1, N 00192 IF( I.EQ.1 ) THEN 00193 WORK( I, I ) = DF( I ) 00194 IF( N.GE.2 ) 00195 $ WORK( I, I+1 ) = DUF( I ) 00196 IF( N.GE.3 ) 00197 $ WORK( I, I+2 ) = DU2( I ) 00198 ELSE IF( I.EQ.N ) THEN 00199 WORK( I, I ) = DF( I ) 00200 ELSE 00201 WORK( I, I ) = DF( I ) 00202 WORK( I, I+1 ) = DUF( I ) 00203 IF( I.LT.N-1 ) 00204 $ WORK( I, I+2 ) = DU2( I ) 00205 END IF 00206 30 CONTINUE 00207 * 00208 * Multiply on the left by L. 00209 * 00210 LASTJ = N 00211 DO 40 I = N - 1, 1, -1 00212 LI = DLF( I ) 00213 CALL DAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK, 00214 $ WORK( I+1, I ), LDWORK ) 00215 IP = IPIV( I ) 00216 IF( IP.EQ.I ) THEN 00217 LASTJ = MIN( I+2, N ) 00218 ELSE 00219 CALL DSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ), 00220 $ LDWORK ) 00221 END IF 00222 40 CONTINUE 00223 * 00224 * Subtract the matrix A. 00225 * 00226 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 ) 00227 IF( N.GT.1 ) THEN 00228 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 ) 00229 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 ) 00230 WORK( N, N ) = WORK( N, N ) - D( N ) 00231 DO 50 I = 2, N - 1 00232 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 ) 00233 WORK( I, I ) = WORK( I, I ) - D( I ) 00234 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I ) 00235 50 CONTINUE 00236 END IF 00237 * 00238 * Compute the 1-norm of the tridiagonal matrix A. 00239 * 00240 ANORM = DLANGT( '1', N, DL, D, DU ) 00241 * 00242 * Compute the 1-norm of WORK, which is only guaranteed to be 00243 * upper Hessenberg. 00244 * 00245 RESID = DLANHS( '1', N, WORK, LDWORK, RWORK ) 00246 * 00247 * Compute norm(L*U - A) / (norm(A) * EPS) 00248 * 00249 IF( ANORM.LE.ZERO ) THEN 00250 IF( RESID.NE.ZERO ) 00251 $ RESID = ONE / EPS 00252 ELSE 00253 RESID = ( RESID / ANORM ) / EPS 00254 END IF 00255 * 00256 RETURN 00257 * 00258 * End of DGTT01 00259 * 00260 END