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