LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtrt01.f
Go to the documentation of this file.
00001 *> \brief \b DTRT01
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 DTRT01( 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 *       DOUBLE PRECISION   RCOND, RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DTRT01 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (N)
00103 *> \endverbatim
00104 *>
00105 *> \param[out] RESID
00106 *> \verbatim
00107 *>          RESID is DOUBLE PRECISION
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 double_lin
00122 *
00123 *  =====================================================================
00124       SUBROUTINE DTRT01( 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       DOUBLE PRECISION   RCOND, RESID
00136 *     ..
00137 *     .. Array Arguments ..
00138       DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), WORK( * )
00139 *     ..
00140 *
00141 *  =====================================================================
00142 *
00143 *     .. Parameters ..
00144       DOUBLE PRECISION   ZERO, ONE
00145       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00146 *     ..
00147 *     .. Local Scalars ..
00148       INTEGER            J
00149       DOUBLE PRECISION   AINVNM, ANORM, EPS
00150 *     ..
00151 *     .. External Functions ..
00152       LOGICAL            LSAME
00153       DOUBLE PRECISION   DLAMCH, DLANTR
00154       EXTERNAL           LSAME, DLAMCH, DLANTR
00155 *     ..
00156 *     .. External Subroutines ..
00157       EXTERNAL           DTRMV
00158 *     ..
00159 *     .. Intrinsic Functions ..
00160       INTRINSIC          DBLE
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 = DLAMCH( 'Epsilon' )
00175       ANORM = DLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK )
00176       AINVNM = DLANTR( '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 DTRMV( '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 DTRMV( '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 = DLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, WORK )
00215 *
00216       RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
00217 *
00218       RETURN
00219 *
00220 *     End of DTRT01
00221 *
00222       END
 All Files Functions