1*> \brief \b SPST01
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 SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
12*                          PIV, RWORK, RESID, RANK )
13*
14*       .. Scalar Arguments ..
15*       REAL               RESID
16*       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
17*       CHARACTER          UPLO
18*       ..
19*       .. Array Arguments ..
20*       REAL               A( LDA, * ), AFAC( LDAFAC, * ),
21*      $                   PERM( LDPERM, * ), RWORK( * )
22*       INTEGER            PIV( * )
23*       ..
24*
25*
26*> \par Purpose:
27*  =============
28*>
29*> \verbatim
30*>
31*> SPST01 reconstructs a symmetric positive semidefinite matrix A
32*> from its L or U factors and the permutation matrix P and computes
33*> the residual
34*>    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
35*>    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
36*> where EPS is the machine epsilon.
37*> \endverbatim
38*
39*  Arguments:
40*  ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*>          UPLO is CHARACTER*1
45*>          Specifies whether the upper or lower triangular part of the
46*>          symmetric matrix A is stored:
47*>          = 'U':  Upper triangular
48*>          = 'L':  Lower triangular
49*> \endverbatim
50*>
51*> \param[in] N
52*> \verbatim
53*>          N is INTEGER
54*>          The number of rows and columns of the matrix A.  N >= 0.
55*> \endverbatim
56*>
57*> \param[in] A
58*> \verbatim
59*>          A is REAL array, dimension (LDA,N)
60*>          The original symmetric matrix A.
61*> \endverbatim
62*>
63*> \param[in] LDA
64*> \verbatim
65*>          LDA is INTEGER
66*>          The leading dimension of the array A.  LDA >= max(1,N)
67*> \endverbatim
68*>
69*> \param[in] AFAC
70*> \verbatim
71*>          AFAC is REAL array, dimension (LDAFAC,N)
72*>          The factor L or U from the L*L' or U'*U
73*>          factorization of A.
74*> \endverbatim
75*>
76*> \param[in] LDAFAC
77*> \verbatim
78*>          LDAFAC is INTEGER
79*>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
80*> \endverbatim
81*>
82*> \param[out] PERM
83*> \verbatim
84*>          PERM is REAL array, dimension (LDPERM,N)
85*>          Overwritten with the reconstructed matrix, and then with the
86*>          difference P*L*L'*P' - A (or P*U'*U*P' - A)
87*> \endverbatim
88*>
89*> \param[in] LDPERM
90*> \verbatim
91*>          LDPERM is INTEGER
92*>          The leading dimension of the array PERM.
93*>          LDAPERM >= max(1,N).
94*> \endverbatim
95*>
96*> \param[in] PIV
97*> \verbatim
98*>          PIV is INTEGER array, dimension (N)
99*>          PIV is such that the nonzero entries are
100*>          P( PIV( K ), K ) = 1.
101*> \endverbatim
102*>
103*> \param[out] RWORK
104*> \verbatim
105*>          RWORK is REAL array, dimension (N)
106*> \endverbatim
107*>
108*> \param[out] RESID
109*> \verbatim
110*>          RESID is REAL
111*>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
112*>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
113*> \endverbatim
114*>
115*> \param[in] RANK
116*> \verbatim
117*>          RANK is INTEGER
118*>          number of nonzero singular values of A.
119*> \endverbatim
120*
121*  Authors:
122*  ========
123*
124*> \author Univ. of Tennessee
125*> \author Univ. of California Berkeley
126*> \author Univ. of Colorado Denver
127*> \author NAG Ltd.
128*
129*> \date December 2016
130*
131*> \ingroup single_lin
132*
133*  =====================================================================
134      SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
135     $                   PIV, RWORK, RESID, RANK )
136*
137*  -- LAPACK test routine (version 3.7.0) --
138*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
139*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*     December 2016
141*
142*     .. Scalar Arguments ..
143      REAL               RESID
144      INTEGER            LDA, LDAFAC, LDPERM, N, RANK
145      CHARACTER          UPLO
146*     ..
147*     .. Array Arguments ..
148      REAL               A( LDA, * ), AFAC( LDAFAC, * ),
149     $                   PERM( LDPERM, * ), RWORK( * )
150      INTEGER            PIV( * )
151*     ..
152*
153*  =====================================================================
154*
155*     .. Parameters ..
156      REAL               ZERO, ONE
157      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
158*     ..
159*     .. Local Scalars ..
160      REAL               ANORM, EPS, T
161      INTEGER            I, J, K
162*     ..
163*     .. External Functions ..
164      REAL               SDOT, SLAMCH, SLANSY
165      LOGICAL            LSAME
166      EXTERNAL           SDOT, SLAMCH, SLANSY, LSAME
167*     ..
168*     .. External Subroutines ..
169      EXTERNAL           SSCAL, SSYR, STRMV
170*     ..
171*     .. Intrinsic Functions ..
172      INTRINSIC          REAL
173*     ..
174*     .. Executable Statements ..
175*
176*     Quick exit if N = 0.
177*
178      IF( N.LE.0 ) THEN
179         RESID = ZERO
180         RETURN
181      END IF
182*
183*     Exit with RESID = 1/EPS if ANORM = 0.
184*
185      EPS = SLAMCH( 'Epsilon' )
186      ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
187      IF( ANORM.LE.ZERO ) THEN
188         RESID = ONE / EPS
189         RETURN
190      END IF
191*
192*     Compute the product U'*U, overwriting U.
193*
194      IF( LSAME( UPLO, 'U' ) ) THEN
195*
196         IF( RANK.LT.N ) THEN
197            DO 110 J = RANK + 1, N
198               DO 100 I = RANK + 1, J
199                  AFAC( I, J ) = ZERO
200  100          CONTINUE
201  110       CONTINUE
202         END IF
203*
204         DO 120 K = N, 1, -1
205*
206*           Compute the (K,K) element of the result.
207*
208            T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
209            AFAC( K, K ) = T
210*
211*           Compute the rest of column K.
212*
213            CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
214     $                  LDAFAC, AFAC( 1, K ), 1 )
215*
216  120    CONTINUE
217*
218*     Compute the product L*L', overwriting L.
219*
220      ELSE
221*
222         IF( RANK.LT.N ) THEN
223            DO 140 J = RANK + 1, N
224               DO 130 I = J, N
225                  AFAC( I, J ) = ZERO
226  130          CONTINUE
227  140       CONTINUE
228         END IF
229*
230         DO 150 K = N, 1, -1
231*           Add a multiple of column K of the factor L to each of
232*           columns K+1 through N.
233*
234            IF( K+1.LE.N )
235     $         CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
236     $                    AFAC( K+1, K+1 ), LDAFAC )
237*
238*           Scale column K by the diagonal element.
239*
240            T = AFAC( K, K )
241            CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
242  150    CONTINUE
243*
244      END IF
245*
246*        Form P*L*L'*P' or P*U'*U*P'
247*
248      IF( LSAME( UPLO, 'U' ) ) THEN
249*
250         DO 170 J = 1, N
251            DO 160 I = 1, N
252               IF( PIV( I ).LE.PIV( J ) ) THEN
253                  IF( I.LE.J ) THEN
254                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
255                  ELSE
256                     PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
257                  END IF
258               END IF
259  160       CONTINUE
260  170    CONTINUE
261*
262*
263      ELSE
264*
265         DO 190 J = 1, N
266            DO 180 I = 1, N
267               IF( PIV( I ).GE.PIV( J ) ) THEN
268                  IF( I.GE.J ) THEN
269                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
270                  ELSE
271                     PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
272                  END IF
273               END IF
274  180       CONTINUE
275  190    CONTINUE
276*
277      END IF
278*
279*     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
280*
281      IF( LSAME( UPLO, 'U' ) ) THEN
282         DO 210 J = 1, N
283            DO 200 I = 1, J
284               PERM( I, J ) = PERM( I, J ) - A( I, J )
285  200       CONTINUE
286  210    CONTINUE
287      ELSE
288         DO 230 J = 1, N
289            DO 220 I = J, N
290               PERM( I, J ) = PERM( I, J ) - A( I, J )
291  220       CONTINUE
292  230    CONTINUE
293      END IF
294*
295*     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
296*     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
297*
298      RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
299*
300      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
301*
302      RETURN
303*
304*     End of SPST01
305*
306      END
307