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*> \ingroup double_lin
115*
116*  =====================================================================
117      SUBROUTINE DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
118     $                   RESID )
119*
120*  -- LAPACK test routine --
121*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
122*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123*
124*     .. Scalar Arguments ..
125      CHARACTER          UPLO
126      INTEGER            KD, LDA, LDAFAC, N
127      DOUBLE PRECISION   RESID
128*     ..
129*     .. Array Arguments ..
130      DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
131*     ..
132*
133*  =====================================================================
134*
135*
136*     .. Parameters ..
137      DOUBLE PRECISION   ZERO, ONE
138      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
139*     ..
140*     .. Local Scalars ..
141      INTEGER            I, J, K, KC, KLEN, ML, MU
142      DOUBLE PRECISION   ANORM, EPS, T
143*     ..
144*     .. External Functions ..
145      LOGICAL            LSAME
146      DOUBLE PRECISION   DDOT, DLAMCH, DLANSB
147      EXTERNAL           LSAME, DDOT, DLAMCH, DLANSB
148*     ..
149*     .. External Subroutines ..
150      EXTERNAL           DSCAL, DSYR, DTRMV
151*     ..
152*     .. Intrinsic Functions ..
153      INTRINSIC          DBLE, MAX, MIN
154*     ..
155*     .. Executable Statements ..
156*
157*     Quick exit if N = 0.
158*
159      IF( N.LE.0 ) THEN
160         RESID = ZERO
161         RETURN
162      END IF
163*
164*     Exit with RESID = 1/EPS if ANORM = 0.
165*
166      EPS = DLAMCH( 'Epsilon' )
167      ANORM = DLANSB( '1', UPLO, N, KD, A, LDA, RWORK )
168      IF( ANORM.LE.ZERO ) THEN
169         RESID = ONE / EPS
170         RETURN
171      END IF
172*
173*     Compute the product U'*U, overwriting U.
174*
175      IF( LSAME( UPLO, 'U' ) ) THEN
176         DO 10 K = N, 1, -1
177            KC = MAX( 1, KD+2-K )
178            KLEN = KD + 1 - KC
179*
180*           Compute the (K,K) element of the result.
181*
182            T = DDOT( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
183            AFAC( KD+1, K ) = T
184*
185*           Compute the rest of column K.
186*
187            IF( KLEN.GT.0 )
188     $         CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', KLEN,
189     $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
190     $                     AFAC( KC, K ), 1 )
191*
192   10    CONTINUE
193*
194*     UPLO = 'L':  Compute the product L*L', overwriting L.
195*
196      ELSE
197         DO 20 K = N, 1, -1
198            KLEN = MIN( KD, N-K )
199*
200*           Add a multiple of column K of the factor L to each of
201*           columns K+1 through N.
202*
203            IF( KLEN.GT.0 )
204     $         CALL DSYR( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
205     $                    AFAC( 1, K+1 ), LDAFAC-1 )
206*
207*           Scale column K by the diagonal element.
208*
209            T = AFAC( 1, K )
210            CALL DSCAL( KLEN+1, T, AFAC( 1, K ), 1 )
211*
212   20    CONTINUE
213      END IF
214*
215*     Compute the difference  L*L' - A  or  U'*U - A.
216*
217      IF( LSAME( UPLO, 'U' ) ) THEN
218         DO 40 J = 1, N
219            MU = MAX( 1, KD+2-J )
220            DO 30 I = MU, KD + 1
221               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
222   30       CONTINUE
223   40    CONTINUE
224      ELSE
225         DO 60 J = 1, N
226            ML = MIN( KD+1, N-J+1 )
227            DO 50 I = 1, ML
228               AFAC( I, J ) = AFAC( I, J ) - A( I, J )
229   50       CONTINUE
230   60    CONTINUE
231      END IF
232*
233*     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
234*
235      RESID = DLANSB( 'I', UPLO, N, KD, AFAC, LDAFAC, RWORK )
236*
237      RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
238*
239      RETURN
240*
241*     End of DPBT01
242*
243      END
244