![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPTT01 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 CPTT01( N, D, E, DF, EF, WORK, RESID ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER N 00015 * REAL RESID 00016 * .. 00017 * .. Array Arguments .. 00018 * REAL D( * ), DF( * ) 00019 * COMPLEX E( * ), EF( * ), WORK( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> CPTT01 reconstructs a tridiagonal matrix A from its L*D*L' 00029 *> factorization and computes the residual 00030 *> norm(L*D*L' - A) / ( n * norm(A) * EPS ), 00031 *> where EPS is the machine epsilon. 00032 *> \endverbatim 00033 * 00034 * Arguments: 00035 * ========== 00036 * 00037 *> \param[in] N 00038 *> \verbatim 00039 *> N is INTEGTER 00040 *> The order of the matrix A. 00041 *> \endverbatim 00042 *> 00043 *> \param[in] D 00044 *> \verbatim 00045 *> D is REAL array, dimension (N) 00046 *> The n diagonal elements of the tridiagonal matrix A. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] E 00050 *> \verbatim 00051 *> E is COMPLEX array, dimension (N-1) 00052 *> The (n-1) subdiagonal elements of the tridiagonal matrix A. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] DF 00056 *> \verbatim 00057 *> DF is REAL array, dimension (N) 00058 *> The n diagonal elements of the factor L from the L*D*L' 00059 *> factorization of A. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] EF 00063 *> \verbatim 00064 *> EF is COMPLEX array, dimension (N-1) 00065 *> The (n-1) subdiagonal elements of the factor L from the 00066 *> L*D*L' factorization of A. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] WORK 00070 *> \verbatim 00071 *> WORK is COMPLEX array, dimension (2*N) 00072 *> \endverbatim 00073 *> 00074 *> \param[out] RESID 00075 *> \verbatim 00076 *> RESID is REAL 00077 *> norm(L*D*L' - A) / (n * norm(A) * EPS) 00078 *> \endverbatim 00079 * 00080 * Authors: 00081 * ======== 00082 * 00083 *> \author Univ. of Tennessee 00084 *> \author Univ. of California Berkeley 00085 *> \author Univ. of Colorado Denver 00086 *> \author NAG Ltd. 00087 * 00088 *> \date November 2011 00089 * 00090 *> \ingroup complex_lin 00091 * 00092 * ===================================================================== 00093 SUBROUTINE CPTT01( N, D, E, DF, EF, WORK, RESID ) 00094 * 00095 * -- LAPACK test routine (version 3.4.0) -- 00096 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00097 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00098 * November 2011 00099 * 00100 * .. Scalar Arguments .. 00101 INTEGER N 00102 REAL RESID 00103 * .. 00104 * .. Array Arguments .. 00105 REAL D( * ), DF( * ) 00106 COMPLEX E( * ), EF( * ), WORK( * ) 00107 * .. 00108 * 00109 * ===================================================================== 00110 * 00111 * .. Parameters .. 00112 REAL ONE, ZERO 00113 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00114 * .. 00115 * .. Local Scalars .. 00116 INTEGER I 00117 REAL ANORM, EPS 00118 COMPLEX DE 00119 * .. 00120 * .. External Functions .. 00121 REAL SLAMCH 00122 EXTERNAL SLAMCH 00123 * .. 00124 * .. Intrinsic Functions .. 00125 INTRINSIC ABS, CONJG, MAX, REAL 00126 * .. 00127 * .. Executable Statements .. 00128 * 00129 * Quick return if possible 00130 * 00131 IF( N.LE.0 ) THEN 00132 RESID = ZERO 00133 RETURN 00134 END IF 00135 * 00136 EPS = SLAMCH( 'Epsilon' ) 00137 * 00138 * Construct the difference L*D*L' - A. 00139 * 00140 WORK( 1 ) = DF( 1 ) - D( 1 ) 00141 DO 10 I = 1, N - 1 00142 DE = DF( I )*EF( I ) 00143 WORK( N+I ) = DE - E( I ) 00144 WORK( 1+I ) = DE*CONJG( EF( I ) ) + DF( I+1 ) - D( I+1 ) 00145 10 CONTINUE 00146 * 00147 * Compute the 1-norms of the tridiagonal matrices A and WORK. 00148 * 00149 IF( N.EQ.1 ) THEN 00150 ANORM = D( 1 ) 00151 RESID = ABS( WORK( 1 ) ) 00152 ELSE 00153 ANORM = MAX( D( 1 )+ABS( E( 1 ) ), D( N )+ABS( E( N-1 ) ) ) 00154 RESID = MAX( ABS( WORK( 1 ) )+ABS( WORK( N+1 ) ), 00155 $ ABS( WORK( N ) )+ABS( WORK( 2*N-1 ) ) ) 00156 DO 20 I = 2, N - 1 00157 ANORM = MAX( ANORM, D( I )+ABS( E( I ) )+ABS( E( I-1 ) ) ) 00158 RESID = MAX( RESID, ABS( WORK( I ) )+ABS( WORK( N+I-1 ) )+ 00159 $ ABS( WORK( N+I ) ) ) 00160 20 CONTINUE 00161 END IF 00162 * 00163 * Compute norm(L*D*L' - A) / (n * norm(A) * EPS) 00164 * 00165 IF( ANORM.LE.ZERO ) THEN 00166 IF( RESID.NE.ZERO ) 00167 $ RESID = ONE / EPS 00168 ELSE 00169 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 00170 END IF 00171 * 00172 RETURN 00173 * 00174 * End of CPTT01 00175 * 00176 END