LAPACK 3.3.1
Linear Algebra PACKage

ctrt01.f

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