1 SUBROUTINE CTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X,
2 $ LDX, B, 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 KD, LDAB, LDB, LDX, N, NRHS
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CTBT02 computes the residual for the computed solution to a
23 * triangular system of linear equations A*x = b, A**T *x = b, or
24 * A**H *x = b when A is a triangular band matrix. Here A**T denotes
25 * the transpose of A, A**H denotes the conjugate transpose of A, and
26 * x and b are N by NRHS matrices. The test ratio is the maximum over
27 * the number of right 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 * KD (input) INTEGER
54 * The number of superdiagonals or subdiagonals of the
55 * triangular band matrix A. KD >= 0.
56 *
57 * NRHS (input) INTEGER
58 * The number of right hand sides, i.e., the number of columns
59 * of the matrices X and B. NRHS >= 0.
60 *
61 * AB (input) COMPLEX array, dimension (LDA,N)
62 * The upper or lower triangular band matrix A, stored in the
63 * first kd+1 rows of the array. The j-th column of A is stored
64 * in the j-th column of the array AB as follows:
65 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
66 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
67 *
68 * LDAB (input) INTEGER
69 * The leading dimension of the array AB. LDAB >= max(1,KD+1).
70 *
71 * X (input) COMPLEX array, dimension (LDX,NRHS)
72 * The computed solution vectors for the system of linear
73 * equations.
74 *
75 * LDX (input) INTEGER
76 * The leading dimension of the array X. LDX >= max(1,N).
77 *
78 * B (input) COMPLEX array, dimension (LDB,NRHS)
79 * The right hand side vectors for the system of linear
80 * equations.
81 *
82 * LDB (input) INTEGER
83 * The leading dimension of the array B. LDB >= max(1,N).
84 *
85 * WORK (workspace) COMPLEX array, dimension (N)
86 *
87 * RWORK (workspace) REAL array, dimension (N)
88 *
89 * RESID (output) REAL
90 * The maximum over the number of right hand sides of
91 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
92 *
93 * =====================================================================
94 *
95 * .. Parameters ..
96 REAL ZERO, ONE
97 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
98 * ..
99 * .. Local Scalars ..
100 INTEGER J
101 REAL ANORM, BNORM, EPS, XNORM
102 * ..
103 * .. External Functions ..
104 LOGICAL LSAME
105 REAL CLANTB, SCASUM, SLAMCH
106 EXTERNAL LSAME, CLANTB, SCASUM, SLAMCH
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL CAXPY, CCOPY, CTBMV
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC CMPLX, MAX
113 * ..
114 * .. Executable Statements ..
115 *
116 * Quick exit if N = 0 or NRHS = 0
117 *
118 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
119 RESID = ZERO
120 RETURN
121 END IF
122 *
123 * Compute the 1-norm of A or A'.
124 *
125 IF( LSAME( TRANS, 'N' ) ) THEN
126 ANORM = CLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
127 ELSE
128 ANORM = CLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
129 END IF
130 *
131 * Exit with RESID = 1/EPS if ANORM = 0.
132 *
133 EPS = SLAMCH( 'Epsilon' )
134 IF( ANORM.LE.ZERO ) THEN
135 RESID = ONE / EPS
136 RETURN
137 END IF
138 *
139 * Compute the maximum over the number of right hand sides of
140 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
141 *
142 RESID = ZERO
143 DO 10 J = 1, NRHS
144 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
145 CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
146 CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
147 BNORM = SCASUM( N, WORK, 1 )
148 XNORM = SCASUM( N, X( 1, J ), 1 )
149 IF( XNORM.LE.ZERO ) THEN
150 RESID = ONE / EPS
151 ELSE
152 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
153 END IF
154 10 CONTINUE
155 *
156 RETURN
157 *
158 * End of CTBT02
159 *
160 END
2 $ LDX, B, 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 KD, LDAB, LDB, LDX, N, NRHS
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CTBT02 computes the residual for the computed solution to a
23 * triangular system of linear equations A*x = b, A**T *x = b, or
24 * A**H *x = b when A is a triangular band matrix. Here A**T denotes
25 * the transpose of A, A**H denotes the conjugate transpose of A, and
26 * x and b are N by NRHS matrices. The test ratio is the maximum over
27 * the number of right 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 * KD (input) INTEGER
54 * The number of superdiagonals or subdiagonals of the
55 * triangular band matrix A. KD >= 0.
56 *
57 * NRHS (input) INTEGER
58 * The number of right hand sides, i.e., the number of columns
59 * of the matrices X and B. NRHS >= 0.
60 *
61 * AB (input) COMPLEX array, dimension (LDA,N)
62 * The upper or lower triangular band matrix A, stored in the
63 * first kd+1 rows of the array. The j-th column of A is stored
64 * in the j-th column of the array AB as follows:
65 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
66 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
67 *
68 * LDAB (input) INTEGER
69 * The leading dimension of the array AB. LDAB >= max(1,KD+1).
70 *
71 * X (input) COMPLEX array, dimension (LDX,NRHS)
72 * The computed solution vectors for the system of linear
73 * equations.
74 *
75 * LDX (input) INTEGER
76 * The leading dimension of the array X. LDX >= max(1,N).
77 *
78 * B (input) COMPLEX array, dimension (LDB,NRHS)
79 * The right hand side vectors for the system of linear
80 * equations.
81 *
82 * LDB (input) INTEGER
83 * The leading dimension of the array B. LDB >= max(1,N).
84 *
85 * WORK (workspace) COMPLEX array, dimension (N)
86 *
87 * RWORK (workspace) REAL array, dimension (N)
88 *
89 * RESID (output) REAL
90 * The maximum over the number of right hand sides of
91 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
92 *
93 * =====================================================================
94 *
95 * .. Parameters ..
96 REAL ZERO, ONE
97 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
98 * ..
99 * .. Local Scalars ..
100 INTEGER J
101 REAL ANORM, BNORM, EPS, XNORM
102 * ..
103 * .. External Functions ..
104 LOGICAL LSAME
105 REAL CLANTB, SCASUM, SLAMCH
106 EXTERNAL LSAME, CLANTB, SCASUM, SLAMCH
107 * ..
108 * .. External Subroutines ..
109 EXTERNAL CAXPY, CCOPY, CTBMV
110 * ..
111 * .. Intrinsic Functions ..
112 INTRINSIC CMPLX, MAX
113 * ..
114 * .. Executable Statements ..
115 *
116 * Quick exit if N = 0 or NRHS = 0
117 *
118 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
119 RESID = ZERO
120 RETURN
121 END IF
122 *
123 * Compute the 1-norm of A or A'.
124 *
125 IF( LSAME( TRANS, 'N' ) ) THEN
126 ANORM = CLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
127 ELSE
128 ANORM = CLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
129 END IF
130 *
131 * Exit with RESID = 1/EPS if ANORM = 0.
132 *
133 EPS = SLAMCH( 'Epsilon' )
134 IF( ANORM.LE.ZERO ) THEN
135 RESID = ONE / EPS
136 RETURN
137 END IF
138 *
139 * Compute the maximum over the number of right hand sides of
140 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
141 *
142 RESID = ZERO
143 DO 10 J = 1, NRHS
144 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
145 CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
146 CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
147 BNORM = SCASUM( N, WORK, 1 )
148 XNORM = SCASUM( N, X( 1, J ), 1 )
149 IF( XNORM.LE.ZERO ) THEN
150 RESID = ONE / EPS
151 ELSE
152 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
153 END IF
154 10 CONTINUE
155 *
156 RETURN
157 *
158 * End of CTBT02
159 *
160 END