1 SUBROUTINE ZTRT03( 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 DOUBLE PRECISION RESID, SCALE, TSCAL
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION CNORM( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZTRT03 computes the residual for the solution to a scaled triangular
23 * system of equations A*x = s*b, A**T *x = s*b, or A**H *x = s*b.
24 * Here A is a triangular matrix, A**T denotes the transpose of A, A**H
25 * denotes the conjugate transpose of A, s is a scalar, and x and b are
26 * N by NRHS matrices. The test ratio is the maximum over the number of
27 * right hand sides of
28 * norm(s*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 = s*b (No transpose)
42 * = 'T': A**T *x = s*b (Transpose)
43 * = 'C': A**H *x = s*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 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * A (input) COMPLEX*16 array, dimension (LDA,N)
58 * The triangular matrix A. If UPLO = 'U', the leading n by n
59 * upper triangular part of the array A contains the upper
60 * triangular matrix, and the strictly lower triangular part of
61 * A is not referenced. If UPLO = 'L', the leading n by n lower
62 * triangular part of the array A contains the lower triangular
63 * matrix, and the strictly upper triangular part of A is not
64 * referenced. If DIAG = 'U', the diagonal elements of A are
65 * also not referenced and are assumed to be 1.
66 *
67 * LDA (input) INTEGER
68 * The leading dimension of the array A. LDA >= max(1,N).
69 *
70 * SCALE (input) DOUBLE PRECISION
71 * The scaling factor s used in solving the triangular system.
72 *
73 * CNORM (input) DOUBLE PRECISION array, dimension (N)
74 * The 1-norms of the columns of A, not counting the diagonal.
75 *
76 * TSCAL (input) DOUBLE PRECISION
77 * The scaling factor used in computing the 1-norms in CNORM.
78 * CNORM actually contains the column norms of TSCAL*A.
79 *
80 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
81 * The computed solution vectors for the system of linear
82 * equations.
83 *
84 * LDX (input) INTEGER
85 * The leading dimension of the array X. LDX >= max(1,N).
86 *
87 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
88 * The right hand side vectors for the system of linear
89 * equations.
90 *
91 * LDB (input) INTEGER
92 * The leading dimension of the array B. LDB >= max(1,N).
93 *
94 * WORK (workspace) COMPLEX*16 array, dimension (N)
95 *
96 * RESID (output) DOUBLE PRECISION
97 * The maximum over the number of right hand sides of
98 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ONE, ZERO
104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER IX, J
108 DOUBLE PRECISION EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
109 * ..
110 * .. External Functions ..
111 LOGICAL LSAME
112 INTEGER IZAMAX
113 DOUBLE PRECISION DLAMCH
114 EXTERNAL LSAME, IZAMAX, DLAMCH
115 * ..
116 * .. External Subroutines ..
117 EXTERNAL ZAXPY, ZCOPY, ZDSCAL, ZTRMV
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC ABS, DBLE, DCMPLX, MAX
121 * ..
122 * .. Executable Statements ..
123 *
124 * Quick exit if N = 0
125 *
126 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
127 RESID = ZERO
128 RETURN
129 END IF
130 EPS = DLAMCH( 'Epsilon' )
131 SMLNUM = DLAMCH( 'Safe minimum' )
132 *
133 * Compute the norm of the triangular matrix A using the column
134 * norms already computed by ZLATRS.
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 ZCOPY( N, X( 1, J ), 1, WORK, 1 )
153 IX = IZAMAX( N, WORK, 1 )
154 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
155 XSCAL = ( ONE / XNORM ) / DBLE( N )
156 CALL ZDSCAL( N, XSCAL, WORK, 1 )
157 CALL ZTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
158 CALL ZAXPY( N, DCMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
159 IX = IZAMAX( N, WORK, 1 )
160 ERR = TSCAL*ABS( WORK( IX ) )
161 IX = IZAMAX( 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 ZTRT03
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 DOUBLE PRECISION RESID, SCALE, TSCAL
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION CNORM( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ),
16 $ X( LDX, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZTRT03 computes the residual for the solution to a scaled triangular
23 * system of equations A*x = s*b, A**T *x = s*b, or A**H *x = s*b.
24 * Here A is a triangular matrix, A**T denotes the transpose of A, A**H
25 * denotes the conjugate transpose of A, s is a scalar, and x and b are
26 * N by NRHS matrices. The test ratio is the maximum over the number of
27 * right hand sides of
28 * norm(s*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 = s*b (No transpose)
42 * = 'T': A**T *x = s*b (Transpose)
43 * = 'C': A**H *x = s*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 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * A (input) COMPLEX*16 array, dimension (LDA,N)
58 * The triangular matrix A. If UPLO = 'U', the leading n by n
59 * upper triangular part of the array A contains the upper
60 * triangular matrix, and the strictly lower triangular part of
61 * A is not referenced. If UPLO = 'L', the leading n by n lower
62 * triangular part of the array A contains the lower triangular
63 * matrix, and the strictly upper triangular part of A is not
64 * referenced. If DIAG = 'U', the diagonal elements of A are
65 * also not referenced and are assumed to be 1.
66 *
67 * LDA (input) INTEGER
68 * The leading dimension of the array A. LDA >= max(1,N).
69 *
70 * SCALE (input) DOUBLE PRECISION
71 * The scaling factor s used in solving the triangular system.
72 *
73 * CNORM (input) DOUBLE PRECISION array, dimension (N)
74 * The 1-norms of the columns of A, not counting the diagonal.
75 *
76 * TSCAL (input) DOUBLE PRECISION
77 * The scaling factor used in computing the 1-norms in CNORM.
78 * CNORM actually contains the column norms of TSCAL*A.
79 *
80 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
81 * The computed solution vectors for the system of linear
82 * equations.
83 *
84 * LDX (input) INTEGER
85 * The leading dimension of the array X. LDX >= max(1,N).
86 *
87 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
88 * The right hand side vectors for the system of linear
89 * equations.
90 *
91 * LDB (input) INTEGER
92 * The leading dimension of the array B. LDB >= max(1,N).
93 *
94 * WORK (workspace) COMPLEX*16 array, dimension (N)
95 *
96 * RESID (output) DOUBLE PRECISION
97 * The maximum over the number of right hand sides of
98 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ONE, ZERO
104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER IX, J
108 DOUBLE PRECISION EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
109 * ..
110 * .. External Functions ..
111 LOGICAL LSAME
112 INTEGER IZAMAX
113 DOUBLE PRECISION DLAMCH
114 EXTERNAL LSAME, IZAMAX, DLAMCH
115 * ..
116 * .. External Subroutines ..
117 EXTERNAL ZAXPY, ZCOPY, ZDSCAL, ZTRMV
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC ABS, DBLE, DCMPLX, MAX
121 * ..
122 * .. Executable Statements ..
123 *
124 * Quick exit if N = 0
125 *
126 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
127 RESID = ZERO
128 RETURN
129 END IF
130 EPS = DLAMCH( 'Epsilon' )
131 SMLNUM = DLAMCH( 'Safe minimum' )
132 *
133 * Compute the norm of the triangular matrix A using the column
134 * norms already computed by ZLATRS.
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 ZCOPY( N, X( 1, J ), 1, WORK, 1 )
153 IX = IZAMAX( N, WORK, 1 )
154 XNORM = MAX( ONE, ABS( X( IX, J ) ) )
155 XSCAL = ( ONE / XNORM ) / DBLE( N )
156 CALL ZDSCAL( N, XSCAL, WORK, 1 )
157 CALL ZTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
158 CALL ZAXPY( N, DCMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
159 IX = IZAMAX( N, WORK, 1 )
160 ERR = TSCAL*ABS( WORK( IX ) )
161 IX = IZAMAX( 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 ZTRT03
183 *
184 END