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