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