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*> \date November 2011
90*
91*> \ingroup single_lin
92*
93*  =====================================================================
94      SUBROUTINE SPPT01( UPLO, N, A, AFAC, RWORK, RESID )
95*
96*  -- LAPACK test routine (version 3.4.0) --
97*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
98*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
99*     November 2011
100*
101*     .. Scalar Arguments ..
102      CHARACTER          UPLO
103      INTEGER            N
104      REAL               RESID
105*     ..
106*     .. Array Arguments ..
107      REAL               A( * ), AFAC( * ), RWORK( * )
108*     ..
109*
110*  =====================================================================
111*
112*     .. Parameters ..
113      REAL               ZERO, ONE
114      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
115*     ..
116*     .. Local Scalars ..
117      INTEGER            I, K, KC, NPP
118      REAL               ANORM, EPS, T
119*     ..
120*     .. External Functions ..
121      LOGICAL            LSAME
122      REAL               SDOT, SLAMCH, SLANSP
123      EXTERNAL           LSAME, SDOT, SLAMCH, SLANSP
124*     ..
125*     .. External Subroutines ..
126      EXTERNAL           SSCAL, SSPR, STPMV
127*     ..
128*     .. Intrinsic Functions ..
129      INTRINSIC          REAL
130*     ..
131*     .. Executable Statements ..
132*
133*     Quick exit if N = 0
134*
135      IF( N.LE.0 ) THEN
136         RESID = ZERO
137         RETURN
138      END IF
139*
140*     Exit with RESID = 1/EPS if ANORM = 0.
141*
142      EPS = SLAMCH( 'Epsilon' )
143      ANORM = SLANSP( '1', UPLO, N, A, RWORK )
144      IF( ANORM.LE.ZERO ) THEN
145         RESID = ONE / EPS
146         RETURN
147      END IF
148*
149*     Compute the product U'*U, overwriting U.
150*
151      IF( LSAME( UPLO, 'U' ) ) THEN
152         KC = ( N*( N-1 ) ) / 2 + 1
153         DO 10 K = N, 1, -1
154*
155*           Compute the (K,K) element of the result.
156*
157            T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 )
158            AFAC( KC+K-1 ) = T
159*
160*           Compute the rest of column K.
161*
162            IF( K.GT.1 ) THEN
163               CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
164     $                     AFAC( KC ), 1 )
165               KC = KC - ( K-1 )
166            END IF
167   10    CONTINUE
168*
169*     Compute the product L*L', overwriting L.
170*
171      ELSE
172         KC = ( N*( N+1 ) ) / 2
173         DO 20 K = N, 1, -1
174*
175*           Add a multiple of column K of the factor L to each of
176*           columns K+1 through N.
177*
178            IF( K.LT.N )
179     $         CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
180     $                    AFAC( KC+N-K+1 ) )
181*
182*           Scale column K by the diagonal element.
183*
184            T = AFAC( KC )
185            CALL SSCAL( N-K+1, T, AFAC( KC ), 1 )
186*
187            KC = KC - ( N-K+2 )
188   20    CONTINUE
189      END IF
190*
191*     Compute the difference  L*L' - A (or U'*U - A).
192*
193      NPP = N*( N+1 ) / 2
194      DO 30 I = 1, NPP
195         AFAC( I ) = AFAC( I ) - A( I )
196   30 CONTINUE
197*
198*     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
199*
200      RESID = SLANSP( '1', UPLO, N, AFAC, RWORK )
201*
202      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
203*
204      RETURN
205*
206*     End of SPPT01
207*
208      END
209