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