LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtpt01.f
Go to the documentation of this file.
1*> \brief \b DTPT01
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 DTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, WORK, RESID )
12*
13* .. Scalar Arguments ..
14* CHARACTER DIAG, UPLO
15* INTEGER N
16* DOUBLE PRECISION RCOND, RESID
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION AINVP( * ), AP( * ), WORK( * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> DTPT01 computes the residual for a triangular matrix A times its
29*> inverse when A is stored in packed format:
30*> RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
31*> where EPS is the machine epsilon.
32*> \endverbatim
33*
34* Arguments:
35* ==========
36*
37*> \param[in] UPLO
38*> \verbatim
39*> UPLO is CHARACTER*1
40*> Specifies whether the matrix A is upper or lower triangular.
41*> = 'U': Upper triangular
42*> = 'L': Lower triangular
43*> \endverbatim
44*>
45*> \param[in] DIAG
46*> \verbatim
47*> DIAG is CHARACTER*1
48*> Specifies whether or not the matrix A is unit triangular.
49*> = 'N': Non-unit triangular
50*> = 'U': Unit triangular
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The order of the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in] AP
60*> \verbatim
61*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
62*> The original upper or lower triangular matrix A, packed
63*> columnwise in a linear array. The j-th column of A is stored
64*> in the array AP as follows:
65*> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
66*> if UPLO = 'L',
67*> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
68*> \endverbatim
69*>
70*> \param[in,out] AINVP
71*> \verbatim
72*> AINVP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
73*> On entry, the (triangular) inverse of the matrix A, packed
74*> columnwise in a linear array as in AP.
75*> On exit, the contents of AINVP are destroyed.
76*> \endverbatim
77*>
78*> \param[out] RCOND
79*> \verbatim
80*> RCOND is DOUBLE PRECISION
81*> The reciprocal condition number of A, computed as
82*> 1/(norm(A) * norm(AINV)).
83*> \endverbatim
84*>
85*> \param[out] WORK
86*> \verbatim
87*> WORK is DOUBLE PRECISION array, dimension (N)
88*> \endverbatim
89*>
90*> \param[out] RESID
91*> \verbatim
92*> RESID is DOUBLE PRECISION
93*> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
94*> \endverbatim
95*
96* Authors:
97* ========
98*
99*> \author Univ. of Tennessee
100*> \author Univ. of California Berkeley
101*> \author Univ. of Colorado Denver
102*> \author NAG Ltd.
103*
104*> \ingroup double_lin
105*
106* =====================================================================
107 SUBROUTINE dtpt01( UPLO, DIAG, N, AP, AINVP, RCOND, WORK, RESID )
108*
109* -- LAPACK test routine --
110* -- LAPACK is a software package provided by Univ. of Tennessee, --
111* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
112*
113* .. Scalar Arguments ..
114 CHARACTER DIAG, UPLO
115 INTEGER N
116 DOUBLE PRECISION RCOND, RESID
117* ..
118* .. Array Arguments ..
119 DOUBLE PRECISION AINVP( * ), AP( * ), WORK( * )
120* ..
121*
122* =====================================================================
123*
124* .. Parameters ..
125 DOUBLE PRECISION ZERO, ONE
126 parameter( zero = 0.0d+0, one = 1.0d+0 )
127* ..
128* .. Local Scalars ..
129 LOGICAL UNITD
130 INTEGER J, JC
131 DOUBLE PRECISION AINVNM, ANORM, EPS
132* ..
133* .. External Functions ..
134 LOGICAL LSAME
135 DOUBLE PRECISION DLAMCH, DLANTP
136 EXTERNAL lsame, dlamch, dlantp
137* ..
138* .. External Subroutines ..
139 EXTERNAL dtpmv
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC dble
143* ..
144* .. Executable Statements ..
145*
146* Quick exit if N = 0.
147*
148 IF( n.LE.0 ) THEN
149 rcond = one
150 resid = zero
151 RETURN
152 END IF
153*
154* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
155*
156 eps = dlamch( 'Epsilon' )
157 anorm = dlantp( '1', uplo, diag, n, ap, work )
158 ainvnm = dlantp( '1', uplo, diag, n, ainvp, work )
159 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
160 rcond = zero
161 resid = one / eps
162 RETURN
163 END IF
164 rcond = ( one / anorm ) / ainvnm
165*
166* Compute A * AINV, overwriting AINV.
167*
168 unitd = lsame( diag, 'U' )
169 IF( lsame( uplo, 'U' ) ) THEN
170 jc = 1
171 DO 10 j = 1, n
172 IF( unitd )
173 $ ainvp( jc+j-1 ) = one
174*
175* Form the j-th column of A*AINV
176*
177 CALL dtpmv( 'Upper', 'No transpose', diag, j, ap,
178 $ ainvp( jc ), 1 )
179*
180* Subtract 1 from the diagonal
181*
182 ainvp( jc+j-1 ) = ainvp( jc+j-1 ) - one
183 jc = jc + j
184 10 CONTINUE
185 ELSE
186 jc = 1
187 DO 20 j = 1, n
188 IF( unitd )
189 $ ainvp( jc ) = one
190*
191* Form the j-th column of A*AINV
192*
193 CALL dtpmv( 'Lower', 'No transpose', diag, n-j+1, ap( jc ),
194 $ ainvp( jc ), 1 )
195*
196* Subtract 1 from the diagonal
197*
198 ainvp( jc ) = ainvp( jc ) - one
199 jc = jc + n - j + 1
200 20 CONTINUE
201 END IF
202*
203* Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
204*
205 resid = dlantp( '1', uplo, 'Non-unit', n, ainvp, work )
206*
207 resid = ( ( resid*rcond ) / dble( n ) ) / eps
208*
209 RETURN
210*
211* End of DTPT01
212*
213 END
subroutine dtpt01(uplo, diag, n, ap, ainvp, rcond, work, resid)
DTPT01
Definition dtpt01.f:108
subroutine dtpmv(uplo, trans, diag, n, ap, x, incx)
DTPMV
Definition dtpmv.f:142