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*> \date December 2016 130* 131*> \ingroup single_lin 132* 133* ===================================================================== 134 SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM, 135 $ PIV, RWORK, RESID, RANK ) 136* 137* -- LAPACK test routine (version 3.7.0) -- 138* -- LAPACK is a software package provided by Univ. of Tennessee, -- 139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 140* December 2016 141* 142* .. Scalar Arguments .. 143 REAL RESID 144 INTEGER LDA, LDAFAC, LDPERM, N, RANK 145 CHARACTER UPLO 146* .. 147* .. Array Arguments .. 148 REAL A( LDA, * ), AFAC( LDAFAC, * ), 149 $ PERM( LDPERM, * ), RWORK( * ) 150 INTEGER PIV( * ) 151* .. 152* 153* ===================================================================== 154* 155* .. Parameters .. 156 REAL ZERO, ONE 157 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 158* .. 159* .. Local Scalars .. 160 REAL ANORM, EPS, T 161 INTEGER I, J, K 162* .. 163* .. External Functions .. 164 REAL SDOT, SLAMCH, SLANSY 165 LOGICAL LSAME 166 EXTERNAL SDOT, SLAMCH, SLANSY, LSAME 167* .. 168* .. External Subroutines .. 169 EXTERNAL SSCAL, SSYR, STRMV 170* .. 171* .. Intrinsic Functions .. 172 INTRINSIC REAL 173* .. 174* .. Executable Statements .. 175* 176* Quick exit if N = 0. 177* 178 IF( N.LE.0 ) THEN 179 RESID = ZERO 180 RETURN 181 END IF 182* 183* Exit with RESID = 1/EPS if ANORM = 0. 184* 185 EPS = SLAMCH( 'Epsilon' ) 186 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) 187 IF( ANORM.LE.ZERO ) THEN 188 RESID = ONE / EPS 189 RETURN 190 END IF 191* 192* Compute the product U'*U, overwriting U. 193* 194 IF( LSAME( UPLO, 'U' ) ) THEN 195* 196 IF( RANK.LT.N ) THEN 197 DO 110 J = RANK + 1, N 198 DO 100 I = RANK + 1, J 199 AFAC( I, J ) = ZERO 200 100 CONTINUE 201 110 CONTINUE 202 END IF 203* 204 DO 120 K = N, 1, -1 205* 206* Compute the (K,K) element of the result. 207* 208 T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 ) 209 AFAC( K, K ) = T 210* 211* Compute the rest of column K. 212* 213 CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC, 214 $ LDAFAC, AFAC( 1, K ), 1 ) 215* 216 120 CONTINUE 217* 218* Compute the product L*L', overwriting L. 219* 220 ELSE 221* 222 IF( RANK.LT.N ) THEN 223 DO 140 J = RANK + 1, N 224 DO 130 I = J, N 225 AFAC( I, J ) = ZERO 226 130 CONTINUE 227 140 CONTINUE 228 END IF 229* 230 DO 150 K = N, 1, -1 231* Add a multiple of column K of the factor L to each of 232* columns K+1 through N. 233* 234 IF( K+1.LE.N ) 235 $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1, 236 $ AFAC( K+1, K+1 ), LDAFAC ) 237* 238* Scale column K by the diagonal element. 239* 240 T = AFAC( K, K ) 241 CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 ) 242 150 CONTINUE 243* 244 END IF 245* 246* Form P*L*L'*P' or P*U'*U*P' 247* 248 IF( LSAME( UPLO, 'U' ) ) THEN 249* 250 DO 170 J = 1, N 251 DO 160 I = 1, N 252 IF( PIV( I ).LE.PIV( J ) ) THEN 253 IF( I.LE.J ) THEN 254 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 255 ELSE 256 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 257 END IF 258 END IF 259 160 CONTINUE 260 170 CONTINUE 261* 262* 263 ELSE 264* 265 DO 190 J = 1, N 266 DO 180 I = 1, N 267 IF( PIV( I ).GE.PIV( J ) ) THEN 268 IF( I.GE.J ) THEN 269 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J ) 270 ELSE 271 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I ) 272 END IF 273 END IF 274 180 CONTINUE 275 190 CONTINUE 276* 277 END IF 278* 279* Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A). 280* 281 IF( LSAME( UPLO, 'U' ) ) THEN 282 DO 210 J = 1, N 283 DO 200 I = 1, J 284 PERM( I, J ) = PERM( I, J ) - A( I, J ) 285 200 CONTINUE 286 210 CONTINUE 287 ELSE 288 DO 230 J = 1, N 289 DO 220 I = J, N 290 PERM( I, J ) = PERM( I, J ) - A( I, J ) 291 220 CONTINUE 292 230 CONTINUE 293 END IF 294* 295* Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or 296* ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ). 297* 298 RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK ) 299* 300 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 301* 302 RETURN 303* 304* End of SPST01 305* 306 END 307