LAPACK 3.3.0

# dtrt01.f

Go to the documentation of this file.
00001       SUBROUTINE DTRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND,
00002      \$                   WORK, RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          DIAG, UPLO
00010       INTEGER            LDA, LDAINV, N
00011       DOUBLE PRECISION   RCOND, RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DTRT01 computes the residual for a triangular matrix A times its
00021 *  inverse:
00022 *     RESID = norm( A*AINV - I ) / ( N * norm(A) * norm(AINV) * EPS ),
00023 *  where EPS is the machine epsilon.
00024 *
00025 *  Arguments
00026 *  ==========
00027 *
00028 *  UPLO    (input) CHARACTER*1
00029 *          Specifies whether the matrix A is upper or lower triangular.
00030 *          = 'U':  Upper triangular
00031 *          = 'L':  Lower triangular
00032 *
00033 *  DIAG    (input) CHARACTER*1
00034 *          Specifies whether or not the matrix A is unit triangular.
00035 *          = 'N':  Non-unit triangular
00036 *          = 'U':  Unit triangular
00037 *
00038 *  N       (input) INTEGER
00039 *          The order of the matrix A.  N >= 0.
00040 *
00041 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00042 *          The triangular matrix A.  If UPLO = 'U', the leading n by n
00043 *          upper triangular part of the array A contains the upper
00044 *          triangular matrix, and the strictly lower triangular part of
00045 *          A is not referenced.  If UPLO = 'L', the leading n by n lower
00046 *          triangular part of the array A contains the lower triangular
00047 *          matrix, and the strictly upper triangular part of A is not
00048 *          referenced.  If DIAG = 'U', the diagonal elements of A are
00049 *          also not referenced and are assumed to be 1.
00050 *
00051 *  LDA     (input) INTEGER
00052 *          The leading dimension of the array A.  LDA >= max(1,N).
00053 *
00054 *  AINV    (input/output) DOUBLE PRECISION array, dimension (LDAINV,N)
00055 *          On entry, the (triangular) inverse of the matrix A, in the
00056 *          same storage format as A.
00057 *          On exit, the contents of AINV are destroyed.
00058 *
00059 *  LDAINV  (input) INTEGER
00060 *          The leading dimension of the array AINV.  LDAINV >= max(1,N).
00061 *
00062 *  RCOND   (output) DOUBLE PRECISION
00063 *          The reciprocal condition number of A, computed as
00064 *          1/(norm(A) * norm(AINV)).
00065 *
00066 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00067 *
00068 *  RESID   (output) DOUBLE PRECISION
00069 *          norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
00070 *
00071 *  =====================================================================
00072 *
00073 *     .. Parameters ..
00074       DOUBLE PRECISION   ZERO, ONE
00075       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00076 *     ..
00077 *     .. Local Scalars ..
00078       INTEGER            J
00079       DOUBLE PRECISION   AINVNM, ANORM, EPS
00080 *     ..
00081 *     .. External Functions ..
00082       LOGICAL            LSAME
00083       DOUBLE PRECISION   DLAMCH, DLANTR
00084       EXTERNAL           LSAME, DLAMCH, DLANTR
00085 *     ..
00086 *     .. External Subroutines ..
00087       EXTERNAL           DTRMV
00088 *     ..
00089 *     .. Intrinsic Functions ..
00090       INTRINSIC          DBLE
00091 *     ..
00092 *     .. Executable Statements ..
00093 *
00094 *     Quick exit if N = 0
00095 *
00096       IF( N.LE.0 ) THEN
00097          RCOND = ONE
00098          RESID = ZERO
00099          RETURN
00100       END IF
00101 *
00102 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00103 *
00104       EPS = DLAMCH( 'Epsilon' )
00105       ANORM = DLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK )
00106       AINVNM = DLANTR( '1', UPLO, DIAG, N, N, AINV, LDAINV, WORK )
00107       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00108          RCOND = ZERO
00109          RESID = ONE / EPS
00110          RETURN
00111       END IF
00112       RCOND = ( ONE / ANORM ) / AINVNM
00113 *
00114 *     Set the diagonal of AINV to 1 if AINV has unit diagonal.
00115 *
00116       IF( LSAME( DIAG, 'U' ) ) THEN
00117          DO 10 J = 1, N
00118             AINV( J, J ) = ONE
00119    10    CONTINUE
00120       END IF
00121 *
00122 *     Compute A * AINV, overwriting AINV.
00123 *
00124       IF( LSAME( UPLO, 'U' ) ) THEN
00125          DO 20 J = 1, N
00126             CALL DTRMV( 'Upper', 'No transpose', DIAG, J, A, LDA,
00127      \$                  AINV( 1, J ), 1 )
00128    20    CONTINUE
00129       ELSE
00130          DO 30 J = 1, N
00131             CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J+1, A( J, J ),
00132      \$                  LDA, AINV( J, J ), 1 )
00133    30    CONTINUE
00134       END IF
00135 *
00136 *     Subtract 1 from each diagonal element to form A*AINV - I.
00137 *
00138       DO 40 J = 1, N
00139          AINV( J, J ) = AINV( J, J ) - ONE
00140    40 CONTINUE
00141 *
00142 *     Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
00143 *
00144       RESID = DLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, WORK )
00145 *
00146       RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
00147 *
00148       RETURN
00149 *
00150 *     End of DTRT01
00151 *
00152       END