1 SUBROUTINE CTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B,
2 $ LDB, WORK, 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, TRANS, UPLO
10 INTEGER LDA, LDB, LDX, N, NRHS
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CTRT02 computes the residual for the computed solution to a
23 * triangular system of linear equations A*x = b, A**T *x = b,
24 * or A**H *x = b. Here A is a triangular matrix, A**T is the transpose
25 * of A, A**H is the conjugate transpose of A, and x and b are N by NRHS
26 * matrices. The test ratio is the maximum over the number of right
27 * hand sides of
28 * norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
29 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the matrix A is upper or lower triangular.
36 * = 'U': Upper triangular
37 * = 'L': Lower triangular
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the operation applied to A.
41 * = 'N': A *x = b (No transpose)
42 * = 'T': A**T *x = b (Transpose)
43 * = 'C': A**H *x = b (Conjugate transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * Specifies whether or not the matrix A is unit triangular.
47 * = 'N': Non-unit triangular
48 * = 'U': Unit triangular
49 *
50 * N (input) INTEGER
51 * The order of the matrix A. N >= 0.
52 *
53 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * A (input) COMPLEX array, dimension (LDA,N)
58 * The triangular matrix A. If UPLO = 'U', the leading n by n
59 * upper triangular part of the array A contains the upper
60 * triangular matrix, and the strictly lower triangular part of
61 * A is not referenced. If UPLO = 'L', the leading n by n lower
62 * triangular part of the array A contains the lower triangular
63 * matrix, and the strictly upper triangular part of A is not
64 * referenced. If DIAG = 'U', the diagonal elements of A are
65 * also not referenced and are assumed to be 1.
66 *
67 * LDA (input) INTEGER
68 * The leading dimension of the array A. LDA >= max(1,N).
69 *
70 * X (input) COMPLEX array, dimension (LDX,NRHS)
71 * The computed solution vectors for the system of linear
72 * equations.
73 *
74 * LDX (input) INTEGER
75 * The leading dimension of the array X. LDX >= max(1,N).
76 *
77 * B (input) COMPLEX array, dimension (LDB,NRHS)
78 * The right hand side vectors for the system of linear
79 * equations.
80 *
81 * LDB (input) INTEGER
82 * The leading dimension of the array B. LDB >= max(1,N).
83 *
84 * WORK (workspace) COMPLEX array, dimension (N)
85 *
86 * RWORK (workspace) REAL array, dimension (N)
87 *
88 * RESID (output) REAL
89 * The maximum over the number of right hand sides of
90 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 REAL ZERO, ONE
96 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
97 * ..
98 * .. Local Scalars ..
99 INTEGER J
100 REAL ANORM, BNORM, EPS, XNORM
101 * ..
102 * .. External Functions ..
103 LOGICAL LSAME
104 REAL CLANTR, SCASUM, SLAMCH
105 EXTERNAL LSAME, CLANTR, SCASUM, SLAMCH
106 * ..
107 * .. External Subroutines ..
108 EXTERNAL CAXPY, CCOPY, CTRMV
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC CMPLX, MAX
112 * ..
113 * .. Executable Statements ..
114 *
115 * Quick exit if N = 0 or NRHS = 0
116 *
117 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
118 RESID = ZERO
119 RETURN
120 END IF
121 *
122 * Compute the 1-norm of A or A**H.
123 *
124 IF( LSAME( TRANS, 'N' ) ) THEN
125 ANORM = CLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK )
126 ELSE
127 ANORM = CLANTR( 'I', UPLO, DIAG, N, N, A, LDA, RWORK )
128 END IF
129 *
130 * Exit with RESID = 1/EPS if ANORM = 0.
131 *
132 EPS = SLAMCH( 'Epsilon' )
133 IF( ANORM.LE.ZERO ) THEN
134 RESID = ONE / EPS
135 RETURN
136 END IF
137 *
138 * Compute the maximum over the number of right hand sides of
139 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS )
140 *
141 RESID = ZERO
142 DO 10 J = 1, NRHS
143 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
144 CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
145 CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
146 BNORM = SCASUM( N, WORK, 1 )
147 XNORM = SCASUM( N, X( 1, J ), 1 )
148 IF( XNORM.LE.ZERO ) THEN
149 RESID = ONE / EPS
150 ELSE
151 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
152 END IF
153 10 CONTINUE
154 *
155 RETURN
156 *
157 * End of CTRT02
158 *
159 END
2 $ LDB, WORK, 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, TRANS, UPLO
10 INTEGER LDA, LDB, LDX, N, NRHS
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CTRT02 computes the residual for the computed solution to a
23 * triangular system of linear equations A*x = b, A**T *x = b,
24 * or A**H *x = b. Here A is a triangular matrix, A**T is the transpose
25 * of A, A**H is the conjugate transpose of A, and x and b are N by NRHS
26 * matrices. The test ratio is the maximum over the number of right
27 * hand sides of
28 * norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
29 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the matrix A is upper or lower triangular.
36 * = 'U': Upper triangular
37 * = 'L': Lower triangular
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the operation applied to A.
41 * = 'N': A *x = b (No transpose)
42 * = 'T': A**T *x = b (Transpose)
43 * = 'C': A**H *x = b (Conjugate transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * Specifies whether or not the matrix A is unit triangular.
47 * = 'N': Non-unit triangular
48 * = 'U': Unit triangular
49 *
50 * N (input) INTEGER
51 * The order of the matrix A. N >= 0.
52 *
53 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * A (input) COMPLEX array, dimension (LDA,N)
58 * The triangular matrix A. If UPLO = 'U', the leading n by n
59 * upper triangular part of the array A contains the upper
60 * triangular matrix, and the strictly lower triangular part of
61 * A is not referenced. If UPLO = 'L', the leading n by n lower
62 * triangular part of the array A contains the lower triangular
63 * matrix, and the strictly upper triangular part of A is not
64 * referenced. If DIAG = 'U', the diagonal elements of A are
65 * also not referenced and are assumed to be 1.
66 *
67 * LDA (input) INTEGER
68 * The leading dimension of the array A. LDA >= max(1,N).
69 *
70 * X (input) COMPLEX array, dimension (LDX,NRHS)
71 * The computed solution vectors for the system of linear
72 * equations.
73 *
74 * LDX (input) INTEGER
75 * The leading dimension of the array X. LDX >= max(1,N).
76 *
77 * B (input) COMPLEX array, dimension (LDB,NRHS)
78 * The right hand side vectors for the system of linear
79 * equations.
80 *
81 * LDB (input) INTEGER
82 * The leading dimension of the array B. LDB >= max(1,N).
83 *
84 * WORK (workspace) COMPLEX array, dimension (N)
85 *
86 * RWORK (workspace) REAL array, dimension (N)
87 *
88 * RESID (output) REAL
89 * The maximum over the number of right hand sides of
90 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 REAL ZERO, ONE
96 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
97 * ..
98 * .. Local Scalars ..
99 INTEGER J
100 REAL ANORM, BNORM, EPS, XNORM
101 * ..
102 * .. External Functions ..
103 LOGICAL LSAME
104 REAL CLANTR, SCASUM, SLAMCH
105 EXTERNAL LSAME, CLANTR, SCASUM, SLAMCH
106 * ..
107 * .. External Subroutines ..
108 EXTERNAL CAXPY, CCOPY, CTRMV
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC CMPLX, MAX
112 * ..
113 * .. Executable Statements ..
114 *
115 * Quick exit if N = 0 or NRHS = 0
116 *
117 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
118 RESID = ZERO
119 RETURN
120 END IF
121 *
122 * Compute the 1-norm of A or A**H.
123 *
124 IF( LSAME( TRANS, 'N' ) ) THEN
125 ANORM = CLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK )
126 ELSE
127 ANORM = CLANTR( 'I', UPLO, DIAG, N, N, A, LDA, RWORK )
128 END IF
129 *
130 * Exit with RESID = 1/EPS if ANORM = 0.
131 *
132 EPS = SLAMCH( 'Epsilon' )
133 IF( ANORM.LE.ZERO ) THEN
134 RESID = ONE / EPS
135 RETURN
136 END IF
137 *
138 * Compute the maximum over the number of right hand sides of
139 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS )
140 *
141 RESID = ZERO
142 DO 10 J = 1, NRHS
143 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
144 CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
145 CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
146 BNORM = SCASUM( N, WORK, 1 )
147 XNORM = SCASUM( N, X( 1, J ), 1 )
148 IF( XNORM.LE.ZERO ) THEN
149 RESID = ONE / EPS
150 ELSE
151 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
152 END IF
153 10 CONTINUE
154 *
155 RETURN
156 *
157 * End of CTRT02
158 *
159 END