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