LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
stpt03.f
Go to the documentation of this file.
1*> \brief \b STPT03
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 STPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM,
12* TSCAL, X, LDX, B, LDB, WORK, RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER DIAG, TRANS, UPLO
16* INTEGER LDB, LDX, N, NRHS
17* REAL RESID, SCALE, TSCAL
18* ..
19* .. Array Arguments ..
20* REAL AP( * ), B( LDB, * ), CNORM( * ), WORK( * ),
21* \$ X( LDX, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> STPT03 computes the residual for the solution to a scaled triangular
31*> system of equations A*x = s*b or A'*x = s*b when the triangular
32*> matrix A is stored in packed format. Here A' is the transpose of A,
33*> s is a scalar, and x and b are N by NRHS matrices. The test ratio is
34*> the 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] AP
81*> \verbatim
82*> AP is REAL 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] SCALE
92*> \verbatim
93*> SCALE is REAL
94*> The scaling factor s used in solving the triangular system.
95*> \endverbatim
96*>
97*> \param[in] CNORM
98*> \verbatim
99*> CNORM is REAL array, dimension (N)
100*> The 1-norms of the columns of A, not counting the diagonal.
101*> \endverbatim
102*>
103*> \param[in] TSCAL
104*> \verbatim
105*> TSCAL is REAL
106*> The scaling factor used in computing the 1-norms in CNORM.
107*> CNORM actually contains the column norms of TSCAL*A.
108*> \endverbatim
109*>
110*> \param[in] X
111*> \verbatim
112*> X is REAL array, dimension (LDX,NRHS)
113*> The computed solution vectors for the system of linear
114*> equations.
115*> \endverbatim
116*>
117*> \param[in] LDX
118*> \verbatim
119*> LDX is INTEGER
120*> The leading dimension of the array X. LDX >= max(1,N).
121*> \endverbatim
122*>
123*> \param[in] B
124*> \verbatim
125*> B is REAL array, dimension (LDB,NRHS)
126*> The right hand side vectors for the system of linear
127*> equations.
128*> \endverbatim
129*>
130*> \param[in] LDB
131*> \verbatim
132*> LDB is INTEGER
133*> The leading dimension of the array B. LDB >= max(1,N).
134*> \endverbatim
135*>
136*> \param[out] WORK
137*> \verbatim
138*> WORK is REAL array, dimension (N)
139*> \endverbatim
140*>
141*> \param[out] RESID
142*> \verbatim
143*> RESID is REAL
144*> The maximum over the number of right hand sides of
145*> norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
146*> \endverbatim
147*
148* Authors:
149* ========
150*
151*> \author Univ. of Tennessee
152*> \author Univ. of California Berkeley
153*> \author Univ. of Colorado Denver
154*> \author NAG Ltd.
155*
156*> \ingroup single_lin
157*
158* =====================================================================
159 SUBROUTINE stpt03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM,
160 \$ TSCAL, X, LDX, B, LDB, WORK, RESID )
161*
162* -- LAPACK test routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 CHARACTER DIAG, TRANS, UPLO
168 INTEGER LDB, LDX, N, NRHS
169 REAL RESID, SCALE, TSCAL
170* ..
171* .. Array Arguments ..
172 REAL AP( * ), B( LDB, * ), CNORM( * ), WORK( * ),
173 \$ x( ldx, * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ONE, ZERO
180 parameter( one = 1.0e+0, zero = 0.0e+0 )
181* ..
182* .. Local Scalars ..
183 INTEGER IX, J, JJ
184 REAL BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
185* ..
186* .. External Functions ..
187 LOGICAL LSAME
188 INTEGER ISAMAX
189 REAL SLAMCH
190 EXTERNAL lsame, isamax, slamch
191* ..
192* .. External Subroutines ..
193 EXTERNAL saxpy, scopy, slabad, sscal, stpmv
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC abs, max, real
197* ..
198* .. Executable Statements ..
199*
200* Quick exit if N = 0.
201*
202 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
203 resid = zero
204 RETURN
205 END IF
206 eps = slamch( 'Epsilon' )
207 smlnum = slamch( 'Safe minimum' )
208 bignum = one / smlnum
209 CALL slabad( smlnum, bignum )
210*
211* Compute the norm of the triangular matrix A using the column
212* norms already computed by SLATPS.
213*
214 tnorm = zero
215 IF( lsame( diag, 'N' ) ) THEN
216 IF( lsame( uplo, 'U' ) ) THEN
217 jj = 1
218 DO 10 j = 1, n
219 tnorm = max( tnorm, tscal*abs( ap( jj ) )+cnorm( j ) )
220 jj = jj + j + 1
221 10 CONTINUE
222 ELSE
223 jj = 1
224 DO 20 j = 1, n
225 tnorm = max( tnorm, tscal*abs( ap( jj ) )+cnorm( j ) )
226 jj = jj + n - j + 1
227 20 CONTINUE
228 END IF
229 ELSE
230 DO 30 j = 1, n
231 tnorm = max( tnorm, tscal+cnorm( j ) )
232 30 CONTINUE
233 END IF
234*
235* Compute the maximum over the number of right hand sides of
236* norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
237*
238 resid = zero
239 DO 40 j = 1, nrhs
240 CALL scopy( n, x( 1, j ), 1, work, 1 )
241 ix = isamax( n, work, 1 )
242 xnorm = max( one, abs( x( ix, j ) ) )
243 xscal = ( one / xnorm ) / real( n )
244 CALL sscal( n, xscal, work, 1 )
245 CALL stpmv( uplo, trans, diag, n, ap, work, 1 )
246 CALL saxpy( n, -scale*xscal, b( 1, j ), 1, work, 1 )
247 ix = isamax( n, work, 1 )
248 err = tscal*abs( work( ix ) )
249 ix = isamax( n, x( 1, j ), 1 )
250 xnorm = abs( x( ix, j ) )
251 IF( err*smlnum.LE.xnorm ) THEN
252 IF( xnorm.GT.zero )
253 \$ err = err / xnorm
254 ELSE
255 IF( err.GT.zero )
256 \$ err = one / eps
257 END IF
258 IF( err*smlnum.LE.tnorm ) THEN
259 IF( tnorm.GT.zero )
260 \$ err = err / tnorm
261 ELSE
262 IF( err.GT.zero )
263 \$ err = one / eps
264 END IF
265 resid = max( resid, err )
266 40 CONTINUE
267*
268 RETURN
269*
270* End of STPT03
271*
272 END
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine stpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
STPMV
Definition: stpmv.f:142
subroutine stpt03(UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID)
STPT03
Definition: stpt03.f:161