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