![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGET01 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 CGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDAFAC, M, N 00016 * REAL RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * REAL RWORK( * ) 00021 * COMPLEX A( LDA, * ), AFAC( LDAFAC, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CGET01 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 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 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 CGETRF. 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 CGETRF. 00082 *> \endverbatim 00083 *> 00084 *> \param[out] RWORK 00085 *> \verbatim 00086 *> RWORK is REAL array, dimension (M) 00087 *> \endverbatim 00088 *> 00089 *> \param[out] RESID 00090 *> \verbatim 00091 *> RESID is REAL 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 complex_lin 00106 * 00107 * ===================================================================== 00108 SUBROUTINE CGET01( 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 REAL RESID 00119 * .. 00120 * .. Array Arguments .. 00121 INTEGER IPIV( * ) 00122 REAL RWORK( * ) 00123 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ) 00124 * .. 00125 * 00126 * ===================================================================== 00127 * 00128 * .. Parameters .. 00129 REAL ONE, ZERO 00130 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00131 COMPLEX CONE 00132 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00133 * .. 00134 * .. Local Scalars .. 00135 INTEGER I, J, K 00136 REAL ANORM, EPS 00137 COMPLEX T 00138 * .. 00139 * .. External Functions .. 00140 REAL CLANGE, SLAMCH 00141 COMPLEX CDOTU 00142 EXTERNAL CLANGE, SLAMCH, CDOTU 00143 * .. 00144 * .. External Subroutines .. 00145 EXTERNAL CGEMV, CLASWP, CSCAL, CTRMV 00146 * .. 00147 * .. Intrinsic Functions .. 00148 INTRINSIC MIN, REAL 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 = SLAMCH( 'Epsilon' ) 00162 ANORM = CLANGE( '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 CTRMV( '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 CSCAL( M-K, T, AFAC( K+1, K ), 1 ) 00179 CALL CGEMV( '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 + CDOTU( K-1, AFAC( K, 1 ), LDAFAC, 00187 $ AFAC( 1, K ), 1 ) 00188 * 00189 * Compute elements (1:K-1,K) 00190 * 00191 CALL CTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC, 00192 $ LDAFAC, AFAC( 1, K ), 1 ) 00193 END IF 00194 10 CONTINUE 00195 CALL CLASWP( 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 = CLANGE( '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/REAL( N ) )/ANORM ) / EPS 00214 END IF 00215 * 00216 RETURN 00217 * 00218 * End of CGET01 00219 * 00220 END