LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 A*x = b or A'*x = b when
31 *> the triangular matrix A is stored in packed format. Here A' is the
32 *> transpose of A and x and b are N by NRHS matrices. The test ratio is
33 *> the maximum over the number of right hand sides of
34 *> norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
35 *> where op(A) denotes A or A' and EPS is the machine epsilon.
36 *> \endverbatim
37 *
38 * Arguments:
39 * ==========
40 *
41 *> \param[in] UPLO
42 *> \verbatim
43 *> UPLO is CHARACTER*1
44 *> Specifies whether the matrix A is upper or lower triangular.
45 *> = 'U': Upper triangular
46 *> = 'L': Lower triangular
47 *> \endverbatim
48 *>
49 *> \param[in] TRANS
50 *> \verbatim
51 *> TRANS is CHARACTER*1
52 *> Specifies the operation applied to A.
53 *> = 'N': A *x = b (No transpose)
54 *> = 'T': A'*x = b (Transpose)
55 *> = 'C': A'*x = b (Conjugate transpose = Transpose)
56 *> \endverbatim
57 *>
58 *> \param[in] DIAG
59 *> \verbatim
60 *> DIAG is CHARACTER*1
61 *> Specifies whether or not the matrix A is unit triangular.
62 *> = 'N': Non-unit triangular
63 *> = 'U': Unit triangular
64 *> \endverbatim
65 *>
66 *> \param[in] N
67 *> \verbatim
68 *> N is INTEGER
69 *> The order of the matrix A. N >= 0.
70 *> \endverbatim
71 *>
72 *> \param[in] NRHS
73 *> \verbatim
74 *> NRHS is INTEGER
75 *> The number of right hand sides, i.e., the number of columns
76 *> of the matrices X and B. NRHS >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in] AP
80 *> \verbatim
81 *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
82 *> The upper or lower triangular matrix A, packed columnwise in
83 *> a linear array. The j-th column of A is stored in the array
84 *> AP as follows:
85 *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
86 *> if UPLO = 'L',
87 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
88 *> \endverbatim
89 *>
90 *> \param[in] X
91 *> \verbatim
92 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
93 *> The computed solution vectors for the system of linear
94 *> equations.
95 *> \endverbatim
96 *>
97 *> \param[in] LDX
98 *> \verbatim
99 *> LDX is INTEGER
100 *> The leading dimension of the array X. LDX >= max(1,N).
101 *> \endverbatim
102 *>
103 *> \param[in] B
104 *> \verbatim
105 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
106 *> The right hand side vectors for the system of linear
107 *> equations.
108 *> \endverbatim
109 *>
110 *> \param[in] LDB
111 *> \verbatim
112 *> LDB is INTEGER
113 *> The leading dimension of the array B. LDB >= max(1,N).
114 *> \endverbatim
115 *>
116 *> \param[out] WORK
117 *> \verbatim
118 *> WORK is DOUBLE PRECISION array, dimension (N)
119 *> \endverbatim
120 *>
121 *> \param[out] RESID
122 *> \verbatim
123 *> RESID is DOUBLE PRECISION
124 *> The maximum over the number of right hand sides of
125 *> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
126 *> \endverbatim
127 *
128 * Authors:
129 * ========
130 *
131 *> \author Univ. of Tennessee
132 *> \author Univ. of California Berkeley
133 *> \author Univ. of Colorado Denver
134 *> \author NAG Ltd.
135 *
136 *> \date November 2011
137 *
138 *> \ingroup double_lin
139 *
140 * =====================================================================
141  SUBROUTINE dtpt02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
142  $ work, resid )
143 *
144 * -- LAPACK test routine (version 3.4.0) --
145 * -- LAPACK is a software package provided by Univ. of Tennessee, --
146 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147 * November 2011
148 *
149 * .. Scalar Arguments ..
150  CHARACTER DIAG, TRANS, UPLO
151  INTEGER LDB, LDX, N, NRHS
152  DOUBLE PRECISION RESID
153 * ..
154 * .. Array Arguments ..
155  DOUBLE PRECISION AP( * ), B( ldb, * ), WORK( * ), X( ldx, * )
156 * ..
157 *
158 * =====================================================================
159 *
160 * .. Parameters ..
161  DOUBLE PRECISION ZERO, ONE
162  parameter ( zero = 0.0d+0, one = 1.0d+0 )
163 * ..
164 * .. Local Scalars ..
165  INTEGER J
166  DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
167 * ..
168 * .. External Functions ..
169  LOGICAL LSAME
170  DOUBLE PRECISION DASUM, DLAMCH, DLANTP
171  EXTERNAL lsame, dasum, dlamch, dlantp
172 * ..
173 * .. External Subroutines ..
174  EXTERNAL daxpy, dcopy, dtpmv
175 * ..
176 * .. Intrinsic Functions ..
177  INTRINSIC max
178 * ..
179 * .. Executable Statements ..
180 *
181 * Quick exit if N = 0 or NRHS = 0
182 *
183  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
184  resid = zero
185  RETURN
186  END IF
187 *
188 * Compute the 1-norm of A or A'.
189 *
190  IF( lsame( trans, 'N' ) ) THEN
191  anorm = dlantp( '1', uplo, diag, n, ap, work )
192  ELSE
193  anorm = dlantp( 'I', uplo, diag, n, ap, work )
194  END IF
195 *
196 * Exit with RESID = 1/EPS if ANORM = 0.
197 *
198  eps = dlamch( 'Epsilon' )
199  IF( anorm.LE.zero ) THEN
200  resid = one / eps
201  RETURN
202  END IF
203 *
204 * Compute the maximum over the number of right hand sides of
205 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
206 *
207  resid = zero
208  DO 10 j = 1, nrhs
209  CALL dcopy( n, x( 1, j ), 1, work, 1 )
210  CALL dtpmv( uplo, trans, diag, n, ap, work, 1 )
211  CALL daxpy( n, -one, b( 1, j ), 1, work, 1 )
212  bnorm = dasum( n, work, 1 )
213  xnorm = dasum( n, x( 1, j ), 1 )
214  IF( xnorm.LE.zero ) THEN
215  resid = one / eps
216  ELSE
217  resid = max( resid, ( ( bnorm / anorm ) / xnorm ) / eps )
218  END IF
219  10 CONTINUE
220 *
221  RETURN
222 *
223 * End of DTPT02
224 *
225  END
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dtpt02(UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB, WORK, RESID)
DTPT02
Definition: dtpt02.f:143
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
subroutine dtpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
DTPMV
Definition: dtpmv.f:144