LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zpbt01.f
Go to the documentation of this file.
1*> \brief \b ZPBT01
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 ZPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
12* RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER KD, LDA, LDAFAC, N
17* DOUBLE PRECISION RESID
18* ..
19* .. Array Arguments ..
20* DOUBLE PRECISION RWORK( * )
21* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> ZPBT01 reconstructs a Hermitian positive definite band matrix A from
31*> its L*L' or U'*U factorization and computes the residual
32*> norm( L*L' - A ) / ( N * norm(A) * EPS ) or
33*> norm( U'*U - A ) / ( N * norm(A) * EPS ),
34*> where EPS is the machine epsilon, L' is the conjugate transpose of
35*> L, and U' is the conjugate transpose of U.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] UPLO
42*> \verbatim
43*> UPLO is CHARACTER*1
44*> Specifies whether the upper or lower triangular part of the
45*> Hermitian matrix A is stored:
46*> = 'U': Upper triangular
47*> = 'L': Lower triangular
48*> \endverbatim
49*>
50*> \param[in] N
51*> \verbatim
52*> N is INTEGER
53*> The number of rows and columns of the matrix A. N >= 0.
54*> \endverbatim
55*>
56*> \param[in] KD
57*> \verbatim
58*> KD is INTEGER
59*> The number of super-diagonals of the matrix A if UPLO = 'U',
60*> or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
61*> \endverbatim
62*>
63*> \param[in] A
64*> \verbatim
65*> A is COMPLEX*16 array, dimension (LDA,N)
66*> The original Hermitian band matrix A. If UPLO = 'U', the
67*> upper triangular part of A is stored as a band matrix; if
68*> UPLO = 'L', the lower triangular part of A is stored. The
69*> columns of the appropriate triangle are stored in the columns
70*> of A and the diagonals of the triangle are stored in the rows
71*> of A. See ZPBTRF for further details.
72*> \endverbatim
73*>
74*> \param[in] LDA
75*> \verbatim
76*> LDA is INTEGER.
77*> The leading dimension of the array A. LDA >= max(1,KD+1).
78*> \endverbatim
79*>
80*> \param[in] AFAC
81*> \verbatim
82*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
83*> The factored form of the matrix A. AFAC contains the factor
84*> L or U from the L*L' or U'*U factorization in band storage
85*> format, as computed by ZPBTRF.
86*> \endverbatim
87*>
88*> \param[in] LDAFAC
89*> \verbatim
90*> LDAFAC is INTEGER
91*> The leading dimension of the array AFAC.
92*> LDAFAC >= max(1,KD+1).
93*> \endverbatim
94*>
95*> \param[out] RWORK
96*> \verbatim
97*> RWORK is DOUBLE PRECISION array, dimension (N)
98*> \endverbatim
99*>
100*> \param[out] RESID
101*> \verbatim
102*> RESID is DOUBLE PRECISION
103*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
104*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
105*> \endverbatim
106*
107* Authors:
108* ========
109*
110*> \author Univ. of Tennessee
111*> \author Univ. of California Berkeley
112*> \author Univ. of Colorado Denver
113*> \author NAG Ltd.
114*
115*> \ingroup complex16_lin
116*
117* =====================================================================
118 SUBROUTINE zpbt01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
119 $ RESID )
120*
121* -- LAPACK test routine --
122* -- LAPACK is a software package provided by Univ. of Tennessee, --
123* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124*
125* .. Scalar Arguments ..
126 CHARACTER UPLO
127 INTEGER KD, LDA, LDAFAC, N
128 DOUBLE PRECISION RESID
129* ..
130* .. Array Arguments ..
131 DOUBLE PRECISION RWORK( * )
132 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * )
133* ..
134*
135* =====================================================================
136*
137*
138* .. Parameters ..
139 DOUBLE PRECISION ZERO, ONE
140 parameter( zero = 0.0d+0, one = 1.0d+0 )
141* ..
142* .. Local Scalars ..
143 INTEGER I, J, K, KC, KLEN, ML, MU
144 DOUBLE PRECISION AKK, ANORM, EPS
145* ..
146* .. External Functions ..
147 LOGICAL LSAME
148 DOUBLE PRECISION DLAMCH, ZLANHB
149 COMPLEX*16 ZDOTC
150 EXTERNAL lsame, dlamch, zlanhb, zdotc
151* ..
152* .. External Subroutines ..
153 EXTERNAL zdscal, zher, ztrmv
154* ..
155* .. Intrinsic Functions ..
156 INTRINSIC dble, dimag, max, min
157* ..
158* .. Executable Statements ..
159*
160* Quick exit if N = 0.
161*
162 IF( n.LE.0 ) THEN
163 resid = zero
164 RETURN
165 END IF
166*
167* Exit with RESID = 1/EPS if ANORM = 0.
168*
169 eps = dlamch( 'Epsilon' )
170 anorm = zlanhb( '1', uplo, n, kd, a, lda, rwork )
171 IF( anorm.LE.zero ) THEN
172 resid = one / eps
173 RETURN
174 END IF
175*
176* Check the imaginary parts of the diagonal elements and return with
177* an error code if any are nonzero.
178*
179 IF( lsame( uplo, 'U' ) ) THEN
180 DO 10 j = 1, n
181 IF( dimag( afac( kd+1, j ) ).NE.zero ) THEN
182 resid = one / eps
183 RETURN
184 END IF
185 10 CONTINUE
186 ELSE
187 DO 20 j = 1, n
188 IF( dimag( afac( 1, j ) ).NE.zero ) THEN
189 resid = one / eps
190 RETURN
191 END IF
192 20 CONTINUE
193 END IF
194*
195* Compute the product U'*U, overwriting U.
196*
197 IF( lsame( uplo, 'U' ) ) THEN
198 DO 30 k = n, 1, -1
199 kc = max( 1, kd+2-k )
200 klen = kd + 1 - kc
201*
202* Compute the (K,K) element of the result.
203*
204 akk = dble(
205 $ zdotc( klen+1, afac( kc, k ), 1, afac( kc, k ), 1 ) )
206 afac( kd+1, k ) = akk
207*
208* Compute the rest of column K.
209*
210 IF( klen.GT.0 )
211 $ CALL ztrmv( 'Upper', 'Conjugate', 'Non-unit', klen,
212 $ afac( kd+1, k-klen ), ldafac-1,
213 $ afac( kc, k ), 1 )
214*
215 30 CONTINUE
216*
217* UPLO = 'L': Compute the product L*L', overwriting L.
218*
219 ELSE
220 DO 40 k = n, 1, -1
221 klen = min( kd, n-k )
222*
223* Add a multiple of column K of the factor L to each of
224* columns K+1 through N.
225*
226 IF( klen.GT.0 )
227 $ CALL zher( 'Lower', klen, one, afac( 2, k ), 1,
228 $ afac( 1, k+1 ), ldafac-1 )
229*
230* Scale column K by the diagonal element.
231*
232 akk = dble( afac( 1, k ) )
233 CALL zdscal( klen+1, akk, afac( 1, k ), 1 )
234*
235 40 CONTINUE
236 END IF
237*
238* Compute the difference L*L' - A or U'*U - A.
239*
240 IF( lsame( uplo, 'U' ) ) THEN
241 DO 60 j = 1, n
242 mu = max( 1, kd+2-j )
243 DO 50 i = mu, kd + 1
244 afac( i, j ) = afac( i, j ) - a( i, j )
245 50 CONTINUE
246 60 CONTINUE
247 ELSE
248 DO 80 j = 1, n
249 ml = min( kd+1, n-j+1 )
250 DO 70 i = 1, ml
251 afac( i, j ) = afac( i, j ) - a( i, j )
252 70 CONTINUE
253 80 CONTINUE
254 END IF
255*
256* Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
257*
258 resid = zlanhb( '1', uplo, n, kd, afac, ldafac, rwork )
259*
260 resid = ( ( resid / dble( n ) ) / anorm ) / eps
261*
262 RETURN
263*
264* End of ZPBT01
265*
266 END
subroutine zher(uplo, n, alpha, x, incx, a, lda)
ZHER
Definition zher.f:135
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine ztrmv(uplo, trans, diag, n, a, lda, x, incx)
ZTRMV
Definition ztrmv.f:147
subroutine zpbt01(uplo, n, kd, a, lda, afac, ldafac, rwork, resid)
ZPBT01
Definition zpbt01.f:120