1 SUBROUTINE CTRT01( UPLO, DIAG, N, A, LDA, AINV, LDAINV, RCOND,
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, UPLO
10 INTEGER LDA, LDAINV, N
11 REAL RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( LDA, * ), AINV( LDAINV, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CTRT01 computes the residual for a triangular matrix A times its
22 * inverse:
23 * RESID = norm( A*AINV - I ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the matrix A is upper or lower triangular.
31 * = 'U': Upper triangular
32 * = 'L': Lower triangular
33 *
34 * DIAG (input) CHARACTER*1
35 * Specifies whether or not the matrix A is unit triangular.
36 * = 'N': Non-unit triangular
37 * = 'U': Unit triangular
38 *
39 * N (input) INTEGER
40 * The order of the matrix A. N >= 0.
41 *
42 * A (input) COMPLEX array, dimension (LDA,N)
43 * The triangular matrix A. If UPLO = 'U', the leading n by n
44 * upper triangular part of the array A contains the upper
45 * triangular matrix, and the strictly lower triangular part of
46 * A is not referenced. If UPLO = 'L', the leading n by n lower
47 * triangular part of the array A contains the lower triangular
48 * matrix, and the strictly upper triangular part of A is not
49 * referenced. If DIAG = 'U', the diagonal elements of A are
50 * also not referenced and are assumed to be 1.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the array A. LDA >= max(1,N).
54 *
55 * AINV (input) COMPLEX array, dimension (LDAINV,N)
56 * On entry, the (triangular) inverse of the matrix A, in the
57 * same storage format as A.
58 * On exit, the contents of AINV are destroyed.
59 *
60 * LDAINV (input) INTEGER
61 * The leading dimension of the array AINV. LDAINV >= max(1,N).
62 *
63 * RCOND (output) REAL
64 * The reciprocal condition number of A, computed as
65 * 1/(norm(A) * norm(AINV)).
66 *
67 * RWORK (workspace) REAL array, dimension (N)
68 *
69 * RESID (output) REAL
70 * norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ZERO, ONE
76 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER J
80 REAL AINVNM, ANORM, EPS
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 REAL CLANTR, SLAMCH
85 EXTERNAL LSAME, CLANTR, SLAMCH
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL CTRMV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC REAL
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0
96 *
97 IF( N.LE.0 ) THEN
98 RCOND = ONE
99 RESID = ZERO
100 RETURN
101 END IF
102 *
103 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
104 *
105 EPS = SLAMCH( 'Epsilon' )
106 ANORM = CLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK )
107 AINVNM = CLANTR( '1', UPLO, DIAG, N, N, AINV, LDAINV, RWORK )
108 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
109 RCOND = ZERO
110 RESID = ONE / EPS
111 RETURN
112 END IF
113 RCOND = ( ONE / ANORM ) / AINVNM
114 *
115 * Set the diagonal of AINV to 1 if AINV has unit diagonal.
116 *
117 IF( LSAME( DIAG, 'U' ) ) THEN
118 DO 10 J = 1, N
119 AINV( J, J ) = ONE
120 10 CONTINUE
121 END IF
122 *
123 * Compute A * AINV, overwriting AINV.
124 *
125 IF( LSAME( UPLO, 'U' ) ) THEN
126 DO 20 J = 1, N
127 CALL CTRMV( 'Upper', 'No transpose', DIAG, J, A, LDA,
128 $ AINV( 1, J ), 1 )
129 20 CONTINUE
130 ELSE
131 DO 30 J = 1, N
132 CALL CTRMV( 'Lower', 'No transpose', DIAG, N-J+1, A( J, J ),
133 $ LDA, AINV( J, J ), 1 )
134 30 CONTINUE
135 END IF
136 *
137 * Subtract 1 from each diagonal element to form A*AINV - I.
138 *
139 DO 40 J = 1, N
140 AINV( J, J ) = AINV( J, J ) - ONE
141 40 CONTINUE
142 *
143 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
144 *
145 RESID = CLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, RWORK )
146 *
147 RESID = ( ( RESID*RCOND ) / REAL( N ) ) / EPS
148 *
149 RETURN
150 *
151 * End of CTRT01
152 *
153 END
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, UPLO
10 INTEGER LDA, LDAINV, N
11 REAL RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( LDA, * ), AINV( LDAINV, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CTRT01 computes the residual for a triangular matrix A times its
22 * inverse:
23 * RESID = norm( A*AINV - I ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the matrix A is upper or lower triangular.
31 * = 'U': Upper triangular
32 * = 'L': Lower triangular
33 *
34 * DIAG (input) CHARACTER*1
35 * Specifies whether or not the matrix A is unit triangular.
36 * = 'N': Non-unit triangular
37 * = 'U': Unit triangular
38 *
39 * N (input) INTEGER
40 * The order of the matrix A. N >= 0.
41 *
42 * A (input) COMPLEX array, dimension (LDA,N)
43 * The triangular matrix A. If UPLO = 'U', the leading n by n
44 * upper triangular part of the array A contains the upper
45 * triangular matrix, and the strictly lower triangular part of
46 * A is not referenced. If UPLO = 'L', the leading n by n lower
47 * triangular part of the array A contains the lower triangular
48 * matrix, and the strictly upper triangular part of A is not
49 * referenced. If DIAG = 'U', the diagonal elements of A are
50 * also not referenced and are assumed to be 1.
51 *
52 * LDA (input) INTEGER
53 * The leading dimension of the array A. LDA >= max(1,N).
54 *
55 * AINV (input) COMPLEX array, dimension (LDAINV,N)
56 * On entry, the (triangular) inverse of the matrix A, in the
57 * same storage format as A.
58 * On exit, the contents of AINV are destroyed.
59 *
60 * LDAINV (input) INTEGER
61 * The leading dimension of the array AINV. LDAINV >= max(1,N).
62 *
63 * RCOND (output) REAL
64 * The reciprocal condition number of A, computed as
65 * 1/(norm(A) * norm(AINV)).
66 *
67 * RWORK (workspace) REAL array, dimension (N)
68 *
69 * RESID (output) REAL
70 * norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ZERO, ONE
76 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER J
80 REAL AINVNM, ANORM, EPS
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 REAL CLANTR, SLAMCH
85 EXTERNAL LSAME, CLANTR, SLAMCH
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL CTRMV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC REAL
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0
96 *
97 IF( N.LE.0 ) THEN
98 RCOND = ONE
99 RESID = ZERO
100 RETURN
101 END IF
102 *
103 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
104 *
105 EPS = SLAMCH( 'Epsilon' )
106 ANORM = CLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK )
107 AINVNM = CLANTR( '1', UPLO, DIAG, N, N, AINV, LDAINV, RWORK )
108 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
109 RCOND = ZERO
110 RESID = ONE / EPS
111 RETURN
112 END IF
113 RCOND = ( ONE / ANORM ) / AINVNM
114 *
115 * Set the diagonal of AINV to 1 if AINV has unit diagonal.
116 *
117 IF( LSAME( DIAG, 'U' ) ) THEN
118 DO 10 J = 1, N
119 AINV( J, J ) = ONE
120 10 CONTINUE
121 END IF
122 *
123 * Compute A * AINV, overwriting AINV.
124 *
125 IF( LSAME( UPLO, 'U' ) ) THEN
126 DO 20 J = 1, N
127 CALL CTRMV( 'Upper', 'No transpose', DIAG, J, A, LDA,
128 $ AINV( 1, J ), 1 )
129 20 CONTINUE
130 ELSE
131 DO 30 J = 1, N
132 CALL CTRMV( 'Lower', 'No transpose', DIAG, N-J+1, A( J, J ),
133 $ LDA, AINV( J, J ), 1 )
134 30 CONTINUE
135 END IF
136 *
137 * Subtract 1 from each diagonal element to form A*AINV - I.
138 *
139 DO 40 J = 1, N
140 AINV( J, J ) = AINV( J, J ) - ONE
141 40 CONTINUE
142 *
143 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
144 *
145 RESID = CLANTR( '1', UPLO, 'Non-unit', N, N, AINV, LDAINV, RWORK )
146 *
147 RESID = ( ( RESID*RCOND ) / REAL( N ) ) / EPS
148 *
149 RETURN
150 *
151 * End of CTRT01
152 *
153 END