LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtrt03.f
Go to the documentation of this file.
1*> \brief \b DTRT03
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DTRT03( UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE,
12* CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER DIAG, TRANS, UPLO
16* INTEGER LDA, LDB, LDX, N, NRHS
17* DOUBLE PRECISION RESID, SCALE, TSCAL
18* ..
19* .. Array Arguments ..
20* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), CNORM( * ),
21* $ WORK( * ), X( LDX, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> DTRT03 computes the residual for the solution to a scaled triangular
31*> system of equations A*x = s*b or A'*x = s*b.
32*> Here A is a triangular matrix, A' is the transpose of A, s is a
33*> scalar, and x and b are N by NRHS matrices. The test ratio is the
34*> maximum over the number of right hand sides of
35*> norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
36*> where op(A) denotes A or A' and EPS is the machine epsilon.
37*> \endverbatim
38*
39* Arguments:
40* ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*> UPLO is CHARACTER*1
45*> Specifies whether the matrix A is upper or lower triangular.
46*> = 'U': Upper triangular
47*> = 'L': Lower triangular
48*> \endverbatim
49*>
50*> \param[in] TRANS
51*> \verbatim
52*> TRANS is CHARACTER*1
53*> Specifies the operation applied to A.
54*> = 'N': A *x = s*b (No transpose)
55*> = 'T': A'*x = s*b (Transpose)
56*> = 'C': A'*x = s*b (Conjugate transpose = Transpose)
57*> \endverbatim
58*>
59*> \param[in] DIAG
60*> \verbatim
61*> DIAG is CHARACTER*1
62*> Specifies whether or not the matrix A is unit triangular.
63*> = 'N': Non-unit triangular
64*> = 'U': Unit triangular
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*> N is INTEGER
70*> The order of the matrix A. N >= 0.
71*> \endverbatim
72*>
73*> \param[in] NRHS
74*> \verbatim
75*> NRHS is INTEGER
76*> The number of right hand sides, i.e., the number of columns
77*> of the matrices X and B. NRHS >= 0.
78*> \endverbatim
79*>
80*> \param[in] A
81*> \verbatim
82*> A is DOUBLE PRECISION array, dimension (LDA,N)
83*> The triangular matrix A. If UPLO = 'U', the leading n by n
84*> upper triangular part of the array A contains the upper
85*> triangular matrix, and the strictly lower triangular part of
86*> A is not referenced. If UPLO = 'L', the leading n by n lower
87*> triangular part of the array A contains the lower triangular
88*> matrix, and the strictly upper triangular part of A is not
89*> referenced. If DIAG = 'U', the diagonal elements of A are
90*> also not referenced and are assumed to be 1.
91*> \endverbatim
92*>
93*> \param[in] LDA
94*> \verbatim
95*> LDA is INTEGER
96*> The leading dimension of the array A. LDA >= max(1,N).
97*> \endverbatim
98*>
99*> \param[in] SCALE
100*> \verbatim
101*> SCALE is DOUBLE PRECISION
102*> The scaling factor s used in solving the triangular system.
103*> \endverbatim
104*>
105*> \param[in] CNORM
106*> \verbatim
107*> CNORM is DOUBLE PRECISION array, dimension (N)
108*> The 1-norms of the columns of A, not counting the diagonal.
109*> \endverbatim
110*>
111*> \param[in] TSCAL
112*> \verbatim
113*> TSCAL is DOUBLE PRECISION
114*> The scaling factor used in computing the 1-norms in CNORM.
115*> CNORM actually contains the column norms of TSCAL*A.
116*> \endverbatim
117*>
118*> \param[in] X
119*> \verbatim
120*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
121*> The computed solution vectors for the system of linear
122*> equations.
123*> \endverbatim
124*>
125*> \param[in] LDX
126*> \verbatim
127*> LDX is INTEGER
128*> The leading dimension of the array X. LDX >= max(1,N).
129*> \endverbatim
130*>
131*> \param[in] B
132*> \verbatim
133*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
134*> The right hand side vectors for the system of linear
135*> equations.
136*> \endverbatim
137*>
138*> \param[in] LDB
139*> \verbatim
140*> LDB is INTEGER
141*> The leading dimension of the array B. LDB >= max(1,N).
142*> \endverbatim
143*>
144*> \param[out] WORK
145*> \verbatim
146*> WORK is DOUBLE PRECISION array, dimension (N)
147*> \endverbatim
148*>
149*> \param[out] RESID
150*> \verbatim
151*> RESID is DOUBLE PRECISION
152*> The maximum over the number of right hand sides of
153*> norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
154*> \endverbatim
155*
156* Authors:
157* ========
158*
159*> \author Univ. of Tennessee
160*> \author Univ. of California Berkeley
161*> \author Univ. of Colorado Denver
162*> \author NAG Ltd.
163*
164*> \ingroup double_lin
165*
166* =====================================================================
167 SUBROUTINE dtrt03( UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE,
168 $ CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID )
169*
170* -- LAPACK test routine --
171* -- LAPACK is a software package provided by Univ. of Tennessee, --
172* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173*
174* .. Scalar Arguments ..
175 CHARACTER DIAG, TRANS, UPLO
176 INTEGER LDA, LDB, LDX, N, NRHS
177 DOUBLE PRECISION RESID, SCALE, TSCAL
178* ..
179* .. Array Arguments ..
180 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), CNORM( * ),
181 $ work( * ), x( ldx, * )
182* ..
183*
184* =====================================================================
185*
186* .. Parameters ..
187 DOUBLE PRECISION ONE, ZERO
188 parameter( one = 1.0d+0, zero = 0.0d+0 )
189* ..
190* .. Local Scalars ..
191 INTEGER IX, J
192 DOUBLE PRECISION BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
193* ..
194* .. External Functions ..
195 LOGICAL LSAME
196 INTEGER IDAMAX
197 DOUBLE PRECISION DLAMCH
198 EXTERNAL lsame, idamax, dlamch
199* ..
200* .. External Subroutines ..
201 EXTERNAL daxpy, dcopy, dlabad, dscal, dtrmv
202* ..
203* .. Intrinsic Functions ..
204 INTRINSIC abs, dble, max
205* ..
206* .. Executable Statements ..
207*
208* Quick exit if N = 0
209*
210 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
211 resid = zero
212 RETURN
213 END IF
214 eps = dlamch( 'Epsilon' )
215 smlnum = dlamch( 'Safe minimum' )
216 bignum = one / smlnum
217 CALL dlabad( smlnum, bignum )
218*
219* Compute the norm of the triangular matrix A using the column
220* norms already computed by DLATRS.
221*
222 tnorm = zero
223 IF( lsame( diag, 'N' ) ) THEN
224 DO 10 j = 1, n
225 tnorm = max( tnorm, tscal*abs( a( j, j ) )+cnorm( j ) )
226 10 CONTINUE
227 ELSE
228 DO 20 j = 1, n
229 tnorm = max( tnorm, tscal+cnorm( j ) )
230 20 CONTINUE
231 END IF
232*
233* Compute the maximum over the number of right hand sides of
234* norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
235*
236 resid = zero
237 DO 30 j = 1, nrhs
238 CALL dcopy( n, x( 1, j ), 1, work, 1 )
239 ix = idamax( n, work, 1 )
240 xnorm = max( one, abs( x( ix, j ) ) )
241 xscal = ( one / xnorm ) / dble( n )
242 CALL dscal( n, xscal, work, 1 )
243 CALL dtrmv( uplo, trans, diag, n, a, lda, work, 1 )
244 CALL daxpy( n, -scale*xscal, b( 1, j ), 1, work, 1 )
245 ix = idamax( n, work, 1 )
246 err = tscal*abs( work( ix ) )
247 ix = idamax( n, x( 1, j ), 1 )
248 xnorm = abs( x( ix, j ) )
249 IF( err*smlnum.LE.xnorm ) THEN
250 IF( xnorm.GT.zero )
251 $ err = err / xnorm
252 ELSE
253 IF( err.GT.zero )
254 $ err = one / eps
255 END IF
256 IF( err*smlnum.LE.tnorm ) THEN
257 IF( tnorm.GT.zero )
258 $ err = err / tnorm
259 ELSE
260 IF( err.GT.zero )
261 $ err = one / eps
262 END IF
263 resid = max( resid, err )
264 30 CONTINUE
265*
266 RETURN
267*
268* End of DTRT03
269*
270 END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:147
subroutine dtrt03(UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID)
DTRT03
Definition: dtrt03.f:169