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 = ZDOTC( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
205            AFAC( KD+1, K ) = AKK
206*
207*           Compute the rest of column K.
208*
209            IF( KLEN.GT.0 )
210     $         CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', KLEN,
211     $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
212     $                     AFAC( KC, K ), 1 )
213*
214   30    CONTINUE
215*
216*     UPLO = 'L':  Compute the product L*L', overwriting L.
217*
218      ELSE
219         DO 40 K = N, 1, -1
220            KLEN = MIN( KD, N-K )
221*
222*           Add a multiple of column K of the factor L to each of
223*           columns K+1 through N.
224*
225            IF( KLEN.GT.0 )
226     $         CALL ZHER( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
227     $                    AFAC( 1, K+1 ), LDAFAC-1 )
228*
229*           Scale column K by the diagonal element.
230*
231            AKK = AFAC( 1, K )
232            CALL ZDSCAL( KLEN+1, AKK, AFAC( 1, K ), 1 )
233*
234   40    CONTINUE
235      END IF
236*
237*     Compute the difference  L*L' - A  or  U'*U - A.
238*
239      IF( LSAME( UPLO, 'U' ) ) THEN
240         DO 60 J = 1, N
241            MU = MAX( 1, KD+2-J )
242            DO 50 I = MU, KD + 1
243               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
244   50       CONTINUE
245   60    CONTINUE
246      ELSE
247         DO 80 J = 1, N
248            ML = MIN( KD+1, N-J+1 )
249            DO 70 I = 1, ML
250               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
251   70       CONTINUE
252   80    CONTINUE
253      END IF
254*
255*     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
256*
257      RESID = ZLANHB( '1', UPLO, N, KD, AFAC, LDAFAC, RWORK )
258*
259      RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
260*
261      RETURN
262*
263*     End of ZPBT01
264*
265      END
266