1*> \brief \b DPBT01
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 DPBT01( 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   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> DPBT01 reconstructs a symmetric positive definite band matrix A from
30*> its L*L' or U'*U factorization and computes the residual
31*>    norm( L*L' - A ) / ( N * norm(A) * EPS ) or
32*>    norm( U'*U - A ) / ( N * norm(A) * EPS ),
33*> where EPS is the machine epsilon, L' is the conjugate transpose of
34*> L, and U' is the conjugate transpose of U.
35*> \endverbatim
36*
37*  Arguments:
38*  ==========
39*
40*> \param[in] UPLO
41*> \verbatim
42*>          UPLO is CHARACTER*1
43*>          Specifies whether the upper or lower triangular part of the
44*>          symmetric matrix A is stored:
45*>          = 'U':  Upper triangular
46*>          = 'L':  Lower triangular
47*> \endverbatim
48*>
49*> \param[in] N
50*> \verbatim
51*>          N is INTEGER
52*>          The number of rows and columns of the matrix A.  N >= 0.
53*> \endverbatim
54*>
55*> \param[in] KD
56*> \verbatim
57*>          KD is INTEGER
58*>          The number of super-diagonals of the matrix A if UPLO = 'U',
59*>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
60*> \endverbatim
61*>
62*> \param[in] A
63*> \verbatim
64*>          A is DOUBLE PRECISION array, dimension (LDA,N)
65*>          The original symmetric band matrix A.  If UPLO = 'U', the
66*>          upper triangular part of A is stored as a band matrix; if
67*>          UPLO = 'L', the lower triangular part of A is stored.  The
68*>          columns of the appropriate triangle are stored in the columns
69*>          of A and the diagonals of the triangle are stored in the rows
70*>          of A.  See DPBTRF for further details.
71*> \endverbatim
72*>
73*> \param[in] LDA
74*> \verbatim
75*>          LDA is INTEGER.
76*>          The leading dimension of the array A.  LDA >= max(1,KD+1).
77*> \endverbatim
78*>
79*> \param[in] AFAC
80*> \verbatim
81*>          AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
82*>          The factored form of the matrix A.  AFAC contains the factor
83*>          L or U from the L*L' or U'*U factorization in band storage
84*>          format, as computed by DPBTRF.
85*> \endverbatim
86*>
87*> \param[in] LDAFAC
88*> \verbatim
89*>          LDAFAC is INTEGER
90*>          The leading dimension of the array AFAC.
91*>          LDAFAC >= max(1,KD+1).
92*> \endverbatim
93*>
94*> \param[out] RWORK
95*> \verbatim
96*>          RWORK is DOUBLE PRECISION array, dimension (N)
97*> \endverbatim
98*>
99*> \param[out] RESID
100*> \verbatim
101*>          RESID is DOUBLE PRECISION
102*>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
103*>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
104*> \endverbatim
105*
106*  Authors:
107*  ========
108*
109*> \author Univ. of Tennessee
110*> \author Univ. of California Berkeley
111*> \author Univ. of Colorado Denver
112*> \author NAG Ltd.
113*
114*> \date November 2011
115*
116*> \ingroup double_lin
117*
118*  =====================================================================
119      SUBROUTINE DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
120     $                   RESID )
121*
122*  -- LAPACK test routine (version 3.4.0) --
123*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
124*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125*     November 2011
126*
127*     .. Scalar Arguments ..
128      CHARACTER          UPLO
129      INTEGER            KD, LDA, LDAFAC, N
130      DOUBLE PRECISION   RESID
131*     ..
132*     .. Array Arguments ..
133      DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
134*     ..
135*
136*  =====================================================================
137*
138*
139*     .. Parameters ..
140      DOUBLE PRECISION   ZERO, ONE
141      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
142*     ..
143*     .. Local Scalars ..
144      INTEGER            I, J, K, KC, KLEN, ML, MU
145      DOUBLE PRECISION   ANORM, EPS, T
146*     ..
147*     .. External Functions ..
148      LOGICAL            LSAME
149      DOUBLE PRECISION   DDOT, DLAMCH, DLANSB
150      EXTERNAL           LSAME, DDOT, DLAMCH, DLANSB
151*     ..
152*     .. External Subroutines ..
153      EXTERNAL           DSCAL, DSYR, DTRMV
154*     ..
155*     .. Intrinsic Functions ..
156      INTRINSIC          DBLE, 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 = DLANSB( '1', UPLO, N, KD, A, LDA, RWORK )
171      IF( ANORM.LE.ZERO ) THEN
172         RESID = ONE / EPS
173         RETURN
174      END IF
175*
176*     Compute the product U'*U, overwriting U.
177*
178      IF( LSAME( UPLO, 'U' ) ) THEN
179         DO 10 K = N, 1, -1
180            KC = MAX( 1, KD+2-K )
181            KLEN = KD + 1 - KC
182*
183*           Compute the (K,K) element of the result.
184*
185            T = DDOT( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
186            AFAC( KD+1, K ) = T
187*
188*           Compute the rest of column K.
189*
190            IF( KLEN.GT.0 )
191     $         CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', KLEN,
192     $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
193     $                     AFAC( KC, K ), 1 )
194*
195   10    CONTINUE
196*
197*     UPLO = 'L':  Compute the product L*L', overwriting L.
198*
199      ELSE
200         DO 20 K = N, 1, -1
201            KLEN = MIN( KD, N-K )
202*
203*           Add a multiple of column K of the factor L to each of
204*           columns K+1 through N.
205*
206            IF( KLEN.GT.0 )
207     $         CALL DSYR( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
208     $                    AFAC( 1, K+1 ), LDAFAC-1 )
209*
210*           Scale column K by the diagonal element.
211*
212            T = AFAC( 1, K )
213            CALL DSCAL( KLEN+1, T, AFAC( 1, K ), 1 )
214*
215   20    CONTINUE
216      END IF
217*
218*     Compute the difference  L*L' - A  or  U'*U - A.
219*
220      IF( LSAME( UPLO, 'U' ) ) THEN
221         DO 40 J = 1, N
222            MU = MAX( 1, KD+2-J )
223            DO 30 I = MU, KD + 1
224               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
225   30       CONTINUE
226   40    CONTINUE
227      ELSE
228         DO 60 J = 1, N
229            ML = MIN( KD+1, N-J+1 )
230            DO 50 I = 1, ML
231               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
232   50       CONTINUE
233   60    CONTINUE
234      END IF
235*
236*     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
237*
238      RESID = DLANSB( 'I', UPLO, N, KD, AFAC, LDAFAC, RWORK )
239*
240      RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
241*
242      RETURN
243*
244*     End of DPBT01
245*
246      END
247