LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
spst01.f
Go to the documentation of this file.
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
subroutine ssyr(uplo, n, alpha, x, incx, a, lda)
SSYR
Definition ssyr.f:132
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine strmv(uplo, trans, diag, n, a, lda, x, incx)
STRMV
Definition strmv.f:147
subroutine spst01(uplo, n, a, lda, afac, ldafac, perm, ldperm, piv, rwork, resid, rank)
SPST01
Definition spst01.f:134