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*> \ingroup single_lin
130*
131*  =====================================================================
132      SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
133     $                   PIV, RWORK, RESID, RANK )
134*
135*  -- LAPACK test routine --
136*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
137*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139*     .. Scalar Arguments ..
140      REAL               RESID
141      INTEGER            LDA, LDAFAC, LDPERM, N, RANK
142      CHARACTER          UPLO
143*     ..
144*     .. Array Arguments ..
145      REAL               A( LDA, * ), AFAC( LDAFAC, * ),
146     $                   PERM( LDPERM, * ), RWORK( * )
147      INTEGER            PIV( * )
148*     ..
149*
150*  =====================================================================
151*
152*     .. Parameters ..
153      REAL               ZERO, ONE
154      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
155*     ..
156*     .. Local Scalars ..
157      REAL               ANORM, EPS, T
158      INTEGER            I, J, K
159*     ..
160*     .. External Functions ..
161      REAL               SDOT, SLAMCH, SLANSY
162      LOGICAL            LSAME
163      EXTERNAL           SDOT, SLAMCH, SLANSY, LSAME
164*     ..
165*     .. External Subroutines ..
166      EXTERNAL           SSCAL, SSYR, STRMV
167*     ..
168*     .. Intrinsic Functions ..
169      INTRINSIC          REAL
170*     ..
171*     .. Executable Statements ..
172*
173*     Quick exit if N = 0.
174*
175      IF( N.LE.0 ) THEN
176         RESID = ZERO
177         RETURN
178      END IF
179*
180*     Exit with RESID = 1/EPS if ANORM = 0.
181*
182      EPS = SLAMCH( 'Epsilon' )
183      ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
184      IF( ANORM.LE.ZERO ) THEN
185         RESID = ONE / EPS
186         RETURN
187      END IF
188*
189*     Compute the product U'*U, overwriting U.
190*
191      IF( LSAME( UPLO, 'U' ) ) THEN
192*
193         IF( RANK.LT.N ) THEN
194            DO 110 J = RANK + 1, N
195               DO 100 I = RANK + 1, J
196                  AFAC( I, J ) = ZERO
197  100          CONTINUE
198  110       CONTINUE
199         END IF
200*
201         DO 120 K = N, 1, -1
202*
203*           Compute the (K,K) element of the result.
204*
205            T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
206            AFAC( K, K ) = T
207*
208*           Compute the rest of column K.
209*
210            CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
211     $                  LDAFAC, AFAC( 1, K ), 1 )
212*
213  120    CONTINUE
214*
215*     Compute the product L*L', overwriting L.
216*
217      ELSE
218*
219         IF( RANK.LT.N ) THEN
220            DO 140 J = RANK + 1, N
221               DO 130 I = J, N
222                  AFAC( I, J ) = ZERO
223  130          CONTINUE
224  140       CONTINUE
225         END IF
226*
227         DO 150 K = N, 1, -1
228*           Add a multiple of column K of the factor L to each of
229*           columns K+1 through N.
230*
231            IF( K+1.LE.N )
232     $         CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
233     $                    AFAC( K+1, K+1 ), LDAFAC )
234*
235*           Scale column K by the diagonal element.
236*
237            T = AFAC( K, K )
238            CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
239  150    CONTINUE
240*
241      END IF
242*
243*        Form P*L*L'*P' or P*U'*U*P'
244*
245      IF( LSAME( UPLO, 'U' ) ) THEN
246*
247         DO 170 J = 1, N
248            DO 160 I = 1, N
249               IF( PIV( I ).LE.PIV( J ) ) THEN
250                  IF( I.LE.J ) THEN
251                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
252                  ELSE
253                     PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
254                  END IF
255               END IF
256  160       CONTINUE
257  170    CONTINUE
258*
259*
260      ELSE
261*
262         DO 190 J = 1, N
263            DO 180 I = 1, N
264               IF( PIV( I ).GE.PIV( J ) ) THEN
265                  IF( I.GE.J ) THEN
266                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
267                  ELSE
268                     PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
269                  END IF
270               END IF
271  180       CONTINUE
272  190    CONTINUE
273*
274      END IF
275*
276*     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
277*
278      IF( LSAME( UPLO, 'U' ) ) THEN
279         DO 210 J = 1, N
280            DO 200 I = 1, J
281               PERM( I, J ) = PERM( I, J ) - A( I, J )
282  200       CONTINUE
283  210    CONTINUE
284      ELSE
285         DO 230 J = 1, N
286            DO 220 I = J, N
287               PERM( I, J ) = PERM( I, J ) - A( I, J )
288  220       CONTINUE
289  230    CONTINUE
290      END IF
291*
292*     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
293*     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
294*
295      RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
296*
297      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
298*
299      RETURN
300*
301*     End of SPST01
302*
303      END
304