1*> \brief \b ZGET01 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 ZGET01( 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 RWORK( * ) 21* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> ZGET01 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*16 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*16 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 ZGETRF. 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 ZGETRF. 82*> \endverbatim 83*> 84*> \param[out] RWORK 85*> \verbatim 86*> RWORK is DOUBLE PRECISION array, dimension (M) 87*> \endverbatim 88*> 89*> \param[out] RESID 90*> \verbatim 91*> RESID is DOUBLE PRECISION 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*> \date November 2011 104* 105*> \ingroup complex16_lin 106* 107* ===================================================================== 108 SUBROUTINE ZGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK, 109 $ RESID ) 110* 111* -- LAPACK test routine (version 3.4.0) -- 112* -- LAPACK is a software package provided by Univ. of Tennessee, -- 113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 114* November 2011 115* 116* .. Scalar Arguments .. 117 INTEGER LDA, LDAFAC, M, N 118 DOUBLE PRECISION RESID 119* .. 120* .. Array Arguments .. 121 INTEGER IPIV( * ) 122 DOUBLE PRECISION RWORK( * ) 123 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 124* .. 125* 126* ===================================================================== 127* 128* .. Parameters .. 129 DOUBLE PRECISION ZERO, ONE 130 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 131 COMPLEX*16 CONE 132 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 133* .. 134* .. Local Scalars .. 135 INTEGER I, J, K 136 DOUBLE PRECISION ANORM, EPS 137 COMPLEX*16 T 138* .. 139* .. External Functions .. 140 DOUBLE PRECISION DLAMCH, ZLANGE 141 COMPLEX*16 ZDOTU 142 EXTERNAL DLAMCH, ZLANGE, ZDOTU 143* .. 144* .. External Subroutines .. 145 EXTERNAL ZGEMV, ZLASWP, ZSCAL, ZTRMV 146* .. 147* .. Intrinsic Functions .. 148 INTRINSIC DBLE, MIN 149* .. 150* .. Executable Statements .. 151* 152* Quick exit if M = 0 or N = 0. 153* 154 IF( M.LE.0 .OR. N.LE.0 ) THEN 155 RESID = ZERO 156 RETURN 157 END IF 158* 159* Determine EPS and the norm of A. 160* 161 EPS = DLAMCH( 'Epsilon' ) 162 ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) 163* 164* Compute the product L*U and overwrite AFAC with the result. 165* A column at a time of the product is obtained, starting with 166* column N. 167* 168 DO 10 K = N, 1, -1 169 IF( K.GT.M ) THEN 170 CALL ZTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC, 171 $ LDAFAC, AFAC( 1, K ), 1 ) 172 ELSE 173* 174* Compute elements (K+1:M,K) 175* 176 T = AFAC( K, K ) 177 IF( K+1.LE.M ) THEN 178 CALL ZSCAL( M-K, T, AFAC( K+1, K ), 1 ) 179 CALL ZGEMV( 'No transpose', M-K, K-1, CONE, 180 $ AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1, 181 $ CONE, AFAC( K+1, K ), 1 ) 182 END IF 183* 184* Compute the (K,K) element 185* 186 AFAC( K, K ) = T + ZDOTU( K-1, AFAC( K, 1 ), LDAFAC, 187 $ AFAC( 1, K ), 1 ) 188* 189* Compute elements (1:K-1,K) 190* 191 CALL ZTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC, 192 $ LDAFAC, AFAC( 1, K ), 1 ) 193 END IF 194 10 CONTINUE 195 CALL ZLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 ) 196* 197* Compute the difference L*U - A and store in AFAC. 198* 199 DO 30 J = 1, N 200 DO 20 I = 1, M 201 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 202 20 CONTINUE 203 30 CONTINUE 204* 205* Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 206* 207 RESID = ZLANGE( '1', M, N, AFAC, LDAFAC, RWORK ) 208* 209 IF( ANORM.LE.ZERO ) THEN 210 IF( RESID.NE.ZERO ) 211 $ RESID = ONE / EPS 212 ELSE 213 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 214 END IF 215* 216 RETURN 217* 218* End of ZGET01 219* 220 END 221