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