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