1*> \brief \b DPBT01 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 DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK, 12* RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER KD, LDA, LDAFAC, N 17* DOUBLE PRECISION RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> DPBT01 reconstructs a symmetric positive definite band matrix A from 30*> 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*> symmetric 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] KD 56*> \verbatim 57*> KD is INTEGER 58*> The number of super-diagonals of the matrix A if UPLO = 'U', 59*> or the number of sub-diagonals if UPLO = 'L'. KD >= 0. 60*> \endverbatim 61*> 62*> \param[in] A 63*> \verbatim 64*> A is DOUBLE PRECISION array, dimension (LDA,N) 65*> The original symmetric band matrix A. If UPLO = 'U', the 66*> upper triangular part of A is stored as a band matrix; if 67*> UPLO = 'L', the lower triangular part of A is stored. The 68*> columns of the appropriate triangle are stored in the columns 69*> of A and the diagonals of the triangle are stored in the rows 70*> of A. See DPBTRF for further details. 71*> \endverbatim 72*> 73*> \param[in] LDA 74*> \verbatim 75*> LDA is INTEGER. 76*> The leading dimension of the array A. LDA >= max(1,KD+1). 77*> \endverbatim 78*> 79*> \param[in] AFAC 80*> \verbatim 81*> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N) 82*> The factored form of the matrix A. AFAC contains the factor 83*> L or U from the L*L' or U'*U factorization in band storage 84*> format, as computed by DPBTRF. 85*> \endverbatim 86*> 87*> \param[in] LDAFAC 88*> \verbatim 89*> LDAFAC is INTEGER 90*> The leading dimension of the array AFAC. 91*> LDAFAC >= max(1,KD+1). 92*> \endverbatim 93*> 94*> \param[out] RWORK 95*> \verbatim 96*> RWORK is DOUBLE PRECISION array, dimension (N) 97*> \endverbatim 98*> 99*> \param[out] RESID 100*> \verbatim 101*> RESID is DOUBLE PRECISION 102*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 103*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 104*> \endverbatim 105* 106* Authors: 107* ======== 108* 109*> \author Univ. of Tennessee 110*> \author Univ. of California Berkeley 111*> \author Univ. of Colorado Denver 112*> \author NAG Ltd. 113* 114*> \date November 2011 115* 116*> \ingroup double_lin 117* 118* ===================================================================== 119 SUBROUTINE DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK, 120 $ RESID ) 121* 122* -- LAPACK test routine (version 3.4.0) -- 123* -- LAPACK is a software package provided by Univ. of Tennessee, -- 124* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 125* November 2011 126* 127* .. Scalar Arguments .. 128 CHARACTER UPLO 129 INTEGER KD, LDA, LDAFAC, N 130 DOUBLE PRECISION RESID 131* .. 132* .. Array Arguments .. 133 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * ) 134* .. 135* 136* ===================================================================== 137* 138* 139* .. Parameters .. 140 DOUBLE PRECISION ZERO, ONE 141 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 142* .. 143* .. Local Scalars .. 144 INTEGER I, J, K, KC, KLEN, ML, MU 145 DOUBLE PRECISION ANORM, EPS, T 146* .. 147* .. External Functions .. 148 LOGICAL LSAME 149 DOUBLE PRECISION DDOT, DLAMCH, DLANSB 150 EXTERNAL LSAME, DDOT, DLAMCH, DLANSB 151* .. 152* .. External Subroutines .. 153 EXTERNAL DSCAL, DSYR, DTRMV 154* .. 155* .. Intrinsic Functions .. 156 INTRINSIC DBLE, MAX, MIN 157* .. 158* .. Executable Statements .. 159* 160* Quick exit if N = 0. 161* 162 IF( N.LE.0 ) THEN 163 RESID = ZERO 164 RETURN 165 END IF 166* 167* Exit with RESID = 1/EPS if ANORM = 0. 168* 169 EPS = DLAMCH( 'Epsilon' ) 170 ANORM = DLANSB( '1', UPLO, N, KD, A, LDA, RWORK ) 171 IF( ANORM.LE.ZERO ) THEN 172 RESID = ONE / EPS 173 RETURN 174 END IF 175* 176* Compute the product U'*U, overwriting U. 177* 178 IF( LSAME( UPLO, 'U' ) ) THEN 179 DO 10 K = N, 1, -1 180 KC = MAX( 1, KD+2-K ) 181 KLEN = KD + 1 - KC 182* 183* Compute the (K,K) element of the result. 184* 185 T = DDOT( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 ) 186 AFAC( KD+1, K ) = T 187* 188* Compute the rest of column K. 189* 190 IF( KLEN.GT.0 ) 191 $ CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', KLEN, 192 $ AFAC( KD+1, K-KLEN ), LDAFAC-1, 193 $ AFAC( KC, K ), 1 ) 194* 195 10 CONTINUE 196* 197* UPLO = 'L': Compute the product L*L', overwriting L. 198* 199 ELSE 200 DO 20 K = N, 1, -1 201 KLEN = MIN( KD, N-K ) 202* 203* Add a multiple of column K of the factor L to each of 204* columns K+1 through N. 205* 206 IF( KLEN.GT.0 ) 207 $ CALL DSYR( 'Lower', KLEN, ONE, AFAC( 2, K ), 1, 208 $ AFAC( 1, K+1 ), LDAFAC-1 ) 209* 210* Scale column K by the diagonal element. 211* 212 T = AFAC( 1, K ) 213 CALL DSCAL( KLEN+1, T, AFAC( 1, K ), 1 ) 214* 215 20 CONTINUE 216 END IF 217* 218* Compute the difference L*L' - A or U'*U - A. 219* 220 IF( LSAME( UPLO, 'U' ) ) THEN 221 DO 40 J = 1, N 222 MU = MAX( 1, KD+2-J ) 223 DO 30 I = MU, KD + 1 224 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 225 30 CONTINUE 226 40 CONTINUE 227 ELSE 228 DO 60 J = 1, N 229 ML = MIN( KD+1, N-J+1 ) 230 DO 50 I = 1, ML 231 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 232 50 CONTINUE 233 60 CONTINUE 234 END IF 235* 236* Compute norm( L*L' - A ) / ( N * norm(A) * EPS ) 237* 238 RESID = DLANSB( 'I', UPLO, N, KD, AFAC, LDAFAC, RWORK ) 239* 240 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 241* 242 RETURN 243* 244* End of DPBT01 245* 246 END 247