1*> \brief \b CSPT01 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 CSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER UPLO 15* INTEGER LDC, N 16* REAL RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* REAL RWORK( * ) 21* COMPLEX A( * ), AFAC( * ), C( LDC, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> CSPT01 reconstructs a symmetric indefinite packed matrix A from its 31*> diagonal pivoting factorization A = U*D*U' or A = L*D*L' and computes 32*> the residual 33*> norm( C - A ) / ( N * norm(A) * EPS ), 34*> where C is the reconstructed matrix and EPS is the machine epsilon. 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 order 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 symmetric matrix A, stored as a packed 59*> triangular matrix. 60*> \endverbatim 61*> 62*> \param[in] AFAC 63*> \verbatim 64*> AFAC is COMPLEX array, dimension (N*(N+1)/2) 65*> The factored form of the matrix A, stored as a packed 66*> triangular matrix. AFAC contains the block diagonal matrix D 67*> and the multipliers used to obtain the factor L or U from the 68*> L*D*L' or U*D*U' factorization as computed by CSPTRF. 69*> \endverbatim 70*> 71*> \param[in] IPIV 72*> \verbatim 73*> IPIV is INTEGER array, dimension (N) 74*> The pivot indices from CSPTRF. 75*> \endverbatim 76*> 77*> \param[out] C 78*> \verbatim 79*> C is COMPLEX array, dimension (LDC,N) 80*> \endverbatim 81*> 82*> \param[in] LDC 83*> \verbatim 84*> LDC is INTEGER 85*> The leading dimension of the array C. LDC >= max(1,N). 86*> \endverbatim 87*> 88*> \param[out] RWORK 89*> \verbatim 90*> RWORK is REAL array, dimension (N) 91*> \endverbatim 92*> 93*> \param[out] RESID 94*> \verbatim 95*> RESID is REAL 96*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 97*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 98*> \endverbatim 99* 100* Authors: 101* ======== 102* 103*> \author Univ. of Tennessee 104*> \author Univ. of California Berkeley 105*> \author Univ. of Colorado Denver 106*> \author NAG Ltd. 107* 108*> \ingroup complex_lin 109* 110* ===================================================================== 111 SUBROUTINE CSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID ) 112* 113* -- LAPACK test routine -- 114* -- LAPACK is a software package provided by Univ. of Tennessee, -- 115* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 116* 117* .. Scalar Arguments .. 118 CHARACTER UPLO 119 INTEGER LDC, N 120 REAL RESID 121* .. 122* .. Array Arguments .. 123 INTEGER IPIV( * ) 124 REAL RWORK( * ) 125 COMPLEX A( * ), AFAC( * ), C( LDC, * ) 126* .. 127* 128* ===================================================================== 129* 130* .. Parameters .. 131 REAL ZERO, ONE 132 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 133 COMPLEX CZERO, CONE 134 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 135 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 136* .. 137* .. Local Scalars .. 138 INTEGER I, INFO, J, JC 139 REAL ANORM, EPS 140* .. 141* .. External Functions .. 142 LOGICAL LSAME 143 REAL CLANSP, CLANSY, SLAMCH 144 EXTERNAL LSAME, CLANSP, CLANSY, SLAMCH 145* .. 146* .. External Subroutines .. 147 EXTERNAL CLAVSP, CLASET 148* .. 149* .. Intrinsic Functions .. 150 INTRINSIC REAL 151* .. 152* .. Executable Statements .. 153* 154* Quick exit if N = 0. 155* 156 IF( N.LE.0 ) THEN 157 RESID = ZERO 158 RETURN 159 END IF 160* 161* Determine EPS and the norm of A. 162* 163 EPS = SLAMCH( 'Epsilon' ) 164 ANORM = CLANSP( '1', UPLO, N, A, RWORK ) 165* 166* Initialize C to the identity matrix. 167* 168 CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC ) 169* 170* Call CLAVSP to form the product D * U' (or D * L' ). 171* 172 CALL CLAVSP( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, IPIV, C, 173 $ LDC, INFO ) 174* 175* Call CLAVSP again to multiply by U ( or L ). 176* 177 CALL CLAVSP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C, 178 $ LDC, INFO ) 179* 180* Compute the difference C - A . 181* 182 IF( LSAME( UPLO, 'U' ) ) THEN 183 JC = 0 184 DO 20 J = 1, N 185 DO 10 I = 1, J 186 C( I, J ) = C( I, J ) - A( JC+I ) 187 10 CONTINUE 188 JC = JC + J 189 20 CONTINUE 190 ELSE 191 JC = 1 192 DO 40 J = 1, N 193 DO 30 I = J, N 194 C( I, J ) = C( I, J ) - A( JC+I-J ) 195 30 CONTINUE 196 JC = JC + N - J + 1 197 40 CONTINUE 198 END IF 199* 200* Compute norm( C - A ) / ( N * norm(A) * EPS ) 201* 202 RESID = CLANSY( '1', UPLO, N, C, LDC, RWORK ) 203* 204 IF( ANORM.LE.ZERO ) THEN 205 IF( RESID.NE.ZERO ) 206 $ RESID = ONE / EPS 207 ELSE 208 RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS 209 END IF 210* 211 RETURN 212* 213* End of CSPT01 214* 215 END 216