1*> \brief \b SPPT01
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 SPPT01( UPLO, N, A, AFAC, RWORK, RESID )
12*
13*       .. Scalar Arguments ..
14*       CHARACTER          UPLO
15*       INTEGER            N
16*       REAL               RESID
17*       ..
18*       .. Array Arguments ..
19*       REAL               A( * ), AFAC( * ), RWORK( * )
20*       ..
21*
22*
23*> \par Purpose:
24*  =============
25*>
26*> \verbatim
27*>
28*> SPPT01 reconstructs a symmetric positive definite packed matrix A
29*> from its L*L' or U'*U factorization and computes the residual
30*>    norm( L*L' - A ) / ( N * norm(A) * EPS ) or
31*>    norm( U'*U - A ) / ( N * norm(A) * EPS ),
32*> where 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 REAL 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,out] AFAC
61*> \verbatim
62*>          AFAC is REAL array, dimension (N*(N+1)/2)
63*>          On entry, the factor L or U from the L*L' or U'*U
64*>          factorization of A, stored as a packed triangular matrix.
65*>          Overwritten with the reconstructed matrix, and then with the
66*>          difference L*L' - A (or U'*U - A).
67*> \endverbatim
68*>
69*> \param[out] RWORK
70*> \verbatim
71*>          RWORK is REAL array, dimension (N)
72*> \endverbatim
73*>
74*> \param[out] RESID
75*> \verbatim
76*>          RESID is REAL
77*>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
78*>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
79*> \endverbatim
80*
81*  Authors:
82*  ========
83*
84*> \author Univ. of Tennessee
85*> \author Univ. of California Berkeley
86*> \author Univ. of Colorado Denver
87*> \author NAG Ltd.
88*
89*> \ingroup single_lin
90*
91*  =====================================================================
92      SUBROUTINE SPPT01( UPLO, N, A, AFAC, RWORK, RESID )
93*
94*  -- LAPACK test routine --
95*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
96*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97*
98*     .. Scalar Arguments ..
99      CHARACTER          UPLO
100      INTEGER            N
101      REAL               RESID
102*     ..
103*     .. Array Arguments ..
104      REAL               A( * ), AFAC( * ), RWORK( * )
105*     ..
106*
107*  =====================================================================
108*
109*     .. Parameters ..
110      REAL               ZERO, ONE
111      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
112*     ..
113*     .. Local Scalars ..
114      INTEGER            I, K, KC, NPP
115      REAL               ANORM, EPS, T
116*     ..
117*     .. External Functions ..
118      LOGICAL            LSAME
119      REAL               SDOT, SLAMCH, SLANSP
120      EXTERNAL           LSAME, SDOT, SLAMCH, SLANSP
121*     ..
122*     .. External Subroutines ..
123      EXTERNAL           SSCAL, SSPR, STPMV
124*     ..
125*     .. Intrinsic Functions ..
126      INTRINSIC          REAL
127*     ..
128*     .. Executable Statements ..
129*
130*     Quick exit if N = 0
131*
132      IF( N.LE.0 ) THEN
133         RESID = ZERO
134         RETURN
135      END IF
136*
137*     Exit with RESID = 1/EPS if ANORM = 0.
138*
139      EPS = SLAMCH( 'Epsilon' )
140      ANORM = SLANSP( '1', UPLO, N, A, RWORK )
141      IF( ANORM.LE.ZERO ) THEN
142         RESID = ONE / EPS
143         RETURN
144      END IF
145*
146*     Compute the product U'*U, overwriting U.
147*
148      IF( LSAME( UPLO, 'U' ) ) THEN
149         KC = ( N*( N-1 ) ) / 2 + 1
150         DO 10 K = N, 1, -1
151*
152*           Compute the (K,K) element of the result.
153*
154            T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 )
155            AFAC( KC+K-1 ) = T
156*
157*           Compute the rest of column K.
158*
159            IF( K.GT.1 ) THEN
160               CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
161     $                     AFAC( KC ), 1 )
162               KC = KC - ( K-1 )
163            END IF
164   10    CONTINUE
165*
166*     Compute the product L*L', overwriting L.
167*
168      ELSE
169         KC = ( N*( N+1 ) ) / 2
170         DO 20 K = N, 1, -1
171*
172*           Add a multiple of column K of the factor L to each of
173*           columns K+1 through N.
174*
175            IF( K.LT.N )
176     $         CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
177     $                    AFAC( KC+N-K+1 ) )
178*
179*           Scale column K by the diagonal element.
180*
181            T = AFAC( KC )
182            CALL SSCAL( N-K+1, T, AFAC( KC ), 1 )
183*
184            KC = KC - ( N-K+2 )
185   20    CONTINUE
186      END IF
187*
188*     Compute the difference  L*L' - A (or U'*U - A).
189*
190      NPP = N*( N+1 ) / 2
191      DO 30 I = 1, NPP
192         AFAC( I ) = AFAC( I ) - A( I )
193   30 CONTINUE
194*
195*     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
196*
197      RESID = SLANSP( '1', UPLO, N, AFAC, RWORK )
198*
199      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
200*
201      RETURN
202*
203*     End of SPPT01
204*
205      END
206