1 SUBROUTINE STRT03( UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE,
2 $ CNORM, TSCAL, X, 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 LDA, LDB, LDX, N, NRHS
11 REAL RESID, SCALE, TSCAL
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), B( LDB, * ), CNORM( * ),
15 $ WORK( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * STRT03 computes the residual for the solution to a scaled triangular
22 * system of equations A*x = s*b or A'*x = s*b.
23 * Here A is a triangular matrix, A' is the transpose of A, s is a
24 * scalar, and x and b are N by NRHS matrices. The test ratio is the
25 * maximum over the number of right hand sides of
26 * norm(s*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 = s*b (No transpose)
40 * = 'T': A'*x = s*b (Transpose)
41 * = 'C': A'*x = s*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) REAL 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 * SCALE (input) REAL
69 * The scaling factor s used in solving the triangular system.
70 *
71 * CNORM (input) REAL array, dimension (N)
72 * The 1-norms of the columns of A, not counting the diagonal.
73 *
74 * TSCAL (input) REAL
75 * The scaling factor used in computing the 1-norms in CNORM.
76 * CNORM actually contains the column norms of TSCAL*A.
77 *
78 * X (input) REAL array, dimension (LDX,NRHS)
79 * The computed solution vectors for the system of linear
80 * equations.
81 *
82 * LDX (input) INTEGER
83 * The leading dimension of the array X. LDX >= max(1,N).
84 *
85 * B (input) REAL array, dimension (LDB,NRHS)
86 * The right hand side vectors for the system of linear
87 * equations.
88 *
89 * LDB (input) INTEGER
90 * The leading dimension of the array B. LDB >= max(1,N).
91 *
92 * WORK (workspace) REAL array, dimension (N)
93 *
94 * RESID (output) REAL
95 * The maximum over the number of right hand sides of
96 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 REAL ONE, ZERO
102 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER IX, J
106 REAL BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 INTEGER ISAMAX
111 REAL SLAMCH
112 EXTERNAL LSAME, ISAMAX, SLAMCH
113 * ..
114 * .. External Subroutines ..
115 EXTERNAL SAXPY, SCOPY, SLABAD, SSCAL, STRMV
116 * ..
117 * .. Intrinsic Functions ..
118 INTRINSIC ABS, MAX, REAL
119 * ..
120 * .. Executable Statements ..
121 *
122 * Quick exit if N = 0
123 *
124 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
125 RESID = ZERO
126 RETURN
127 END IF
128 EPS = SLAMCH( 'Epsilon' )
129 SMLNUM = SLAMCH( 'Safe minimum' )
130 BIGNUM = ONE / SMLNUM
131 CALL SLABAD( SMLNUM, BIGNUM )
132 *
133 * Compute the norm of the triangular matrix A using the column
134 * norms already computed by SLATRS.
135 *
136 TNORM = ZERO
137 IF( LSAME( DIAG, 'N' ) ) THEN
138 DO 10 J = 1, N
139 TNORM = MAX( TNORM, TSCAL*ABS( A( J, J ) )+CNORM( J ) )
140 10 CONTINUE
141 ELSE
142 DO 20 J = 1, N
143 TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
144 20 CONTINUE
145 END IF
146 *
147 * Compute the maximum over the number of right hand sides of
148 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
149 *
150 RESID = ZERO
151 DO 30 J = 1, NRHS
152 CALL SCOPY( N, X( 1, J ), 1, WORK, 1 )
153 IX = ISAMAX( N, WORK, 1 )
154 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
155 XSCAL = ( ONE / XNORM ) / REAL( N )
156 CALL SSCAL( N, XSCAL, WORK, 1 )
157 CALL STRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
158 CALL SAXPY( N, -SCALE*XSCAL, B( 1, J ), 1, WORK, 1 )
159 IX = ISAMAX( N, WORK, 1 )
160 ERR = TSCAL*ABS( WORK( IX ) )
161 IX = ISAMAX( N, X( 1, J ), 1 )
162 XNORM = ABS( X( IX, J ) )
163 IF( ERR*SMLNUM.LE.XNORM ) THEN
164 IF( XNORM.GT.ZERO )
165 $ ERR = ERR / XNORM
166 ELSE
167 IF( ERR.GT.ZERO )
168 $ ERR = ONE / EPS
169 END IF
170 IF( ERR*SMLNUM.LE.TNORM ) THEN
171 IF( TNORM.GT.ZERO )
172 $ ERR = ERR / TNORM
173 ELSE
174 IF( ERR.GT.ZERO )
175 $ ERR = ONE / EPS
176 END IF
177 RESID = MAX( RESID, ERR )
178 30 CONTINUE
179 *
180 RETURN
181 *
182 * End of STRT03
183 *
184 END
2 $ CNORM, TSCAL, X, 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 LDA, LDB, LDX, N, NRHS
11 REAL RESID, SCALE, TSCAL
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), B( LDB, * ), CNORM( * ),
15 $ WORK( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * STRT03 computes the residual for the solution to a scaled triangular
22 * system of equations A*x = s*b or A'*x = s*b.
23 * Here A is a triangular matrix, A' is the transpose of A, s is a
24 * scalar, and x and b are N by NRHS matrices. The test ratio is the
25 * maximum over the number of right hand sides of
26 * norm(s*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 = s*b (No transpose)
40 * = 'T': A'*x = s*b (Transpose)
41 * = 'C': A'*x = s*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) REAL 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 * SCALE (input) REAL
69 * The scaling factor s used in solving the triangular system.
70 *
71 * CNORM (input) REAL array, dimension (N)
72 * The 1-norms of the columns of A, not counting the diagonal.
73 *
74 * TSCAL (input) REAL
75 * The scaling factor used in computing the 1-norms in CNORM.
76 * CNORM actually contains the column norms of TSCAL*A.
77 *
78 * X (input) REAL array, dimension (LDX,NRHS)
79 * The computed solution vectors for the system of linear
80 * equations.
81 *
82 * LDX (input) INTEGER
83 * The leading dimension of the array X. LDX >= max(1,N).
84 *
85 * B (input) REAL array, dimension (LDB,NRHS)
86 * The right hand side vectors for the system of linear
87 * equations.
88 *
89 * LDB (input) INTEGER
90 * The leading dimension of the array B. LDB >= max(1,N).
91 *
92 * WORK (workspace) REAL array, dimension (N)
93 *
94 * RESID (output) REAL
95 * The maximum over the number of right hand sides of
96 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 REAL ONE, ZERO
102 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER IX, J
106 REAL BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 INTEGER ISAMAX
111 REAL SLAMCH
112 EXTERNAL LSAME, ISAMAX, SLAMCH
113 * ..
114 * .. External Subroutines ..
115 EXTERNAL SAXPY, SCOPY, SLABAD, SSCAL, STRMV
116 * ..
117 * .. Intrinsic Functions ..
118 INTRINSIC ABS, MAX, REAL
119 * ..
120 * .. Executable Statements ..
121 *
122 * Quick exit if N = 0
123 *
124 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
125 RESID = ZERO
126 RETURN
127 END IF
128 EPS = SLAMCH( 'Epsilon' )
129 SMLNUM = SLAMCH( 'Safe minimum' )
130 BIGNUM = ONE / SMLNUM
131 CALL SLABAD( SMLNUM, BIGNUM )
132 *
133 * Compute the norm of the triangular matrix A using the column
134 * norms already computed by SLATRS.
135 *
136 TNORM = ZERO
137 IF( LSAME( DIAG, 'N' ) ) THEN
138 DO 10 J = 1, N
139 TNORM = MAX( TNORM, TSCAL*ABS( A( J, J ) )+CNORM( J ) )
140 10 CONTINUE
141 ELSE
142 DO 20 J = 1, N
143 TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
144 20 CONTINUE
145 END IF
146 *
147 * Compute the maximum over the number of right hand sides of
148 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
149 *
150 RESID = ZERO
151 DO 30 J = 1, NRHS
152 CALL SCOPY( N, X( 1, J ), 1, WORK, 1 )
153 IX = ISAMAX( N, WORK, 1 )
154 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
155 XSCAL = ( ONE / XNORM ) / REAL( N )
156 CALL SSCAL( N, XSCAL, WORK, 1 )
157 CALL STRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
158 CALL SAXPY( N, -SCALE*XSCAL, B( 1, J ), 1, WORK, 1 )
159 IX = ISAMAX( N, WORK, 1 )
160 ERR = TSCAL*ABS( WORK( IX ) )
161 IX = ISAMAX( N, X( 1, J ), 1 )
162 XNORM = ABS( X( IX, J ) )
163 IF( ERR*SMLNUM.LE.XNORM ) THEN
164 IF( XNORM.GT.ZERO )
165 $ ERR = ERR / XNORM
166 ELSE
167 IF( ERR.GT.ZERO )
168 $ ERR = ONE / EPS
169 END IF
170 IF( ERR*SMLNUM.LE.TNORM ) THEN
171 IF( TNORM.GT.ZERO )
172 $ ERR = ERR / TNORM
173 ELSE
174 IF( ERR.GT.ZERO )
175 $ ERR = ONE / EPS
176 END IF
177 RESID = MAX( RESID, ERR )
178 30 CONTINUE
179 *
180 RETURN
181 *
182 * End of STRT03
183 *
184 END