![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b STRT01 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 STRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND, 00012 * WORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER DIAG, UPLO 00016 * INTEGER LDA, LDAINV, N 00017 * REAL RCOND, RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * REAL A( LDA, * ), AINV( LDAINV, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> STRT01 computes the residual for a triangular matrix A times its 00030 *> inverse: 00031 *> RESID = norm( A*AINV - I ) / ( N * norm(A) * norm(AINV) * EPS ), 00032 *> where EPS is the machine epsilon. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] UPLO 00039 *> \verbatim 00040 *> UPLO is CHARACTER*1 00041 *> Specifies whether the matrix A is upper or lower triangular. 00042 *> = 'U': Upper triangular 00043 *> = 'L': Lower triangular 00044 *> \endverbatim 00045 *> 00046 *> \param[in] DIAG 00047 *> \verbatim 00048 *> DIAG is CHARACTER*1 00049 *> Specifies whether or not the matrix A is unit triangular. 00050 *> = 'N': Non-unit triangular 00051 *> = 'U': Unit triangular 00052 *> \endverbatim 00053 *> 00054 *> \param[in] N 00055 *> \verbatim 00056 *> N is INTEGER 00057 *> The order of the matrix A. N >= 0. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] A 00061 *> \verbatim 00062 *> A is REAL array, dimension (LDA,N) 00063 *> The triangular matrix A. If UPLO = 'U', the leading n by n 00064 *> upper triangular part of the array A contains the upper 00065 *> triangular matrix, and the strictly lower triangular part of 00066 *> A is not referenced. If UPLO = 'L', the leading n by n lower 00067 *> triangular part of the array A contains the lower triangular 00068 *> matrix, and the strictly upper triangular part of A is not 00069 *> referenced. If DIAG = 'U', the diagonal elements of A are 00070 *> also not referenced and are assumed to be 1. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDA 00074 *> \verbatim 00075 *> LDA is INTEGER 00076 *> The leading dimension of the array A. LDA >= max(1,N). 00077 *> \endverbatim 00078 *> 00079 *> \param[in,out] AINV 00080 *> \verbatim 00081 *> AINV is REAL array, dimension (LDAINV,N) 00082 *> On entry, the (triangular) inverse of the matrix A, in the 00083 *> same storage format as A. 00084 *> On exit, the contents of AINV are destroyed. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] LDAINV 00088 *> \verbatim 00089 *> LDAINV is INTEGER 00090 *> The leading dimension of the array AINV. LDAINV >= max(1,N). 00091 *> \endverbatim 00092 *> 00093 *> \param[out] RCOND 00094 *> \verbatim 00095 *> RCOND is REAL 00096 *> The reciprocal condition number of A, computed as 00097 *> 1/(norm(A) * norm(AINV)). 00098 *> \endverbatim 00099 *> 00100 *> \param[out] WORK 00101 *> \verbatim 00102 *> WORK is REAL array, dimension (N) 00103 *> \endverbatim 00104 *> 00105 *> \param[out] RESID 00106 *> \verbatim 00107 *> RESID is REAL 00108 *> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ) 00109 *> \endverbatim 00110 * 00111 * Authors: 00112 * ======== 00113 * 00114 *> \author Univ. of Tennessee 00115 *> \author Univ. of California Berkeley 00116 *> \author Univ. of Colorado Denver 00117 *> \author NAG Ltd. 00118 * 00119 *> \date November 2011 00120 * 00121 *> \ingroup single_lin 00122 * 00123 * ===================================================================== 00124 SUBROUTINE STRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND, 00125 $ WORK, RESID ) 00126 * 00127 * -- LAPACK test routine (version 3.4.0) -- 00128 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00130 * November 2011 00131 * 00132 * .. Scalar Arguments .. 00133 CHARACTER DIAG, UPLO 00134 INTEGER LDA, LDAINV, N 00135 REAL RCOND, RESID 00136 * .. 00137 * .. Array Arguments .. 00138 REAL A( LDA, * ), AINV( LDAINV, * ), WORK( * ) 00139 * .. 00140 * 00141 * ===================================================================== 00142 * 00143 * .. Parameters .. 00144 REAL ZERO, ONE 00145 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00146 * .. 00147 * .. Local Scalars .. 00148 INTEGER J 00149 REAL AINVNM, ANORM, EPS 00150 * .. 00151 * .. External Functions .. 00152 LOGICAL LSAME 00153 REAL SLAMCH, SLANTR 00154 EXTERNAL LSAME, SLAMCH, SLANTR 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL STRMV 00158 * .. 00159 * .. Intrinsic Functions .. 00160 INTRINSIC REAL 00161 * .. 00162 * .. Executable Statements .. 00163 * 00164 * Quick exit if N = 0 00165 * 00166 IF( N.LE.0 ) THEN 00167 RCOND = ONE 00168 RESID = ZERO 00169 RETURN 00170 END IF 00171 * 00172 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00173 * 00174 EPS = SLAMCH( 'Epsilon' ) 00175 ANORM = SLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK ) 00176 AINVNM = SLANTR( '1', UPLO, DIAG, N, N, AINV, LDAINV, WORK ) 00177 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00178 RCOND = ZERO 00179 RESID = ONE / EPS 00180 RETURN 00181 END IF 00182 RCOND = ( ONE / ANORM ) / AINVNM 00183 * 00184 * Set the diagonal of AINV to 1 if AINV has unit diagonal. 00185 * 00186 IF( LSAME( DIAG, 'U' ) ) THEN 00187 DO 10 J = 1, N 00188 AINV( J, J ) = ONE 00189 10 CONTINUE 00190 END IF 00191 * 00192 * Compute A * AINV, overwriting AINV. 00193 * 00194 IF( LSAME( UPLO, 'U' ) ) THEN 00195 DO 20 J = 1, N 00196 CALL STRMV( 'Upper', 'No transpose', DIAG, J, A, LDA, 00197 $ AINV( 1, J ), 1 ) 00198 20 CONTINUE 00199 ELSE 00200 DO 30 J = 1, N 00201 CALL STRMV( 'Lower', 'No transpose', DIAG, N-J+1, A( J, J ), 00202 $ LDA, AINV( J, J ), 1 ) 00203 30 CONTINUE 00204 END IF 00205 * 00206 * Subtract 1 from each diagonal element to form A*AINV - I. 00207 * 00208 DO 40 J = 1, N 00209 AINV( J, J ) = AINV( J, J ) - ONE 00210 40 CONTINUE 00211 * 00212 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS) 00213 * 00214 RESID = SLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, WORK ) 00215 * 00216 RESID = ( ( RESID*RCOND ) / REAL( N ) ) / EPS 00217 * 00218 RETURN 00219 * 00220 * End of STRT01 00221 * 00222 END