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*> \ingroup single_lin 90* 91* ===================================================================== 92 SUBROUTINE SPPT01( UPLO, N, A, AFAC, RWORK, RESID ) 93* 94* -- LAPACK test routine -- 95* -- LAPACK is a software package provided by Univ. of Tennessee, -- 96* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 97* 98* .. Scalar Arguments .. 99 CHARACTER UPLO 100 INTEGER N 101 REAL RESID 102* .. 103* .. Array Arguments .. 104 REAL A( * ), AFAC( * ), RWORK( * ) 105* .. 106* 107* ===================================================================== 108* 109* .. Parameters .. 110 REAL ZERO, ONE 111 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 112* .. 113* .. Local Scalars .. 114 INTEGER I, K, KC, NPP 115 REAL ANORM, EPS, T 116* .. 117* .. External Functions .. 118 LOGICAL LSAME 119 REAL SDOT, SLAMCH, SLANSP 120 EXTERNAL LSAME, SDOT, SLAMCH, SLANSP 121* .. 122* .. External Subroutines .. 123 EXTERNAL SSCAL, SSPR, STPMV 124* .. 125* .. Intrinsic Functions .. 126 INTRINSIC REAL 127* .. 128* .. Executable Statements .. 129* 130* Quick exit if N = 0 131* 132 IF( N.LE.0 ) THEN 133 RESID = ZERO 134 RETURN 135 END IF 136* 137* Exit with RESID = 1/EPS if ANORM = 0. 138* 139 EPS = SLAMCH( 'Epsilon' ) 140 ANORM = SLANSP( '1', UPLO, N, A, RWORK ) 141 IF( ANORM.LE.ZERO ) THEN 142 RESID = ONE / EPS 143 RETURN 144 END IF 145* 146* Compute the product U'*U, overwriting U. 147* 148 IF( LSAME( UPLO, 'U' ) ) THEN 149 KC = ( N*( N-1 ) ) / 2 + 1 150 DO 10 K = N, 1, -1 151* 152* Compute the (K,K) element of the result. 153* 154 T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 ) 155 AFAC( KC+K-1 ) = T 156* 157* Compute the rest of column K. 158* 159 IF( K.GT.1 ) THEN 160 CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 161 $ AFAC( KC ), 1 ) 162 KC = KC - ( K-1 ) 163 END IF 164 10 CONTINUE 165* 166* Compute the product L*L', overwriting L. 167* 168 ELSE 169 KC = ( N*( N+1 ) ) / 2 170 DO 20 K = N, 1, -1 171* 172* Add a multiple of column K of the factor L to each of 173* columns K+1 through N. 174* 175 IF( K.LT.N ) 176 $ CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1, 177 $ AFAC( KC+N-K+1 ) ) 178* 179* Scale column K by the diagonal element. 180* 181 T = AFAC( KC ) 182 CALL SSCAL( N-K+1, T, AFAC( KC ), 1 ) 183* 184 KC = KC - ( N-K+2 ) 185 20 CONTINUE 186 END IF 187* 188* Compute the difference L*L' - A (or U'*U - A). 189* 190 NPP = N*( N+1 ) / 2 191 DO 30 I = 1, NPP 192 AFAC( I ) = AFAC( I ) - A( I ) 193 30 CONTINUE 194* 195* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 196* 197 RESID = SLANSP( '1', UPLO, N, AFAC, RWORK ) 198* 199 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 200* 201 RETURN 202* 203* End of SPPT01 204* 205 END 206