1*> \brief \b DTRT01 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DTRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND, 12* WORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, UPLO 16* INTEGER LDA, LDAINV, N 17* DOUBLE PRECISION RCOND, RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> DTRT01 computes the residual for a triangular matrix A times its 30*> inverse: 31*> RESID = norm( A*AINV - I ) / ( N * norm(A) * norm(AINV) * EPS ), 32*> where EPS is the machine epsilon. 33*> \endverbatim 34* 35* Arguments: 36* ========== 37* 38*> \param[in] UPLO 39*> \verbatim 40*> UPLO is CHARACTER*1 41*> Specifies whether the matrix A is upper or lower triangular. 42*> = 'U': Upper triangular 43*> = 'L': Lower triangular 44*> \endverbatim 45*> 46*> \param[in] DIAG 47*> \verbatim 48*> DIAG is CHARACTER*1 49*> Specifies whether or not the matrix A is unit triangular. 50*> = 'N': Non-unit triangular 51*> = 'U': Unit triangular 52*> \endverbatim 53*> 54*> \param[in] N 55*> \verbatim 56*> N is INTEGER 57*> The order of the matrix A. N >= 0. 58*> \endverbatim 59*> 60*> \param[in] A 61*> \verbatim 62*> A is DOUBLE PRECISION array, dimension (LDA,N) 63*> The triangular matrix A. If UPLO = 'U', the leading n by n 64*> upper triangular part of the array A contains the upper 65*> triangular matrix, and the strictly lower triangular part of 66*> A is not referenced. If UPLO = 'L', the leading n by n lower 67*> triangular part of the array A contains the lower triangular 68*> matrix, and the strictly upper triangular part of A is not 69*> referenced. If DIAG = 'U', the diagonal elements of A are 70*> also not referenced and are assumed to be 1. 71*> \endverbatim 72*> 73*> \param[in] LDA 74*> \verbatim 75*> LDA is INTEGER 76*> The leading dimension of the array A. LDA >= max(1,N). 77*> \endverbatim 78*> 79*> \param[in,out] AINV 80*> \verbatim 81*> AINV is DOUBLE PRECISION array, dimension (LDAINV,N) 82*> On entry, the (triangular) inverse of the matrix A, in the 83*> same storage format as A. 84*> On exit, the contents of AINV are destroyed. 85*> \endverbatim 86*> 87*> \param[in] LDAINV 88*> \verbatim 89*> LDAINV is INTEGER 90*> The leading dimension of the array AINV. LDAINV >= max(1,N). 91*> \endverbatim 92*> 93*> \param[out] RCOND 94*> \verbatim 95*> RCOND is DOUBLE PRECISION 96*> The reciprocal condition number of A, computed as 97*> 1/(norm(A) * norm(AINV)). 98*> \endverbatim 99*> 100*> \param[out] WORK 101*> \verbatim 102*> WORK is DOUBLE PRECISION array, dimension (N) 103*> \endverbatim 104*> 105*> \param[out] RESID 106*> \verbatim 107*> RESID is DOUBLE PRECISION 108*> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ) 109*> \endverbatim 110* 111* Authors: 112* ======== 113* 114*> \author Univ. of Tennessee 115*> \author Univ. of California Berkeley 116*> \author Univ. of Colorado Denver 117*> \author NAG Ltd. 118* 119*> \date December 2016 120* 121*> \ingroup double_lin 122* 123* ===================================================================== 124 SUBROUTINE DTRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND, 125 $ WORK, RESID ) 126* 127* -- LAPACK test routine (version 3.7.0) -- 128* -- LAPACK is a software package provided by Univ. of Tennessee, -- 129* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 130* December 2016 131* 132* .. Scalar Arguments .. 133 CHARACTER DIAG, UPLO 134 INTEGER LDA, LDAINV, N 135 DOUBLE PRECISION RCOND, RESID 136* .. 137* .. Array Arguments .. 138 DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), WORK( * ) 139* .. 140* 141* ===================================================================== 142* 143* .. Parameters .. 144 DOUBLE PRECISION ZERO, ONE 145 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 146* .. 147* .. Local Scalars .. 148 INTEGER J 149 DOUBLE PRECISION AINVNM, ANORM, EPS 150* .. 151* .. External Functions .. 152 LOGICAL LSAME 153 DOUBLE PRECISION DLAMCH, DLANTR 154 EXTERNAL LSAME, DLAMCH, DLANTR 155* .. 156* .. External Subroutines .. 157 EXTERNAL DTRMV 158* .. 159* .. Intrinsic Functions .. 160 INTRINSIC DBLE 161* .. 162* .. Executable Statements .. 163* 164* Quick exit if N = 0 165* 166 IF( N.LE.0 ) THEN 167 RCOND = ONE 168 RESID = ZERO 169 RETURN 170 END IF 171* 172* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 173* 174 EPS = DLAMCH( 'Epsilon' ) 175 ANORM = DLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK ) 176 AINVNM = DLANTR( '1', UPLO, DIAG, N, N, AINV, LDAINV, WORK ) 177 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 178 RCOND = ZERO 179 RESID = ONE / EPS 180 RETURN 181 END IF 182 RCOND = ( ONE / ANORM ) / AINVNM 183* 184* Set the diagonal of AINV to 1 if AINV has unit diagonal. 185* 186 IF( LSAME( DIAG, 'U' ) ) THEN 187 DO 10 J = 1, N 188 AINV( J, J ) = ONE 189 10 CONTINUE 190 END IF 191* 192* Compute A * AINV, overwriting AINV. 193* 194 IF( LSAME( UPLO, 'U' ) ) THEN 195 DO 20 J = 1, N 196 CALL DTRMV( 'Upper', 'No transpose', DIAG, J, A, LDA, 197 $ AINV( 1, J ), 1 ) 198 20 CONTINUE 199 ELSE 200 DO 30 J = 1, N 201 CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J+1, A( J, J ), 202 $ LDA, AINV( J, J ), 1 ) 203 30 CONTINUE 204 END IF 205* 206* Subtract 1 from each diagonal element to form A*AINV - I. 207* 208 DO 40 J = 1, N 209 AINV( J, J ) = AINV( J, J ) - ONE 210 40 CONTINUE 211* 212* Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS) 213* 214 RESID = DLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, WORK ) 215* 216 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS 217* 218 RETURN 219* 220* End of DTRT01 221* 222 END 223