1*> \brief \b DSPT01
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 DSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
12*
13*       .. Scalar Arguments ..
14*       CHARACTER          UPLO
15*       INTEGER            LDC, N
16*       DOUBLE PRECISION   RESID
17*       ..
18*       .. Array Arguments ..
19*       INTEGER            IPIV( * )
20*       DOUBLE PRECISION   A( * ), AFAC( * ), C( LDC, * ), RWORK( * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> DSPT01 reconstructs a symmetric indefinite packed matrix A from its
30*> block L*D*L' or U*D*U' factorization and computes the residual
31*>      norm( C - A ) / ( N * norm(A) * EPS ),
32*> where C is the reconstructed matrix and EPS is the machine epsilon.
33*> \endverbatim
34*
35*  Arguments:
36*  ==========
37*
38*> \param[in] UPLO
39*> \verbatim
40*>          UPLO is CHARACTER*1
41*>          Specifies whether the upper or lower triangular part of the
42*>          symmetric matrix A is stored:
43*>          = 'U':  Upper triangular
44*>          = 'L':  Lower triangular
45*> \endverbatim
46*>
47*> \param[in] N
48*> \verbatim
49*>          N is INTEGER
50*>          The number of rows and columns of the matrix A.  N >= 0.
51*> \endverbatim
52*>
53*> \param[in] A
54*> \verbatim
55*>          A is DOUBLE PRECISION array, dimension (N*(N+1)/2)
56*>          The original symmetric matrix A, stored as a packed
57*>          triangular matrix.
58*> \endverbatim
59*>
60*> \param[in] AFAC
61*> \verbatim
62*>          AFAC is DOUBLE PRECISION array, dimension (N*(N+1)/2)
63*>          The factored form of the matrix A, stored as a packed
64*>          triangular matrix.  AFAC contains the block diagonal matrix D
65*>          and the multipliers used to obtain the factor L or U from the
66*>          block L*D*L' or U*D*U' factorization as computed by DSPTRF.
67*> \endverbatim
68*>
69*> \param[in] IPIV
70*> \verbatim
71*>          IPIV is INTEGER array, dimension (N)
72*>          The pivot indices from DSPTRF.
73*> \endverbatim
74*>
75*> \param[out] C
76*> \verbatim
77*>          C is DOUBLE PRECISION array, dimension (LDC,N)
78*> \endverbatim
79*>
80*> \param[in] LDC
81*> \verbatim
82*>          LDC is INTEGER
83*>          The leading dimension of the array C.  LDC >= max(1,N).
84*> \endverbatim
85*>
86*> \param[out] RWORK
87*> \verbatim
88*>          RWORK is DOUBLE PRECISION array, dimension (N)
89*> \endverbatim
90*>
91*> \param[out] RESID
92*> \verbatim
93*>          RESID is DOUBLE PRECISION
94*>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
95*>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
96*> \endverbatim
97*
98*  Authors:
99*  ========
100*
101*> \author Univ. of Tennessee
102*> \author Univ. of California Berkeley
103*> \author Univ. of Colorado Denver
104*> \author NAG Ltd.
105*
106*> \date November 2011
107*
108*> \ingroup double_lin
109*
110*  =====================================================================
111      SUBROUTINE DSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
112*
113*  -- LAPACK test routine (version 3.4.0) --
114*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
115*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
116*     November 2011
117*
118*     .. Scalar Arguments ..
119      CHARACTER          UPLO
120      INTEGER            LDC, N
121      DOUBLE PRECISION   RESID
122*     ..
123*     .. Array Arguments ..
124      INTEGER            IPIV( * )
125      DOUBLE PRECISION   A( * ), AFAC( * ), C( LDC, * ), RWORK( * )
126*     ..
127*
128*  =====================================================================
129*
130*     .. Parameters ..
131      DOUBLE PRECISION   ZERO, ONE
132      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
133*     ..
134*     .. Local Scalars ..
135      INTEGER            I, INFO, J, JC
136      DOUBLE PRECISION   ANORM, EPS
137*     ..
138*     .. External Functions ..
139      LOGICAL            LSAME
140      DOUBLE PRECISION   DLAMCH, DLANSP, DLANSY
141      EXTERNAL           LSAME, DLAMCH, DLANSP, DLANSY
142*     ..
143*     .. External Subroutines ..
144      EXTERNAL           DLASET, DLAVSP
145*     ..
146*     .. Intrinsic Functions ..
147      INTRINSIC          DBLE
148*     ..
149*     .. Executable Statements ..
150*
151*     Quick exit if N = 0.
152*
153      IF( N.LE.0 ) THEN
154         RESID = ZERO
155         RETURN
156      END IF
157*
158*     Determine EPS and the norm of A.
159*
160      EPS = DLAMCH( 'Epsilon' )
161      ANORM = DLANSP( '1', UPLO, N, A, RWORK )
162*
163*     Initialize C to the identity matrix.
164*
165      CALL DLASET( 'Full', N, N, ZERO, ONE, C, LDC )
166*
167*     Call DLAVSP to form the product D * U' (or D * L' ).
168*
169      CALL DLAVSP( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, IPIV, C,
170     $             LDC, INFO )
171*
172*     Call DLAVSP again to multiply by U ( or L ).
173*
174      CALL DLAVSP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C,
175     $             LDC, INFO )
176*
177*     Compute the difference  C - A .
178*
179      IF( LSAME( UPLO, 'U' ) ) THEN
180         JC = 0
181         DO 20 J = 1, N
182            DO 10 I = 1, J
183               C( I, J ) = C( I, J ) - A( JC+I )
184   10       CONTINUE
185            JC = JC + J
186   20    CONTINUE
187      ELSE
188         JC = 1
189         DO 40 J = 1, N
190            DO 30 I = J, N
191               C( I, J ) = C( I, J ) - A( JC+I-J )
192   30       CONTINUE
193            JC = JC + N - J + 1
194   40    CONTINUE
195      END IF
196*
197*     Compute norm( C - A ) / ( N * norm(A) * EPS )
198*
199      RESID = DLANSY( '1', UPLO, N, C, LDC, RWORK )
200*
201      IF( ANORM.LE.ZERO ) THEN
202         IF( RESID.NE.ZERO )
203     $      RESID = ONE / EPS
204      ELSE
205         RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
206      END IF
207*
208      RETURN
209*
210*     End of DSPT01
211*
212      END
213