![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSYT01 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 SSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, 00012 * RWORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER LDA, LDAFAC, LDC, N 00017 * REAL RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * INTEGER IPIV( * ) 00021 * REAL A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 00022 * $ RWORK( * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> SSYT01 reconstructs a symmetric indefinite matrix A from its 00032 *> block L*D*L' or U*D*U' factorization and computes the residual 00033 *> norm( C - A ) / ( N * norm(A) * EPS ), 00034 *> where C is the reconstructed matrix and EPS is the machine epsilon. 00035 *> \endverbatim 00036 * 00037 * Arguments: 00038 * ========== 00039 * 00040 *> \param[in] UPLO 00041 *> \verbatim 00042 *> UPLO is CHARACTER*1 00043 *> Specifies whether the upper or lower triangular part of the 00044 *> symmetric matrix A is stored: 00045 *> = 'U': Upper triangular 00046 *> = 'L': Lower triangular 00047 *> \endverbatim 00048 *> 00049 *> \param[in] N 00050 *> \verbatim 00051 *> N is INTEGER 00052 *> The number of rows and columns of the matrix A. N >= 0. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] A 00056 *> \verbatim 00057 *> A is REAL array, dimension (LDA,N) 00058 *> The original symmetric matrix A. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] LDA 00062 *> \verbatim 00063 *> LDA is INTEGER 00064 *> The leading dimension of the array A. LDA >= max(1,N) 00065 *> \endverbatim 00066 *> 00067 *> \param[in] AFAC 00068 *> \verbatim 00069 *> AFAC is REAL array, dimension (LDAFAC,N) 00070 *> The factored form of the matrix A. AFAC contains the block 00071 *> diagonal matrix D and the multipliers used to obtain the 00072 *> factor L or U from the block L*D*L' or U*D*U' factorization 00073 *> as computed by SSYTRF. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] LDAFAC 00077 *> \verbatim 00078 *> LDAFAC is INTEGER 00079 *> The leading dimension of the array AFAC. LDAFAC >= max(1,N). 00080 *> \endverbatim 00081 *> 00082 *> \param[in] IPIV 00083 *> \verbatim 00084 *> IPIV is INTEGER array, dimension (N) 00085 *> The pivot indices from SSYTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] C 00089 *> \verbatim 00090 *> C is REAL array, dimension (LDC,N) 00091 *> \endverbatim 00092 *> 00093 *> \param[in] LDC 00094 *> \verbatim 00095 *> LDC is INTEGER 00096 *> The leading dimension of the array C. LDC >= max(1,N). 00097 *> \endverbatim 00098 *> 00099 *> \param[out] RWORK 00100 *> \verbatim 00101 *> RWORK is REAL array, dimension (N) 00102 *> \endverbatim 00103 *> 00104 *> \param[out] RESID 00105 *> \verbatim 00106 *> RESID is REAL 00107 *> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 00108 *> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 00109 *> \endverbatim 00110 * 00111 * Authors: 00112 * ======== 00113 * 00114 *> \author Univ. of Tennessee 00115 *> \author Univ. of California Berkeley 00116 *> \author Univ. of Colorado Denver 00117 *> \author NAG Ltd. 00118 * 00119 *> \date April 2012 00120 * 00121 *> \ingroup single_lin 00122 * 00123 * ===================================================================== 00124 SUBROUTINE SSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, 00125 $ RWORK, RESID ) 00126 * 00127 * -- LAPACK test routine (version 3.4.1) -- 00128 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00130 * April 2012 00131 * 00132 * .. Scalar Arguments .. 00133 CHARACTER UPLO 00134 INTEGER LDA, LDAFAC, LDC, N 00135 REAL RESID 00136 * .. 00137 * .. Array Arguments .. 00138 INTEGER IPIV( * ) 00139 REAL A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 00140 $ RWORK( * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 REAL ZERO, ONE 00147 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00148 * .. 00149 * .. Local Scalars .. 00150 INTEGER I, INFO, J 00151 REAL ANORM, EPS 00152 * .. 00153 * .. External Functions .. 00154 LOGICAL LSAME 00155 REAL SLAMCH, SLANSY 00156 EXTERNAL LSAME, SLAMCH, SLANSY 00157 * .. 00158 * .. External Subroutines .. 00159 EXTERNAL SLASET, SLAVSY 00160 * .. 00161 * .. Intrinsic Functions .. 00162 INTRINSIC REAL 00163 * .. 00164 * .. Executable Statements .. 00165 * 00166 * Quick exit if N = 0. 00167 * 00168 IF( N.LE.0 ) THEN 00169 RESID = ZERO 00170 RETURN 00171 END IF 00172 * 00173 * Determine EPS and the norm of A. 00174 * 00175 EPS = SLAMCH( 'Epsilon' ) 00176 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) 00177 * 00178 * Initialize C to the identity matrix. 00179 * 00180 CALL SLASET( 'Full', N, N, ZERO, ONE, C, LDC ) 00181 * 00182 * Call SLAVSY to form the product D * U' (or D * L' ). 00183 * 00184 CALL SLAVSY( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, LDAFAC, 00185 $ IPIV, C, LDC, INFO ) 00186 * 00187 * Call SLAVSY again to multiply by U (or L ). 00188 * 00189 CALL SLAVSY( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC, 00190 $ IPIV, C, LDC, INFO ) 00191 * 00192 * Compute the difference C - A . 00193 * 00194 IF( LSAME( UPLO, 'U' ) ) THEN 00195 DO 20 J = 1, N 00196 DO 10 I = 1, J 00197 C( I, J ) = C( I, J ) - A( I, J ) 00198 10 CONTINUE 00199 20 CONTINUE 00200 ELSE 00201 DO 40 J = 1, N 00202 DO 30 I = J, N 00203 C( I, J ) = C( I, J ) - A( I, J ) 00204 30 CONTINUE 00205 40 CONTINUE 00206 END IF 00207 * 00208 * Compute norm( C - A ) / ( N * norm(A) * EPS ) 00209 * 00210 RESID = SLANSY( '1', UPLO, N, C, LDC, RWORK ) 00211 * 00212 IF( ANORM.LE.ZERO ) THEN 00213 IF( RESID.NE.ZERO ) 00214 $ RESID = ONE / EPS 00215 ELSE 00216 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 00217 END IF 00218 * 00219 RETURN 00220 * 00221 * End of SSYT01 00222 * 00223 END