![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGTT01 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 CGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK, 00012 * LDWORK, RWORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDWORK, N 00016 * REAL RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * REAL RWORK( * ) 00021 * COMPLEX D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 00022 * $ DU2( * ), DUF( * ), WORK( LDWORK, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> CGTT01 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 COMPLEX 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 COMPLEX array, dimension (N) 00055 *> The diagonal elements of A. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] DU 00059 *> \verbatim 00060 *> DU is COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 REAL array, dimension (N) 00113 *> \endverbatim 00114 *> 00115 *> \param[out] RESID 00116 *> \verbatim 00117 *> RESID is REAL 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 complex_lin 00132 * 00133 * ===================================================================== 00134 SUBROUTINE CGTT01( 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 REAL RESID 00145 * .. 00146 * .. Array Arguments .. 00147 INTEGER IPIV( * ) 00148 REAL RWORK( * ) 00149 COMPLEX D( * ), DF( * ), DL( * ), DLF( * ), DU( * ), 00150 $ DU2( * ), DUF( * ), WORK( LDWORK, * ) 00151 * .. 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 REAL ONE, ZERO 00157 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 INTEGER I, IP, J, LASTJ 00161 REAL ANORM, EPS 00162 COMPLEX LI 00163 * .. 00164 * .. External Functions .. 00165 REAL CLANGT, CLANHS, SLAMCH 00166 EXTERNAL CLANGT, CLANHS, SLAMCH 00167 * .. 00168 * .. Intrinsic Functions .. 00169 INTRINSIC MIN 00170 * .. 00171 * .. External Subroutines .. 00172 EXTERNAL CAXPY, CSWAP 00173 * .. 00174 * .. Executable Statements .. 00175 * 00176 * Quick return if possible 00177 * 00178 IF( N.LE.0 ) THEN 00179 RESID = ZERO 00180 RETURN 00181 END IF 00182 * 00183 EPS = SLAMCH( 'Epsilon' ) 00184 * 00185 * Copy the matrix U to WORK. 00186 * 00187 DO 20 J = 1, N 00188 DO 10 I = 1, N 00189 WORK( I, J ) = ZERO 00190 10 CONTINUE 00191 20 CONTINUE 00192 DO 30 I = 1, N 00193 IF( I.EQ.1 ) THEN 00194 WORK( I, I ) = DF( I ) 00195 IF( N.GE.2 ) 00196 $ WORK( I, I+1 ) = DUF( I ) 00197 IF( N.GE.3 ) 00198 $ WORK( I, I+2 ) = DU2( I ) 00199 ELSE IF( I.EQ.N ) THEN 00200 WORK( I, I ) = DF( I ) 00201 ELSE 00202 WORK( I, I ) = DF( I ) 00203 WORK( I, I+1 ) = DUF( I ) 00204 IF( I.LT.N-1 ) 00205 $ WORK( I, I+2 ) = DU2( I ) 00206 END IF 00207 30 CONTINUE 00208 * 00209 * Multiply on the left by L. 00210 * 00211 LASTJ = N 00212 DO 40 I = N - 1, 1, -1 00213 LI = DLF( I ) 00214 CALL CAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK, 00215 $ WORK( I+1, I ), LDWORK ) 00216 IP = IPIV( I ) 00217 IF( IP.EQ.I ) THEN 00218 LASTJ = MIN( I+2, N ) 00219 ELSE 00220 CALL CSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ), 00221 $ LDWORK ) 00222 END IF 00223 40 CONTINUE 00224 * 00225 * Subtract the matrix A. 00226 * 00227 WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 ) 00228 IF( N.GT.1 ) THEN 00229 WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 ) 00230 WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 ) 00231 WORK( N, N ) = WORK( N, N ) - D( N ) 00232 DO 50 I = 2, N - 1 00233 WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 ) 00234 WORK( I, I ) = WORK( I, I ) - D( I ) 00235 WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I ) 00236 50 CONTINUE 00237 END IF 00238 * 00239 * Compute the 1-norm of the tridiagonal matrix A. 00240 * 00241 ANORM = CLANGT( '1', N, DL, D, DU ) 00242 * 00243 * Compute the 1-norm of WORK, which is only guaranteed to be 00244 * upper Hessenberg. 00245 * 00246 RESID = CLANHS( '1', N, WORK, LDWORK, RWORK ) 00247 * 00248 * Compute norm(L*U - A) / (norm(A) * EPS) 00249 * 00250 IF( ANORM.LE.ZERO ) THEN 00251 IF( RESID.NE.ZERO ) 00252 $ RESID = ONE / EPS 00253 ELSE 00254 RESID = ( RESID / ANORM ) / EPS 00255 END IF 00256 * 00257 RETURN 00258 * 00259 * End of CGTT01 00260 * 00261 END