LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
116 *
117 *> \ingroup complex16_lin
118 *
119 * =====================================================================
120  SUBROUTINE zpbt01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
121  $ resid )
122 *
123 * -- LAPACK test routine (version 3.4.0) --
124 * -- LAPACK is a software package provided by Univ. of Tennessee, --
125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126 * November 2011
127 *
128 * .. Scalar Arguments ..
129  CHARACTER uplo
130  INTEGER kd, lda, ldafac, n
131  DOUBLE PRECISION resid
132 * ..
133 * .. Array Arguments ..
134  DOUBLE PRECISION rwork( * )
135  COMPLEX*16 a( lda, * ), afac( ldafac, * )
136 * ..
137 *
138 * =====================================================================
139 *
140 *
141 * .. Parameters ..
142  DOUBLE PRECISION zero, one
143  parameter( zero = 0.0d+0, one = 1.0d+0 )
144 * ..
145 * .. Local Scalars ..
146  INTEGER i, j, k, kc, klen, ml, mu
147  DOUBLE PRECISION akk, anorm, eps
148 * ..
149 * .. External Functions ..
150  LOGICAL lsame
151  DOUBLE PRECISION dlamch, zlanhb
152  COMPLEX*16 zdotc
153  EXTERNAL lsame, dlamch, zlanhb, zdotc
154 * ..
155 * .. External Subroutines ..
156  EXTERNAL zdscal, zher, ztrmv
157 * ..
158 * .. Intrinsic Functions ..
159  INTRINSIC dble, dimag, max, min
160 * ..
161 * .. Executable Statements ..
162 *
163 * Quick exit if N = 0.
164 *
165  IF( n.LE.0 ) THEN
166  resid = zero
167  return
168  END IF
169 *
170 * Exit with RESID = 1/EPS if ANORM = 0.
171 *
172  eps = dlamch( 'Epsilon' )
173  anorm = zlanhb( '1', uplo, n, kd, a, lda, rwork )
174  IF( anorm.LE.zero ) THEN
175  resid = one / eps
176  return
177  END IF
178 *
179 * Check the imaginary parts of the diagonal elements and return with
180 * an error code if any are nonzero.
181 *
182  IF( lsame( uplo, 'U' ) ) THEN
183  DO 10 j = 1, n
184  IF( dimag( afac( kd+1, j ) ).NE.zero ) THEN
185  resid = one / eps
186  return
187  END IF
188  10 continue
189  ELSE
190  DO 20 j = 1, n
191  IF( dimag( afac( 1, j ) ).NE.zero ) THEN
192  resid = one / eps
193  return
194  END IF
195  20 continue
196  END IF
197 *
198 * Compute the product U'*U, overwriting U.
199 *
200  IF( lsame( uplo, 'U' ) ) THEN
201  DO 30 k = n, 1, -1
202  kc = max( 1, kd+2-k )
203  klen = kd + 1 - kc
204 *
205 * Compute the (K,K) element of the result.
206 *
207  akk = zdotc( klen+1, afac( kc, k ), 1, afac( kc, k ), 1 )
208  afac( kd+1, k ) = akk
209 *
210 * Compute the rest of column K.
211 *
212  IF( klen.GT.0 )
213  $ CALL ztrmv( 'Upper', 'Conjugate', 'Non-unit', klen,
214  $ afac( kd+1, k-klen ), ldafac-1,
215  $ afac( kc, k ), 1 )
216 *
217  30 continue
218 *
219 * UPLO = 'L': Compute the product L*L', overwriting L.
220 *
221  ELSE
222  DO 40 k = n, 1, -1
223  klen = min( kd, n-k )
224 *
225 * Add a multiple of column K of the factor L to each of
226 * columns K+1 through N.
227 *
228  IF( klen.GT.0 )
229  $ CALL zher( 'Lower', klen, one, afac( 2, k ), 1,
230  $ afac( 1, k+1 ), ldafac-1 )
231 *
232 * Scale column K by the diagonal element.
233 *
234  akk = afac( 1, k )
235  CALL zdscal( klen+1, akk, afac( 1, k ), 1 )
236 *
237  40 continue
238  END IF
239 *
240 * Compute the difference L*L' - A or U'*U - A.
241 *
242  IF( lsame( uplo, 'U' ) ) THEN
243  DO 60 j = 1, n
244  mu = max( 1, kd+2-j )
245  DO 50 i = mu, kd + 1
246  afac( i, j ) = afac( i, j ) - a( i, j )
247  50 continue
248  60 continue
249  ELSE
250  DO 80 j = 1, n
251  ml = min( kd+1, n-j+1 )
252  DO 70 i = 1, ml
253  afac( i, j ) = afac( i, j ) - a( i, j )
254  70 continue
255  80 continue
256  END IF
257 *
258 * Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
259 *
260  resid = zlanhb( '1', uplo, n, kd, afac, ldafac, rwork )
261 *
262  resid = ( ( resid / dble( n ) ) / anorm ) / eps
263 *
264  return
265 *
266 * End of ZPBT01
267 *
268  END