LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtpt02.f
Go to the documentation of this file.
1*> \brief \b DTPT02
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 DTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
12* WORK, RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER DIAG, TRANS, UPLO
16* INTEGER LDB, LDX, N, NRHS
17* DOUBLE PRECISION RESID
18* ..
19* .. Array Arguments ..
20* DOUBLE PRECISION AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DTPT02 computes the residual for the computed solution to a
30*> triangular system of linear equations op(A)*X = B, when the
31*> triangular matrix A is stored in packed format. The test ratio is
32*> the maximum over
33*> norm(b - op(A)*x) / ( ||op(A)||_1 * norm(x) * EPS ),
34*> where op(A) = A or A**T, b is the column of B, x is the solution
35*> vector, and EPS is the machine epsilon.
36*> The norm used is the 1-norm.
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 = B (No transpose)
55*> = 'T': A**T * X = B (Transpose)
56*> = 'C': A**H * X = 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] AP
81*> \verbatim
82*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
83*> The upper or lower triangular matrix A, packed columnwise in
84*> a linear array. The j-th column of A is stored in the array
85*> AP as follows:
86*> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
87*> if UPLO = 'L',
88*> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
89*> \endverbatim
90*>
91*> \param[in] X
92*> \verbatim
93*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
94*> The computed solution vectors for the system of linear
95*> equations.
96*> \endverbatim
97*>
98*> \param[in] LDX
99*> \verbatim
100*> LDX is INTEGER
101*> The leading dimension of the array X. LDX >= max(1,N).
102*> \endverbatim
103*>
104*> \param[in] B
105*> \verbatim
106*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
107*> The right hand side vectors for the system of linear
108*> equations.
109*> \endverbatim
110*>
111*> \param[in] LDB
112*> \verbatim
113*> LDB is INTEGER
114*> The leading dimension of the array B. LDB >= max(1,N).
115*> \endverbatim
116*>
117*> \param[out] WORK
118*> \verbatim
119*> WORK is DOUBLE PRECISION array, dimension (N)
120*> \endverbatim
121*>
122*> \param[out] RESID
123*> \verbatim
124*> RESID is DOUBLE PRECISION
125*> The maximum over the number of right hand sides of
126*> norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ).
127*> \endverbatim
128*
129* Authors:
130* ========
131*
132*> \author Univ. of Tennessee
133*> \author Univ. of California Berkeley
134*> \author Univ. of Colorado Denver
135*> \author NAG Ltd.
136*
137*> \ingroup double_lin
138*
139* =====================================================================
140 SUBROUTINE dtpt02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
141 $ WORK, RESID )
142*
143* -- LAPACK test routine --
144* -- LAPACK is a software package provided by Univ. of Tennessee, --
145* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146*
147* .. Scalar Arguments ..
148 CHARACTER DIAG, TRANS, UPLO
149 INTEGER LDB, LDX, N, NRHS
150 DOUBLE PRECISION RESID
151* ..
152* .. Array Arguments ..
153 DOUBLE PRECISION AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 DOUBLE PRECISION ZERO, ONE
160 parameter( zero = 0.0d+0, one = 1.0d+0 )
161* ..
162* .. Local Scalars ..
163 INTEGER J
164 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
165* ..
166* .. External Functions ..
167 LOGICAL LSAME
168 DOUBLE PRECISION DASUM, DLAMCH, DLANTP
169 EXTERNAL lsame, dasum, dlamch, dlantp
170* ..
171* .. External Subroutines ..
172 EXTERNAL daxpy, dcopy, dtpmv
173* ..
174* .. Intrinsic Functions ..
175 INTRINSIC max
176* ..
177* .. Executable Statements ..
178*
179* Quick exit if N = 0 or NRHS = 0
180*
181 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
182 resid = zero
183 RETURN
184 END IF
185*
186* Compute the 1-norm of op(A).
187*
188 IF( lsame( trans, 'N' ) ) THEN
189 anorm = dlantp( '1', uplo, diag, n, ap, work )
190 ELSE
191 anorm = dlantp( 'I', uplo, diag, n, ap, work )
192 END IF
193*
194* Exit with RESID = 1/EPS if ANORM = 0.
195*
196 eps = dlamch( 'Epsilon' )
197 IF( anorm.LE.zero ) THEN
198 resid = one / eps
199 RETURN
200 END IF
201*
202* Compute the maximum over the number of right hand sides of
203* norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ).
204*
205 resid = zero
206 DO 10 j = 1, nrhs
207 CALL dcopy( n, x( 1, j ), 1, work, 1 )
208 CALL dtpmv( uplo, trans, diag, n, ap, work, 1 )
209 CALL daxpy( n, -one, b( 1, j ), 1, work, 1 )
210 bnorm = dasum( n, work, 1 )
211 xnorm = dasum( n, x( 1, j ), 1 )
212 IF( xnorm.LE.zero ) THEN
213 resid = one / eps
214 ELSE
215 resid = max( resid, ( ( bnorm / anorm ) / xnorm ) / eps )
216 END IF
217 10 CONTINUE
218*
219 RETURN
220*
221* End of DTPT02
222*
223 END
subroutine dtpt02(uplo, trans, diag, n, nrhs, ap, x, ldx, b, ldb, work, resid)
DTPT02
Definition dtpt02.f:142
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dtpmv(uplo, trans, diag, n, ap, x, incx)
DTPMV
Definition dtpmv.f:142