1*> \brief \b SPST01 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 SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 12* PIV, RWORK, RESID, RANK ) 13* 14* .. Scalar Arguments .. 15* REAL RESID 16* INTEGER LDA, LDAFAC, LDPERM, N, RANK 17* CHARACTER UPLO 18* .. 19* .. Array Arguments .. 20* REAL A( LDA, * ), AFAC( LDAFAC, * ), 21* $ PERM( LDPERM, * ), RWORK( * ) 22* INTEGER PIV( * ) 23* .. 24* 25* 26*> \par Purpose: 27* ============= 28*> 29*> \verbatim 30*> 31*> SPST01 reconstructs a symmetric positive semidefinite matrix A 32*> from its L or U factors and the permutation matrix P and computes 33*> the residual 34*> norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or 35*> norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ), 36*> where EPS is the machine epsilon. 37*> \endverbatim 38* 39* Arguments: 40* ========== 41* 42*> \param[in] UPLO 43*> \verbatim 44*> UPLO is CHARACTER*1 45*> Specifies whether the upper or lower triangular part of the 46*> symmetric matrix A is stored: 47*> = 'U': Upper triangular 48*> = 'L': Lower triangular 49*> \endverbatim 50*> 51*> \param[in] N 52*> \verbatim 53*> N is INTEGER 54*> The number of rows and columns of the matrix A. N >= 0. 55*> \endverbatim 56*> 57*> \param[in] A 58*> \verbatim 59*> A is REAL array, dimension (LDA,N) 60*> The original symmetric matrix A. 61*> \endverbatim 62*> 63*> \param[in] LDA 64*> \verbatim 65*> LDA is INTEGER 66*> The leading dimension of the array A. LDA >= max(1,N) 67*> \endverbatim 68*> 69*> \param[in] AFAC 70*> \verbatim 71*> AFAC is REAL array, dimension (LDAFAC,N) 72*> The factor L or U from the L*L' or U'*U 73*> factorization of A. 74*> \endverbatim 75*> 76*> \param[in] LDAFAC 77*> \verbatim 78*> LDAFAC is INTEGER 79*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). 80*> \endverbatim 81*> 82*> \param[out] PERM 83*> \verbatim 84*> PERM is REAL array, dimension (LDPERM,N) 85*> Overwritten with the reconstructed matrix, and then with the 86*> difference P*L*L'*P' - A (or P*U'*U*P' - A) 87*> \endverbatim 88*> 89*> \param[in] LDPERM 90*> \verbatim 91*> LDPERM is INTEGER 92*> The leading dimension of the array PERM. 93*> LDAPERM >= max(1,N). 94*> \endverbatim 95*> 96*> \param[in] PIV 97*> \verbatim 98*> PIV is INTEGER array, dimension (N) 99*> PIV is such that the nonzero entries are 100*> P( PIV( K ), K ) = 1. 101*> \endverbatim 102*> 103*> \param[out] RWORK 104*> \verbatim 105*> RWORK is REAL array, dimension (N) 106*> \endverbatim 107*> 108*> \param[out] RESID 109*> \verbatim 110*> RESID is REAL 111*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 112*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 113*> \endverbatim 114*> 115*> \param[in] RANK 116*> \verbatim 117*> RANK is INTEGER 118*> number of nonzero singular values of A. 119*> \endverbatim 120* 121* Authors: 122* ======== 123* 124*> \author Univ. of Tennessee 125*> \author Univ. of California Berkeley 126*> \author Univ. of Colorado Denver 127*> \author NAG Ltd. 128* 129*> \ingroup single_lin 130* 131* ===================================================================== 132 SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 133 $ PIV, RWORK, RESID, RANK ) 134* 135* -- LAPACK test routine -- 136* -- LAPACK is a software package provided by Univ. of Tennessee, -- 137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 138* 139* .. Scalar Arguments .. 140 REAL RESID 141 INTEGER LDA, LDAFAC, LDPERM, N, RANK 142 CHARACTER UPLO 143* .. 144* .. Array Arguments .. 145 REAL A( LDA, * ), AFAC( LDAFAC, * ), 146 $ PERM( LDPERM, * ), RWORK( * ) 147 INTEGER PIV( * ) 148* .. 149* 150* ===================================================================== 151* 152* .. Parameters .. 153 REAL ZERO, ONE 154 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 155* .. 156* .. Local Scalars .. 157 REAL ANORM, EPS, T 158 INTEGER I, J, K 159* .. 160* .. External Functions .. 161 REAL SDOT, SLAMCH, SLANSY 162 LOGICAL LSAME 163 EXTERNAL SDOT, SLAMCH, SLANSY, LSAME 164* .. 165* .. External Subroutines .. 166 EXTERNAL SSCAL, SSYR, STRMV 167* .. 168* .. Intrinsic Functions .. 169 INTRINSIC REAL 170* .. 171* .. Executable Statements .. 172* 173* Quick exit if N = 0. 174* 175 IF( N.LE.0 ) THEN 176 RESID = ZERO 177 RETURN 178 END IF 179* 180* Exit with RESID = 1/EPS if ANORM = 0. 181* 182 EPS = SLAMCH( 'Epsilon' ) 183 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) 184 IF( ANORM.LE.ZERO ) THEN 185 RESID = ONE / EPS 186 RETURN 187 END IF 188* 189* Compute the product U'*U, overwriting U. 190* 191 IF( LSAME( UPLO, 'U' ) ) THEN 192* 193 IF( RANK.LT.N ) THEN 194 DO 110 J = RANK + 1, N 195 DO 100 I = RANK + 1, J 196 AFAC( I, J ) = ZERO 197 100 CONTINUE 198 110 CONTINUE 199 END IF 200* 201 DO 120 K = N, 1, -1 202* 203* Compute the (K,K) element of the result. 204* 205 T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 206 AFAC( K, K ) = T 207* 208* Compute the rest of column K. 209* 210 CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 211 $ LDAFAC, AFAC( 1, K ), 1 ) 212* 213 120 CONTINUE 214* 215* Compute the product L*L', overwriting L. 216* 217 ELSE 218* 219 IF( RANK.LT.N ) THEN 220 DO 140 J = RANK + 1, N 221 DO 130 I = J, N 222 AFAC( I, J ) = ZERO 223 130 CONTINUE 224 140 CONTINUE 225 END IF 226* 227 DO 150 K = N, 1, -1 228* Add a multiple of column K of the factor L to each of 229* columns K+1 through N. 230* 231 IF( K+1.LE.N ) 232 $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 233 $ AFAC( K+1, K+1 ), LDAFAC ) 234* 235* Scale column K by the diagonal element. 236* 237 T = AFAC( K, K ) 238 CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 ) 239 150 CONTINUE 240* 241 END IF 242* 243* Form P*L*L'*P' or P*U'*U*P' 244* 245 IF( LSAME( UPLO, 'U' ) ) THEN 246* 247 DO 170 J = 1, N 248 DO 160 I = 1, N 249 IF( PIV( I ).LE.PIV( J ) ) THEN 250 IF( I.LE.J ) THEN 251 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 252 ELSE 253 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 254 END IF 255 END IF 256 160 CONTINUE 257 170 CONTINUE 258* 259* 260 ELSE 261* 262 DO 190 J = 1, N 263 DO 180 I = 1, N 264 IF( PIV( I ).GE.PIV( J ) ) THEN 265 IF( I.GE.J ) THEN 266 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 267 ELSE 268 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 269 END IF 270 END IF 271 180 CONTINUE 272 190 CONTINUE 273* 274 END IF 275* 276* Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A). 277* 278 IF( LSAME( UPLO, 'U' ) ) THEN 279 DO 210 J = 1, N 280 DO 200 I = 1, J 281 PERM( I, J ) = PERM( I, J ) - A( I, J ) 282 200 CONTINUE 283 210 CONTINUE 284 ELSE 285 DO 230 J = 1, N 286 DO 220 I = J, N 287 PERM( I, J ) = PERM( I, J ) - A( I, J ) 288 220 CONTINUE 289 230 CONTINUE 290 END IF 291* 292* Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or 293* ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ). 294* 295 RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK ) 296* 297 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 298* 299 RETURN 300* 301* End of SPST01 302* 303 END 304