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