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*> \date November 2011 92* 93*> \ingroup complex_lin 94* 95* ===================================================================== 96 SUBROUTINE CPPT01( UPLO, N, A, AFAC, RWORK, RESID ) 97* 98* -- LAPACK test routine (version 3.4.0) -- 99* -- LAPACK is a software package provided by Univ. of Tennessee, -- 100* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 101* November 2011 102* 103* .. Scalar Arguments .. 104 CHARACTER UPLO 105 INTEGER N 106 REAL RESID 107* .. 108* .. Array Arguments .. 109 REAL RWORK( * ) 110 COMPLEX A( * ), AFAC( * ) 111* .. 112* 113* ===================================================================== 114* 115* .. Parameters .. 116 REAL ZERO, ONE 117 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 118* .. 119* .. Local Scalars .. 120 INTEGER I, K, KC 121 REAL ANORM, EPS, TR 122 COMPLEX TC 123* .. 124* .. External Functions .. 125 LOGICAL LSAME 126 REAL CLANHP, SLAMCH 127 COMPLEX CDOTC 128 EXTERNAL LSAME, CLANHP, SLAMCH, CDOTC 129* .. 130* .. External Subroutines .. 131 EXTERNAL CHPR, CSCAL, CTPMV 132* .. 133* .. Intrinsic Functions .. 134 INTRINSIC AIMAG, REAL 135* .. 136* .. Executable Statements .. 137* 138* Quick exit if N = 0 139* 140 IF( N.LE.0 ) THEN 141 RESID = ZERO 142 RETURN 143 END IF 144* 145* Exit with RESID = 1/EPS if ANORM = 0. 146* 147 EPS = SLAMCH( 'Epsilon' ) 148 ANORM = CLANHP( '1', UPLO, N, A, RWORK ) 149 IF( ANORM.LE.ZERO ) THEN 150 RESID = ONE / EPS 151 RETURN 152 END IF 153* 154* Check the imaginary parts of the diagonal elements and return with 155* an error code if any are nonzero. 156* 157 KC = 1 158 IF( LSAME( UPLO, 'U' ) ) THEN 159 DO 10 K = 1, N 160 IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN 161 RESID = ONE / EPS 162 RETURN 163 END IF 164 KC = KC + K + 1 165 10 CONTINUE 166 ELSE 167 DO 20 K = 1, N 168 IF( AIMAG( AFAC( KC ) ).NE.ZERO ) THEN 169 RESID = ONE / EPS 170 RETURN 171 END IF 172 KC = KC + N - K + 1 173 20 CONTINUE 174 END IF 175* 176* Compute the product U'*U, overwriting U. 177* 178 IF( LSAME( UPLO, 'U' ) ) THEN 179 KC = ( N*( N-1 ) ) / 2 + 1 180 DO 30 K = N, 1, -1 181* 182* Compute the (K,K) element of the result. 183* 184 TR = CDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 ) 185 AFAC( KC+K-1 ) = TR 186* 187* Compute the rest of column K. 188* 189 IF( K.GT.1 ) THEN 190 CALL CTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC, 191 $ AFAC( KC ), 1 ) 192 KC = KC - ( K-1 ) 193 END IF 194 30 CONTINUE 195* 196* Compute the difference L*L' - A 197* 198 KC = 1 199 DO 50 K = 1, N 200 DO 40 I = 1, K - 1 201 AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 ) 202 40 CONTINUE 203 AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - REAL( A( KC+K-1 ) ) 204 KC = KC + K 205 50 CONTINUE 206* 207* Compute the product L*L', overwriting L. 208* 209 ELSE 210 KC = ( N*( N+1 ) ) / 2 211 DO 60 K = N, 1, -1 212* 213* Add a multiple of column K of the factor L to each of 214* columns K+1 through N. 215* 216 IF( K.LT.N ) 217 $ CALL CHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1, 218 $ AFAC( KC+N-K+1 ) ) 219* 220* Scale column K by the diagonal element. 221* 222 TC = AFAC( KC ) 223 CALL CSCAL( N-K+1, TC, AFAC( KC ), 1 ) 224* 225 KC = KC - ( N-K+2 ) 226 60 CONTINUE 227* 228* Compute the difference U'*U - A 229* 230 KC = 1 231 DO 80 K = 1, N 232 AFAC( KC ) = AFAC( KC ) - REAL( A( KC ) ) 233 DO 70 I = K + 1, N 234 AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K ) 235 70 CONTINUE 236 KC = KC + N - K + 1 237 80 CONTINUE 238 END IF 239* 240* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 241* 242 RESID = CLANHP( '1', UPLO, N, AFAC, RWORK ) 243* 244 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 245* 246 RETURN 247* 248* End of CPPT01 249* 250 END 251