1*> \brief \b SPPT01 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 SPPT01( UPLO, N, A, AFAC, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER UPLO 15* INTEGER N 16* REAL RESID 17* .. 18* .. Array Arguments .. 19* REAL A( * ), AFAC( * ), RWORK( * ) 20* .. 21* 22* 23*> \par Purpose: 24* ============= 25*> 26*> \verbatim 27*> 28*> SPPT01 reconstructs a symmetric positive definite packed matrix A 29*> from 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 (N*(N+1)/2) 56*> The original symmetric matrix A, stored as a packed 57*> triangular matrix. 58*> \endverbatim 59*> 60*> \param[in,out] AFAC 61*> \verbatim 62*> AFAC is REAL array, dimension (N*(N+1)/2) 63*> On entry, the factor L or U from the L*L' or U'*U 64*> factorization of A, stored as a packed triangular matrix. 65*> Overwritten with the reconstructed matrix, and then with the 66*> difference L*L' - A (or U'*U - A). 67*> \endverbatim 68*> 69*> \param[out] RWORK 70*> \verbatim 71*> RWORK is REAL array, dimension (N) 72*> \endverbatim 73*> 74*> \param[out] RESID 75*> \verbatim 76*> RESID is REAL 77*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 78*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 79*> \endverbatim 80* 81* Authors: 82* ======== 83* 84*> \author Univ. of Tennessee 85*> \author Univ. of California Berkeley 86*> \author Univ. of Colorado Denver 87*> \author NAG Ltd. 88* 89*> \date November 2011 90* 91*> \ingroup single_lin 92* 93* ===================================================================== 94 SUBROUTINE SPPT01( UPLO, N, A, AFAC, RWORK, RESID ) 95* 96* -- LAPACK test routine (version 3.4.0) -- 97* -- LAPACK is a software package provided by Univ. of Tennessee, -- 98* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 99* November 2011 100* 101* .. Scalar Arguments .. 102 CHARACTER UPLO 103 INTEGER N 104 REAL RESID 105* .. 106* .. Array Arguments .. 107 REAL A( * ), AFAC( * ), RWORK( * ) 108* .. 109* 110* ===================================================================== 111* 112* .. Parameters .. 113 REAL ZERO, ONE 114 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 115* .. 116* .. Local Scalars .. 117 INTEGER I, K, KC, NPP 118 REAL ANORM, EPS, T 119* .. 120* .. External Functions .. 121 LOGICAL LSAME 122 REAL SDOT, SLAMCH, SLANSP 123 EXTERNAL LSAME, SDOT, SLAMCH, SLANSP 124* .. 125* .. External Subroutines .. 126 EXTERNAL SSCAL, SSPR, STPMV 127* .. 128* .. Intrinsic Functions .. 129 INTRINSIC REAL 130* .. 131* .. Executable Statements .. 132* 133* Quick exit if N = 0 134* 135 IF( N.LE.0 ) THEN 136 RESID = ZERO 137 RETURN 138 END IF 139* 140* Exit with RESID = 1/EPS if ANORM = 0. 141* 142 EPS = SLAMCH( 'Epsilon' ) 143 ANORM = SLANSP( '1', UPLO, N, A, RWORK ) 144 IF( ANORM.LE.ZERO ) THEN 145 RESID = ONE / EPS 146 RETURN 147 END IF 148* 149* Compute the product U'*U, overwriting U. 150* 151 IF( LSAME( UPLO, 'U' ) ) THEN 152 KC = ( N*( N-1 ) ) / 2 + 1 153 DO 10 K = N, 1, -1 154* 155* Compute the (K,K) element of the result. 156* 157 T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 ) 158 AFAC( KC+K-1 ) = T 159* 160* Compute the rest of column K. 161* 162 IF( K.GT.1 ) THEN 163 CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 164 $ AFAC( KC ), 1 ) 165 KC = KC - ( K-1 ) 166 END IF 167 10 CONTINUE 168* 169* Compute the product L*L', overwriting L. 170* 171 ELSE 172 KC = ( N*( N+1 ) ) / 2 173 DO 20 K = N, 1, -1 174* 175* Add a multiple of column K of the factor L to each of 176* columns K+1 through N. 177* 178 IF( K.LT.N ) 179 $ CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1, 180 $ AFAC( KC+N-K+1 ) ) 181* 182* Scale column K by the diagonal element. 183* 184 T = AFAC( KC ) 185 CALL SSCAL( N-K+1, T, AFAC( KC ), 1 ) 186* 187 KC = KC - ( N-K+2 ) 188 20 CONTINUE 189 END IF 190* 191* Compute the difference L*L' - A (or U'*U - A). 192* 193 NPP = N*( N+1 ) / 2 194 DO 30 I = 1, NPP 195 AFAC( I ) = AFAC( I ) - A( I ) 196 30 CONTINUE 197* 198* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 199* 200 RESID = SLANSP( '1', UPLO, N, AFAC, RWORK ) 201* 202 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 203* 204 RETURN 205* 206* End of SPPT01 207* 208 END 209