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