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