1 SUBROUTINE CTBT03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
2 $ SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
3 $ RESID )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER DIAG, TRANS, UPLO
11 INTEGER KD, LDAB, LDB, LDX, N, NRHS
12 REAL RESID, SCALE, TSCAL
13 * ..
14 * .. Array Arguments ..
15 REAL CNORM( * )
16 COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
17 $ X( LDX, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * CTBT03 computes the residual for the solution to a scaled triangular
24 * system of equations A*x = s*b, A**T *x = s*b, or A**H *x = s*b
25 * when A is a triangular band matrix. Here A**T denotes the transpose
26 * of A, A**H denotes the conjugate transpose of A, s is a scalar, and
27 * x and b are N by NRHS matrices. The test ratio is the maximum over
28 * the number of right hand sides of
29 * norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
30 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * Specifies whether the matrix A is upper or lower triangular.
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * TRANS (input) CHARACTER*1
41 * Specifies the operation applied to A.
42 * = 'N': A *x = s*b (No transpose)
43 * = 'T': A**T *x = s*b (Transpose)
44 * = 'C': A**H *x = s*b (Conjugate transpose)
45 *
46 * DIAG (input) CHARACTER*1
47 * Specifies whether or not the matrix A is unit triangular.
48 * = 'N': Non-unit triangular
49 * = 'U': Unit triangular
50 *
51 * N (input) INTEGER
52 * The order of the matrix A. N >= 0.
53 *
54 * KD (input) INTEGER
55 * The number of superdiagonals or subdiagonals of the
56 * triangular band matrix A. KD >= 0.
57 *
58 * NRHS (input) INTEGER
59 * The number of right hand sides, i.e., the number of columns
60 * of the matrices X and B. NRHS >= 0.
61 *
62 * AB (input) COMPLEX array, dimension (LDAB,N)
63 * The upper or lower triangular band matrix A, stored in the
64 * first kd+1 rows of the array. The j-th column of A is stored
65 * in the j-th column of the array AB as follows:
66 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
67 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
68 *
69 * LDAB (input) INTEGER
70 * The leading dimension of the array AB. LDAB >= KD+1.
71 *
72 * SCALE (input) REAL
73 * The scaling factor s used in solving the triangular system.
74 *
75 * CNORM (input) REAL array, dimension (N)
76 * The 1-norms of the columns of A, not counting the diagonal.
77 *
78 * TSCAL (input) REAL
79 * The scaling factor used in computing the 1-norms in CNORM.
80 * CNORM actually contains the column norms of TSCAL*A.
81 *
82 * X (input) COMPLEX array, dimension (LDX,NRHS)
83 * The computed solution vectors for the system of linear
84 * equations.
85 *
86 * LDX (input) INTEGER
87 * The leading dimension of the array X. LDX >= max(1,N).
88 *
89 * B (input) COMPLEX array, dimension (LDB,NRHS)
90 * The right hand side vectors for the system of linear
91 * equations.
92 *
93 * LDB (input) INTEGER
94 * The leading dimension of the array B. LDB >= max(1,N).
95 *
96 * WORK (workspace) COMPLEX array, dimension (N)
97 *
98 * RESID (output) REAL
99 * The maximum over the number of right hand sides of
100 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
101 *
102 * =====================================================================
103 *
104 *
105 * .. Parameters ..
106 REAL ONE, ZERO
107 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
108 * ..
109 * .. Local Scalars ..
110 INTEGER IX, J
111 REAL EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
112 * ..
113 * .. External Functions ..
114 LOGICAL LSAME
115 INTEGER ICAMAX
116 REAL SLAMCH
117 EXTERNAL LSAME, ICAMAX, SLAMCH
118 * ..
119 * .. External Subroutines ..
120 EXTERNAL CAXPY, CCOPY, CSSCAL, CTBMV
121 * ..
122 * .. Intrinsic Functions ..
123 INTRINSIC ABS, CMPLX, MAX, REAL
124 * ..
125 * .. Executable Statements ..
126 *
127 * Quick exit if N = 0
128 *
129 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
130 RESID = ZERO
131 RETURN
132 END IF
133 EPS = SLAMCH( 'Epsilon' )
134 SMLNUM = SLAMCH( 'Safe minimum' )
135 *
136 * Compute the norm of the triangular matrix A using the column
137 * norms already computed by CLATBS.
138 *
139 TNORM = ZERO
140 IF( LSAME( DIAG, 'N' ) ) THEN
141 IF( LSAME( UPLO, 'U' ) ) THEN
142 DO 10 J = 1, N
143 TNORM = MAX( TNORM, TSCAL*ABS( AB( KD+1, J ) )+
144 $ CNORM( J ) )
145 10 CONTINUE
146 ELSE
147 DO 20 J = 1, N
148 TNORM = MAX( TNORM, TSCAL*ABS( AB( 1, J ) )+CNORM( J ) )
149 20 CONTINUE
150 END IF
151 ELSE
152 DO 30 J = 1, N
153 TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
154 30 CONTINUE
155 END IF
156 *
157 * Compute the maximum over the number of right hand sides of
158 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
159 *
160 RESID = ZERO
161 DO 40 J = 1, NRHS
162 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
163 IX = ICAMAX( N, WORK, 1 )
164 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
165 XSCAL = ( ONE / XNORM ) / REAL( KD+1 )
166 CALL CSSCAL( N, XSCAL, WORK, 1 )
167 CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
168 CALL CAXPY( N, CMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
169 IX = ICAMAX( N, WORK, 1 )
170 ERR = TSCAL*ABS( WORK( IX ) )
171 IX = ICAMAX( N, X( 1, J ), 1 )
172 XNORM = ABS( X( IX, J ) )
173 IF( ERR*SMLNUM.LE.XNORM ) THEN
174 IF( XNORM.GT.ZERO )
175 $ ERR = ERR / XNORM
176 ELSE
177 IF( ERR.GT.ZERO )
178 $ ERR = ONE / EPS
179 END IF
180 IF( ERR*SMLNUM.LE.TNORM ) THEN
181 IF( TNORM.GT.ZERO )
182 $ ERR = ERR / TNORM
183 ELSE
184 IF( ERR.GT.ZERO )
185 $ ERR = ONE / EPS
186 END IF
187 RESID = MAX( RESID, ERR )
188 40 CONTINUE
189 *
190 RETURN
191 *
192 * End of CTBT03
193 *
194 END
2 $ SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
3 $ RESID )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER DIAG, TRANS, UPLO
11 INTEGER KD, LDAB, LDB, LDX, N, NRHS
12 REAL RESID, SCALE, TSCAL
13 * ..
14 * .. Array Arguments ..
15 REAL CNORM( * )
16 COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
17 $ X( LDX, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * CTBT03 computes the residual for the solution to a scaled triangular
24 * system of equations A*x = s*b, A**T *x = s*b, or A**H *x = s*b
25 * when A is a triangular band matrix. Here A**T denotes the transpose
26 * of A, A**H denotes the conjugate transpose of A, s is a scalar, and
27 * x and b are N by NRHS matrices. The test ratio is the maximum over
28 * the number of right hand sides of
29 * norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
30 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * Specifies whether the matrix A is upper or lower triangular.
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * TRANS (input) CHARACTER*1
41 * Specifies the operation applied to A.
42 * = 'N': A *x = s*b (No transpose)
43 * = 'T': A**T *x = s*b (Transpose)
44 * = 'C': A**H *x = s*b (Conjugate transpose)
45 *
46 * DIAG (input) CHARACTER*1
47 * Specifies whether or not the matrix A is unit triangular.
48 * = 'N': Non-unit triangular
49 * = 'U': Unit triangular
50 *
51 * N (input) INTEGER
52 * The order of the matrix A. N >= 0.
53 *
54 * KD (input) INTEGER
55 * The number of superdiagonals or subdiagonals of the
56 * triangular band matrix A. KD >= 0.
57 *
58 * NRHS (input) INTEGER
59 * The number of right hand sides, i.e., the number of columns
60 * of the matrices X and B. NRHS >= 0.
61 *
62 * AB (input) COMPLEX array, dimension (LDAB,N)
63 * The upper or lower triangular band matrix A, stored in the
64 * first kd+1 rows of the array. The j-th column of A is stored
65 * in the j-th column of the array AB as follows:
66 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
67 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
68 *
69 * LDAB (input) INTEGER
70 * The leading dimension of the array AB. LDAB >= KD+1.
71 *
72 * SCALE (input) REAL
73 * The scaling factor s used in solving the triangular system.
74 *
75 * CNORM (input) REAL array, dimension (N)
76 * The 1-norms of the columns of A, not counting the diagonal.
77 *
78 * TSCAL (input) REAL
79 * The scaling factor used in computing the 1-norms in CNORM.
80 * CNORM actually contains the column norms of TSCAL*A.
81 *
82 * X (input) COMPLEX array, dimension (LDX,NRHS)
83 * The computed solution vectors for the system of linear
84 * equations.
85 *
86 * LDX (input) INTEGER
87 * The leading dimension of the array X. LDX >= max(1,N).
88 *
89 * B (input) COMPLEX array, dimension (LDB,NRHS)
90 * The right hand side vectors for the system of linear
91 * equations.
92 *
93 * LDB (input) INTEGER
94 * The leading dimension of the array B. LDB >= max(1,N).
95 *
96 * WORK (workspace) COMPLEX array, dimension (N)
97 *
98 * RESID (output) REAL
99 * The maximum over the number of right hand sides of
100 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
101 *
102 * =====================================================================
103 *
104 *
105 * .. Parameters ..
106 REAL ONE, ZERO
107 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
108 * ..
109 * .. Local Scalars ..
110 INTEGER IX, J
111 REAL EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
112 * ..
113 * .. External Functions ..
114 LOGICAL LSAME
115 INTEGER ICAMAX
116 REAL SLAMCH
117 EXTERNAL LSAME, ICAMAX, SLAMCH
118 * ..
119 * .. External Subroutines ..
120 EXTERNAL CAXPY, CCOPY, CSSCAL, CTBMV
121 * ..
122 * .. Intrinsic Functions ..
123 INTRINSIC ABS, CMPLX, MAX, REAL
124 * ..
125 * .. Executable Statements ..
126 *
127 * Quick exit if N = 0
128 *
129 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
130 RESID = ZERO
131 RETURN
132 END IF
133 EPS = SLAMCH( 'Epsilon' )
134 SMLNUM = SLAMCH( 'Safe minimum' )
135 *
136 * Compute the norm of the triangular matrix A using the column
137 * norms already computed by CLATBS.
138 *
139 TNORM = ZERO
140 IF( LSAME( DIAG, 'N' ) ) THEN
141 IF( LSAME( UPLO, 'U' ) ) THEN
142 DO 10 J = 1, N
143 TNORM = MAX( TNORM, TSCAL*ABS( AB( KD+1, J ) )+
144 $ CNORM( J ) )
145 10 CONTINUE
146 ELSE
147 DO 20 J = 1, N
148 TNORM = MAX( TNORM, TSCAL*ABS( AB( 1, J ) )+CNORM( J ) )
149 20 CONTINUE
150 END IF
151 ELSE
152 DO 30 J = 1, N
153 TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
154 30 CONTINUE
155 END IF
156 *
157 * Compute the maximum over the number of right hand sides of
158 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
159 *
160 RESID = ZERO
161 DO 40 J = 1, NRHS
162 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
163 IX = ICAMAX( N, WORK, 1 )
164 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
165 XSCAL = ( ONE / XNORM ) / REAL( KD+1 )
166 CALL CSSCAL( N, XSCAL, WORK, 1 )
167 CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
168 CALL CAXPY( N, CMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
169 IX = ICAMAX( N, WORK, 1 )
170 ERR = TSCAL*ABS( WORK( IX ) )
171 IX = ICAMAX( N, X( 1, J ), 1 )
172 XNORM = ABS( X( IX, J ) )
173 IF( ERR*SMLNUM.LE.XNORM ) THEN
174 IF( XNORM.GT.ZERO )
175 $ ERR = ERR / XNORM
176 ELSE
177 IF( ERR.GT.ZERO )
178 $ ERR = ONE / EPS
179 END IF
180 IF( ERR*SMLNUM.LE.TNORM ) THEN
181 IF( TNORM.GT.ZERO )
182 $ ERR = ERR / TNORM
183 ELSE
184 IF( ERR.GT.ZERO )
185 $ ERR = ONE / EPS
186 END IF
187 RESID = MAX( RESID, ERR )
188 40 CONTINUE
189 *
190 RETURN
191 *
192 * End of CTBT03
193 *
194 END