1*> \brief \b SPOT01 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER UPLO 15* INTEGER LDA, LDAFAC, N 16* REAL RESID 17* .. 18* .. Array Arguments .. 19* REAL A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * ) 20* .. 21* 22* 23*> \par Purpose: 24* ============= 25*> 26*> \verbatim 27*> 28*> SPOT01 reconstructs a symmetric positive definite matrix A from 29*> its L*L' or U'*U factorization and computes the residual 30*> norm( L*L' - A ) / ( N * norm(A) * EPS ) or 31*> norm( U'*U - A ) / ( N * norm(A) * EPS ), 32*> where EPS is the machine epsilon. 33*> \endverbatim 34* 35* Arguments: 36* ========== 37* 38*> \param[in] UPLO 39*> \verbatim 40*> UPLO is CHARACTER*1 41*> Specifies whether the upper or lower triangular part of the 42*> symmetric matrix A is stored: 43*> = 'U': Upper triangular 44*> = 'L': Lower triangular 45*> \endverbatim 46*> 47*> \param[in] N 48*> \verbatim 49*> N is INTEGER 50*> The number of rows and columns of the matrix A. N >= 0. 51*> \endverbatim 52*> 53*> \param[in] A 54*> \verbatim 55*> A is REAL array, dimension (LDA,N) 56*> The original symmetric matrix A. 57*> \endverbatim 58*> 59*> \param[in] LDA 60*> \verbatim 61*> LDA is INTEGER 62*> The leading dimension of the array A. LDA >= max(1,N) 63*> \endverbatim 64*> 65*> \param[in,out] AFAC 66*> \verbatim 67*> AFAC is REAL array, dimension (LDAFAC,N) 68*> On entry, the factor L or U from the L * L**T or U**T * U 69*> factorization of A. 70*> Overwritten with the reconstructed matrix, and then with 71*> the difference L * L**T - A (or U**T * U - A). 72*> \endverbatim 73*> 74*> \param[in] LDAFAC 75*> \verbatim 76*> LDAFAC is INTEGER 77*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). 78*> \endverbatim 79*> 80*> \param[out] RWORK 81*> \verbatim 82*> RWORK is REAL array, dimension (N) 83*> \endverbatim 84*> 85*> \param[out] RESID 86*> \verbatim 87*> RESID is REAL 88*> If UPLO = 'L', norm(L * L**T - A) / ( N * norm(A) * EPS ) 89*> If UPLO = 'U', norm(U**T * U - A) / ( N * norm(A) * EPS ) 90*> \endverbatim 91* 92* Authors: 93* ======== 94* 95*> \author Univ. of Tennessee 96*> \author Univ. of California Berkeley 97*> \author Univ. of Colorado Denver 98*> \author NAG Ltd. 99* 100*> \ingroup single_lin 101* 102* ===================================================================== 103 SUBROUTINE SPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID ) 104* 105* -- LAPACK test routine -- 106* -- LAPACK is a software package provided by Univ. of Tennessee, -- 107* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 108* 109* .. Scalar Arguments .. 110 CHARACTER UPLO 111 INTEGER LDA, LDAFAC, N 112 REAL RESID 113* .. 114* .. Array Arguments .. 115 REAL A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * ) 116* .. 117* 118* ===================================================================== 119* 120* .. Parameters .. 121 REAL ZERO, ONE 122 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 123* .. 124* .. Local Scalars .. 125 INTEGER I, J, K 126 REAL ANORM, EPS, T 127* .. 128* .. External Functions .. 129 LOGICAL LSAME 130 REAL SDOT, SLAMCH, SLANSY 131 EXTERNAL LSAME, SDOT, SLAMCH, SLANSY 132* .. 133* .. External Subroutines .. 134 EXTERNAL SSCAL, SSYR, STRMV 135* .. 136* .. Intrinsic Functions .. 137 INTRINSIC REAL 138* .. 139* .. Executable Statements .. 140* 141* Quick exit if N = 0. 142* 143 IF( N.LE.0 ) THEN 144 RESID = ZERO 145 RETURN 146 END IF 147* 148* Exit with RESID = 1/EPS if ANORM = 0. 149* 150 EPS = SLAMCH( 'Epsilon' ) 151 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) 152 IF( ANORM.LE.ZERO ) THEN 153 RESID = ONE / EPS 154 RETURN 155 END IF 156* 157* Compute the product U**T * U, overwriting U. 158* 159 IF( LSAME( UPLO, 'U' ) ) THEN 160 DO 10 K = N, 1, -1 161* 162* Compute the (K,K) element of the result. 163* 164 T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 165 AFAC( K, K ) = T 166* 167* Compute the rest of column K. 168* 169 CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 170 $ LDAFAC, AFAC( 1, K ), 1 ) 171* 172 10 CONTINUE 173* 174* Compute the product L * L**T, overwriting L. 175* 176 ELSE 177 DO 20 K = N, 1, -1 178* 179* Add a multiple of column K of the factor L to each of 180* columns K+1 through N. 181* 182 IF( K+1.LE.N ) 183 $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 184 $ AFAC( K+1, K+1 ), LDAFAC ) 185* 186* Scale column K by the diagonal element. 187* 188 T = AFAC( K, K ) 189 CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 ) 190* 191 20 CONTINUE 192 END IF 193* 194* Compute the difference L * L**T - A (or U**T * U - A). 195* 196 IF( LSAME( UPLO, 'U' ) ) THEN 197 DO 40 J = 1, N 198 DO 30 I = 1, J 199 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 200 30 CONTINUE 201 40 CONTINUE 202 ELSE 203 DO 60 J = 1, N 204 DO 50 I = J, N 205 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 206 50 CONTINUE 207 60 CONTINUE 208 END IF 209* 210* Compute norm(L*U - A) / ( N * norm(A) * EPS ) 211* 212 RESID = SLANSY( '1', UPLO, N, AFAC, LDAFAC, RWORK ) 213* 214 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 215* 216 RETURN 217* 218* End of SPOT01 219* 220 END 221