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