1*> \brief \b CPPT01 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 CPPT01( 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 RWORK( * ) 20* COMPLEX A( * ), AFAC( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> CPPT01 reconstructs a Hermitian positive definite packed matrix A 30*> from its L*L' or U'*U factorization and computes the residual 31*> norm( L*L' - A ) / ( N * norm(A) * EPS ) or 32*> norm( U'*U - A ) / ( N * norm(A) * EPS ), 33*> where EPS is the machine epsilon, L' is the conjugate transpose of 34*> L, and U' is the conjugate transpose of U. 35*> \endverbatim 36* 37* Arguments: 38* ========== 39* 40*> \param[in] UPLO 41*> \verbatim 42*> UPLO is CHARACTER*1 43*> Specifies whether the upper or lower triangular part of the 44*> Hermitian matrix A is stored: 45*> = 'U': Upper triangular 46*> = 'L': Lower triangular 47*> \endverbatim 48*> 49*> \param[in] N 50*> \verbatim 51*> N is INTEGER 52*> The number of rows and columns of the matrix A. N >= 0. 53*> \endverbatim 54*> 55*> \param[in] A 56*> \verbatim 57*> A is COMPLEX array, dimension (N*(N+1)/2) 58*> The original Hermitian matrix A, stored as a packed 59*> triangular matrix. 60*> \endverbatim 61*> 62*> \param[in,out] AFAC 63*> \verbatim 64*> AFAC is COMPLEX array, dimension (N*(N+1)/2) 65*> On entry, the factor L or U from the L*L' or U'*U 66*> factorization of A, stored as a packed triangular matrix. 67*> Overwritten with the reconstructed matrix, and then with the 68*> difference L*L' - A (or U'*U - A). 69*> \endverbatim 70*> 71*> \param[out] RWORK 72*> \verbatim 73*> RWORK is REAL array, dimension (N) 74*> \endverbatim 75*> 76*> \param[out] RESID 77*> \verbatim 78*> RESID is REAL 79*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 80*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 81*> \endverbatim 82* 83* Authors: 84* ======== 85* 86*> \author Univ. of Tennessee 87*> \author Univ. of California Berkeley 88*> \author Univ. of Colorado Denver 89*> \author NAG Ltd. 90* 91*> \ingroup complex_lin 92* 93* ===================================================================== 94 SUBROUTINE CPPT01( UPLO, N, A, AFAC, RWORK, RESID ) 95* 96* -- LAPACK test routine -- 97* -- LAPACK is a software package provided by Univ. of Tennessee, -- 98* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 99* 100* .. Scalar Arguments .. 101 CHARACTER UPLO 102 INTEGER N 103 REAL RESID 104* .. 105* .. Array Arguments .. 106 REAL RWORK( * ) 107 COMPLEX A( * ), AFAC( * ) 108* .. 109* 110* ===================================================================== 111* 112* .. Parameters .. 113 REAL ZERO, ONE 114 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 115* .. 116* .. Local Scalars .. 117 INTEGER I, K, KC 118 REAL ANORM, EPS, TR 119 COMPLEX TC 120* .. 121* .. External Functions .. 122 LOGICAL LSAME 123 REAL CLANHP, SLAMCH 124 COMPLEX CDOTC 125 EXTERNAL LSAME, CLANHP, SLAMCH, CDOTC 126* .. 127* .. External Subroutines .. 128 EXTERNAL CHPR, CSCAL, CTPMV 129* .. 130* .. Intrinsic Functions .. 131 INTRINSIC AIMAG, REAL 132* .. 133* .. Executable Statements .. 134* 135* Quick exit if N = 0 136* 137 IF( N.LE.0 ) THEN 138 RESID = ZERO 139 RETURN 140 END IF 141* 142* Exit with RESID = 1/EPS if ANORM = 0. 143* 144 EPS = SLAMCH( 'Epsilon' ) 145 ANORM = CLANHP( '1', UPLO, N, A, RWORK ) 146 IF( ANORM.LE.ZERO ) THEN 147 RESID = ONE / EPS 148 RETURN 149 END IF 150* 151* Check the imaginary parts of the diagonal elements and return with 152* an error code if any are nonzero. 153* 154 KC = 1 155 IF( LSAME( UPLO, 'U' ) ) THEN 156 DO 10 K = 1, N 157 IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN 158 RESID = ONE / EPS 159 RETURN 160 END IF 161 KC = KC + K + 1 162 10 CONTINUE 163 ELSE 164 DO 20 K = 1, N 165 IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN 166 RESID = ONE / EPS 167 RETURN 168 END IF 169 KC = KC + N - K + 1 170 20 CONTINUE 171 END IF 172* 173* Compute the product U'*U, overwriting U. 174* 175 IF( LSAME( UPLO, 'U' ) ) THEN 176 KC = ( N*( N-1 ) ) / 2 + 1 177 DO 30 K = N, 1, -1 178* 179* Compute the (K,K) element of the result. 180* 181 TR = CDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 ) 182 AFAC( KC+K-1 ) = TR 183* 184* Compute the rest of column K. 185* 186 IF( K.GT.1 ) THEN 187 CALL CTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC, 188 $ AFAC( KC ), 1 ) 189 KC = KC - ( K-1 ) 190 END IF 191 30 CONTINUE 192* 193* Compute the difference L*L' - A 194* 195 KC = 1 196 DO 50 K = 1, N 197 DO 40 I = 1, K - 1 198 AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 ) 199 40 CONTINUE 200 AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - REAL( A( KC+K-1 ) ) 201 KC = KC + K 202 50 CONTINUE 203* 204* Compute the product L*L', overwriting L. 205* 206 ELSE 207 KC = ( N*( N+1 ) ) / 2 208 DO 60 K = N, 1, -1 209* 210* Add a multiple of column K of the factor L to each of 211* columns K+1 through N. 212* 213 IF( K.LT.N ) 214 $ CALL CHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1, 215 $ AFAC( KC+N-K+1 ) ) 216* 217* Scale column K by the diagonal element. 218* 219 TC = AFAC( KC ) 220 CALL CSCAL( N-K+1, TC, AFAC( KC ), 1 ) 221* 222 KC = KC - ( N-K+2 ) 223 60 CONTINUE 224* 225* Compute the difference U'*U - A 226* 227 KC = 1 228 DO 80 K = 1, N 229 AFAC( KC ) = AFAC( KC ) - REAL( A( KC ) ) 230 DO 70 I = K + 1, N 231 AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K ) 232 70 CONTINUE 233 KC = KC + N - K + 1 234 80 CONTINUE 235 END IF 236* 237* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 238* 239 RESID = CLANHP( '1', UPLO, N, AFAC, RWORK ) 240* 241 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 242* 243 RETURN 244* 245* End of CPPT01 246* 247 END 248