![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGET01 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDAFAC, M, N 00016 * DOUBLE PRECISION RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * DOUBLE PRECISION RWORK( * ) 00021 * COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> ZGET01 reconstructs a matrix A from its L*U factorization and 00031 *> computes the residual 00032 *> norm(L*U - A) / ( N * norm(A) * EPS ), 00033 *> where EPS is the machine epsilon. 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] M 00040 *> \verbatim 00041 *> M is INTEGER 00042 *> The number of rows of the matrix A. M >= 0. 00043 *> \endverbatim 00044 *> 00045 *> \param[in] N 00046 *> \verbatim 00047 *> N is INTEGER 00048 *> The number of columns of the matrix A. N >= 0. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] A 00052 *> \verbatim 00053 *> A is COMPLEX*16 array, dimension (LDA,N) 00054 *> The original M x N matrix A. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] LDA 00058 *> \verbatim 00059 *> LDA is INTEGER 00060 *> The leading dimension of the array A. LDA >= max(1,M). 00061 *> \endverbatim 00062 *> 00063 *> \param[in,out] AFAC 00064 *> \verbatim 00065 *> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) 00066 *> The factored form of the matrix A. AFAC contains the factors 00067 *> L and U from the L*U factorization as computed by ZGETRF. 00068 *> Overwritten with the reconstructed matrix, and then with the 00069 *> difference L*U - A. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] LDAFAC 00073 *> \verbatim 00074 *> LDAFAC is INTEGER 00075 *> The leading dimension of the array AFAC. LDAFAC >= max(1,M). 00076 *> \endverbatim 00077 *> 00078 *> \param[in] IPIV 00079 *> \verbatim 00080 *> IPIV is INTEGER array, dimension (N) 00081 *> The pivot indices from ZGETRF. 00082 *> \endverbatim 00083 *> 00084 *> \param[out] RWORK 00085 *> \verbatim 00086 *> RWORK is DOUBLE PRECISION array, dimension (M) 00087 *> \endverbatim 00088 *> 00089 *> \param[out] RESID 00090 *> \verbatim 00091 *> RESID is DOUBLE PRECISION 00092 *> norm(L*U - A) / ( N * norm(A) * EPS ) 00093 *> \endverbatim 00094 * 00095 * Authors: 00096 * ======== 00097 * 00098 *> \author Univ. of Tennessee 00099 *> \author Univ. of California Berkeley 00100 *> \author Univ. of Colorado Denver 00101 *> \author NAG Ltd. 00102 * 00103 *> \date November 2011 00104 * 00105 *> \ingroup complex16_lin 00106 * 00107 * ===================================================================== 00108 SUBROUTINE ZGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK, 00109 $ RESID ) 00110 * 00111 * -- LAPACK test routine (version 3.4.0) -- 00112 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00113 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00114 * November 2011 00115 * 00116 * .. Scalar Arguments .. 00117 INTEGER LDA, LDAFAC, M, N 00118 DOUBLE PRECISION RESID 00119 * .. 00120 * .. Array Arguments .. 00121 INTEGER IPIV( * ) 00122 DOUBLE PRECISION RWORK( * ) 00123 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ) 00124 * .. 00125 * 00126 * ===================================================================== 00127 * 00128 * .. Parameters .. 00129 DOUBLE PRECISION ZERO, ONE 00130 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00131 COMPLEX*16 CONE 00132 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00133 * .. 00134 * .. Local Scalars .. 00135 INTEGER I, J, K 00136 DOUBLE PRECISION ANORM, EPS 00137 COMPLEX*16 T 00138 * .. 00139 * .. External Functions .. 00140 DOUBLE PRECISION DLAMCH, ZLANGE 00141 COMPLEX*16 ZDOTU 00142 EXTERNAL DLAMCH, ZLANGE, ZDOTU 00143 * .. 00144 * .. External Subroutines .. 00145 EXTERNAL ZGEMV, ZLASWP, ZSCAL, ZTRMV 00146 * .. 00147 * .. Intrinsic Functions .. 00148 INTRINSIC DBLE, MIN 00149 * .. 00150 * .. Executable Statements .. 00151 * 00152 * Quick exit if M = 0 or N = 0. 00153 * 00154 IF( M.LE.0 .OR. N.LE.0 ) THEN 00155 RESID = ZERO 00156 RETURN 00157 END IF 00158 * 00159 * Determine EPS and the norm of A. 00160 * 00161 EPS = DLAMCH( 'Epsilon' ) 00162 ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) 00163 * 00164 * Compute the product L*U and overwrite AFAC with the result. 00165 * A column at a time of the product is obtained, starting with 00166 * column N. 00167 * 00168 DO 10 K = N, 1, -1 00169 IF( K.GT.M ) THEN 00170 CALL ZTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC, 00171 $ LDAFAC, AFAC( 1, K ), 1 ) 00172 ELSE 00173 * 00174 * Compute elements (K+1:M,K) 00175 * 00176 T = AFAC( K, K ) 00177 IF( K+1.LE.M ) THEN 00178 CALL ZSCAL( M-K, T, AFAC( K+1, K ), 1 ) 00179 CALL ZGEMV( 'No transpose', M-K, K-1, CONE, 00180 $ AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1, 00181 $ CONE, AFAC( K+1, K ), 1 ) 00182 END IF 00183 * 00184 * Compute the (K,K) element 00185 * 00186 AFAC( K, K ) = T + ZDOTU( K-1, AFAC( K, 1 ), LDAFAC, 00187 $ AFAC( 1, K ), 1 ) 00188 * 00189 * Compute elements (1:K-1,K) 00190 * 00191 CALL ZTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC, 00192 $ LDAFAC, AFAC( 1, K ), 1 ) 00193 END IF 00194 10 CONTINUE 00195 CALL ZLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 ) 00196 * 00197 * Compute the difference L*U - A and store in AFAC. 00198 * 00199 DO 30 J = 1, N 00200 DO 20 I = 1, M 00201 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 00202 20 CONTINUE 00203 30 CONTINUE 00204 * 00205 * Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 00206 * 00207 RESID = ZLANGE( '1', M, N, AFAC, LDAFAC, RWORK ) 00208 * 00209 IF( ANORM.LE.ZERO ) THEN 00210 IF( RESID.NE.ZERO ) 00211 $ RESID = ONE / EPS 00212 ELSE 00213 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 00214 END IF 00215 * 00216 RETURN 00217 * 00218 * End of ZGET01 00219 * 00220 END