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