![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGET03 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 SGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, 00012 * RCOND, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDAINV, LDWORK, N 00016 * REAL RCOND, RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 00020 * $ WORK( LDWORK, * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SGET03 computes the residual for a general matrix times its inverse: 00030 *> norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ), 00031 *> where EPS is the machine epsilon. 00032 *> \endverbatim 00033 * 00034 * Arguments: 00035 * ========== 00036 * 00037 *> \param[in] N 00038 *> \verbatim 00039 *> N is INTEGER 00040 *> The number of rows and columns of the matrix A. N >= 0. 00041 *> \endverbatim 00042 *> 00043 *> \param[in] A 00044 *> \verbatim 00045 *> A is REAL array, dimension (LDA,N) 00046 *> The original N x N matrix A. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] LDA 00050 *> \verbatim 00051 *> LDA is INTEGER 00052 *> The leading dimension of the array A. LDA >= max(1,N). 00053 *> \endverbatim 00054 *> 00055 *> \param[in] AINV 00056 *> \verbatim 00057 *> AINV is REAL array, dimension (LDAINV,N) 00058 *> The inverse of the matrix A. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] LDAINV 00062 *> \verbatim 00063 *> LDAINV is INTEGER 00064 *> The leading dimension of the array AINV. LDAINV >= max(1,N). 00065 *> \endverbatim 00066 *> 00067 *> \param[out] WORK 00068 *> \verbatim 00069 *> WORK is REAL array, dimension (LDWORK,N) 00070 *> \endverbatim 00071 *> 00072 *> \param[in] LDWORK 00073 *> \verbatim 00074 *> LDWORK is INTEGER 00075 *> The leading dimension of the array WORK. LDWORK >= max(1,N). 00076 *> \endverbatim 00077 *> 00078 *> \param[out] RWORK 00079 *> \verbatim 00080 *> RWORK is REAL array, dimension (N) 00081 *> \endverbatim 00082 *> 00083 *> \param[out] RCOND 00084 *> \verbatim 00085 *> RCOND is REAL 00086 *> The reciprocal of the condition number of A, computed as 00087 *> ( 1/norm(A) ) / norm(AINV). 00088 *> \endverbatim 00089 *> 00090 *> \param[out] RESID 00091 *> \verbatim 00092 *> RESID is REAL 00093 *> norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS ) 00094 *> \endverbatim 00095 * 00096 * Authors: 00097 * ======== 00098 * 00099 *> \author Univ. of Tennessee 00100 *> \author Univ. of California Berkeley 00101 *> \author Univ. of Colorado Denver 00102 *> \author NAG Ltd. 00103 * 00104 *> \date November 2011 00105 * 00106 *> \ingroup single_lin 00107 * 00108 * ===================================================================== 00109 SUBROUTINE SGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, 00110 $ RCOND, RESID ) 00111 * 00112 * -- LAPACK test routine (version 3.4.0) -- 00113 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00114 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00115 * November 2011 00116 * 00117 * .. Scalar Arguments .. 00118 INTEGER LDA, LDAINV, LDWORK, N 00119 REAL RCOND, RESID 00120 * .. 00121 * .. Array Arguments .. 00122 REAL A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 00123 $ WORK( LDWORK, * ) 00124 * .. 00125 * 00126 * ===================================================================== 00127 * 00128 * .. Parameters .. 00129 REAL ZERO, ONE 00130 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00131 * .. 00132 * .. Local Scalars .. 00133 INTEGER I 00134 REAL AINVNM, ANORM, EPS 00135 * .. 00136 * .. External Functions .. 00137 REAL SLAMCH, SLANGE 00138 EXTERNAL SLAMCH, SLANGE 00139 * .. 00140 * .. External Subroutines .. 00141 EXTERNAL SGEMM 00142 * .. 00143 * .. Intrinsic Functions .. 00144 INTRINSIC REAL 00145 * .. 00146 * .. Executable Statements .. 00147 * 00148 * Quick exit if N = 0. 00149 * 00150 IF( N.LE.0 ) THEN 00151 RCOND = ONE 00152 RESID = ZERO 00153 RETURN 00154 END IF 00155 * 00156 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00157 * 00158 EPS = SLAMCH( 'Epsilon' ) 00159 ANORM = SLANGE( '1', N, N, A, LDA, RWORK ) 00160 AINVNM = SLANGE( '1', N, N, AINV, LDAINV, RWORK ) 00161 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00162 RCOND = ZERO 00163 RESID = ONE / EPS 00164 RETURN 00165 END IF 00166 RCOND = ( ONE / ANORM ) / AINVNM 00167 * 00168 * Compute I - A * AINV 00169 * 00170 CALL SGEMM( 'No transpose', 'No transpose', N, N, N, -ONE, 00171 $ AINV, LDAINV, A, LDA, ZERO, WORK, LDWORK ) 00172 DO 10 I = 1, N 00173 WORK( I, I ) = ONE + WORK( I, I ) 00174 10 CONTINUE 00175 * 00176 * Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS) 00177 * 00178 RESID = SLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00179 * 00180 RESID = ( ( RESID*RCOND ) / EPS ) / REAL( N ) 00181 * 00182 RETURN 00183 * 00184 * End of SGET03 00185 * 00186 END