LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zppt03.f
Go to the documentation of this file.
1*> \brief \b ZPPT03
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 ZPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
12* RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER LDWORK, N
17* DOUBLE PRECISION RCOND, RESID
18* ..
19* .. Array Arguments ..
20* DOUBLE PRECISION RWORK( * )
21* COMPLEX*16 A( * ), AINV( * ), WORK( LDWORK, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> ZPPT03 computes the residual for a Hermitian packed matrix times its
31*> inverse:
32*> norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
33*> where EPS is the machine epsilon.
34*> \endverbatim
35*
36* Arguments:
37* ==========
38*
39*> \param[in] UPLO
40*> \verbatim
41*> UPLO is CHARACTER*1
42*> Specifies whether the upper or lower triangular part of the
43*> Hermitian matrix A is stored:
44*> = 'U': Upper triangular
45*> = 'L': Lower triangular
46*> \endverbatim
47*>
48*> \param[in] N
49*> \verbatim
50*> N is INTEGER
51*> The number of rows and columns of the matrix A. N >= 0.
52*> \endverbatim
53*>
54*> \param[in] A
55*> \verbatim
56*> A is COMPLEX*16 array, dimension (N*(N+1)/2)
57*> The original Hermitian matrix A, stored as a packed
58*> triangular matrix.
59*> \endverbatim
60*>
61*> \param[in] AINV
62*> \verbatim
63*> AINV is COMPLEX*16 array, dimension (N*(N+1)/2)
64*> The (Hermitian) inverse of the matrix A, stored as a packed
65*> triangular matrix.
66*> \endverbatim
67*>
68*> \param[out] WORK
69*> \verbatim
70*> WORK is COMPLEX*16 array, dimension (LDWORK,N)
71*> \endverbatim
72*>
73*> \param[in] LDWORK
74*> \verbatim
75*> LDWORK is INTEGER
76*> The leading dimension of the array WORK. LDWORK >= max(1,N).
77*> \endverbatim
78*>
79*> \param[out] RWORK
80*> \verbatim
81*> RWORK is DOUBLE PRECISION array, dimension (N)
82*> \endverbatim
83*>
84*> \param[out] RCOND
85*> \verbatim
86*> RCOND is DOUBLE PRECISION
87*> The reciprocal of the condition number of A, computed as
88*> ( 1/norm(A) ) / norm(AINV).
89*> \endverbatim
90*>
91*> \param[out] RESID
92*> \verbatim
93*> RESID is DOUBLE PRECISION
94*> norm(I - A*AINV) / ( 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 zppt03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
109 $ RESID )
110*
111* -- LAPACK test routine --
112* -- LAPACK is a software package provided by Univ. of Tennessee, --
113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114*
115* .. Scalar Arguments ..
116 CHARACTER UPLO
117 INTEGER LDWORK, N
118 DOUBLE PRECISION RCOND, RESID
119* ..
120* .. Array Arguments ..
121 DOUBLE PRECISION RWORK( * )
122 COMPLEX*16 A( * ), AINV( * ), WORK( LDWORK, * )
123* ..
124*
125* =====================================================================
126*
127* .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE
129 parameter( zero = 0.0d+0, one = 1.0d+0 )
130 COMPLEX*16 CZERO, CONE
131 parameter( czero = ( 0.0d+0, 0.0d+0 ),
132 $ cone = ( 1.0d+0, 0.0d+0 ) )
133* ..
134* .. Local Scalars ..
135 INTEGER I, J, JJ
136 DOUBLE PRECISION AINVNM, ANORM, EPS
137* ..
138* .. External Functions ..
139 LOGICAL LSAME
140 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP
141 EXTERNAL lsame, dlamch, zlange, zlanhp
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC dble, dconjg
145* ..
146* .. External Subroutines ..
147 EXTERNAL zcopy, zhpmv
148* ..
149* .. Executable Statements ..
150*
151* Quick exit if N = 0.
152*
153 IF( n.LE.0 ) THEN
154 rcond = one
155 resid = zero
156 RETURN
157 END IF
158*
159* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
160*
161 eps = dlamch( 'Epsilon' )
162 anorm = zlanhp( '1', uplo, n, a, rwork )
163 ainvnm = zlanhp( '1', uplo, n, ainv, rwork )
164 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
165 rcond = zero
166 resid = one / eps
167 RETURN
168 END IF
169 rcond = ( one / anorm ) / ainvnm
170*
171* UPLO = 'U':
172* Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
173* expand it to a full matrix, then multiply by A one column at a
174* time, moving the result one column to the left.
175*
176 IF( lsame( uplo, 'U' ) ) THEN
177*
178* Copy AINV
179*
180 jj = 1
181 DO 20 j = 1, n - 1
182 CALL zcopy( j, ainv( jj ), 1, work( 1, j+1 ), 1 )
183 DO 10 i = 1, j - 1
184 work( j, i+1 ) = dconjg( ainv( jj+i-1 ) )
185 10 CONTINUE
186 jj = jj + j
187 20 CONTINUE
188 jj = ( ( n-1 )*n ) / 2 + 1
189 DO 30 i = 1, n - 1
190 work( n, i+1 ) = dconjg( ainv( jj+i-1 ) )
191 30 CONTINUE
192*
193* Multiply by A
194*
195 DO 40 j = 1, n - 1
196 CALL zhpmv( 'Upper', n, -cone, a, work( 1, j+1 ), 1, czero,
197 $ work( 1, j ), 1 )
198 40 CONTINUE
199 CALL zhpmv( 'Upper', n, -cone, a, ainv( jj ), 1, czero,
200 $ work( 1, n ), 1 )
201*
202* UPLO = 'L':
203* Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
204* and multiply by A, moving each column to the right.
205*
206 ELSE
207*
208* Copy AINV
209*
210 DO 50 i = 1, n - 1
211 work( 1, i ) = dconjg( ainv( i+1 ) )
212 50 CONTINUE
213 jj = n + 1
214 DO 70 j = 2, n
215 CALL zcopy( n-j+1, ainv( jj ), 1, work( j, j-1 ), 1 )
216 DO 60 i = 1, n - j
217 work( j, j+i-1 ) = dconjg( ainv( jj+i ) )
218 60 CONTINUE
219 jj = jj + n - j + 1
220 70 CONTINUE
221*
222* Multiply by A
223*
224 DO 80 j = n, 2, -1
225 CALL zhpmv( 'Lower', n, -cone, a, work( 1, j-1 ), 1, czero,
226 $ work( 1, j ), 1 )
227 80 CONTINUE
228 CALL zhpmv( 'Lower', n, -cone, a, ainv( 1 ), 1, czero,
229 $ work( 1, 1 ), 1 )
230*
231 END IF
232*
233* Add the identity matrix to WORK .
234*
235 DO 90 i = 1, n
236 work( i, i ) = work( i, i ) + cone
237 90 CONTINUE
238*
239* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
240*
241 resid = zlange( '1', n, n, work, ldwork, rwork )
242*
243 resid = ( ( resid*rcond ) / eps ) / dble( n )
244*
245 RETURN
246*
247* End of ZPPT03
248*
249 END
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zhpmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
ZHPMV
Definition zhpmv.f:149
subroutine zppt03(uplo, n, a, ainv, work, ldwork, rwork, rcond, resid)
ZPPT03
Definition zppt03.f:110