LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ spst01()

subroutine spst01 ( character  uplo,
integer  n,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( ldafac, * )  afac,
integer  ldafac,
real, dimension( ldperm, * )  perm,
integer  ldperm,
integer, dimension( * )  piv,
real, dimension( * )  rwork,
real  resid,
integer  rank 
)

SPST01

Purpose:
 SPST01 reconstructs a symmetric positive semidefinite matrix A
 from its L or U factors and the permutation matrix P and computes
 the residual
    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
 where EPS is the machine epsilon.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The number of rows and columns of the matrix A.  N >= 0.
[in]A
          A is REAL array, dimension (LDA,N)
          The original symmetric matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N)
[in]AFAC
          AFAC is REAL array, dimension (LDAFAC,N)
          The factor L or U from the L*L' or U'*U
          factorization of A.
[in]LDAFAC
          LDAFAC is INTEGER
          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
[out]PERM
          PERM is REAL array, dimension (LDPERM,N)
          Overwritten with the reconstructed matrix, and then with the
          difference P*L*L'*P' - A (or P*U'*U*P' - A)
[in]LDPERM
          LDPERM is INTEGER
          The leading dimension of the array PERM.
          LDAPERM >= max(1,N).
[in]PIV
          PIV is INTEGER array, dimension (N)
          PIV is such that the nonzero entries are
          P( PIV( K ), K ) = 1.
[out]RWORK
          RWORK is REAL array, dimension (N)
[out]RESID
          RESID is REAL
          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
[in]RANK
          RANK is INTEGER
          number of nonzero singular values of A.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 132 of file spst01.f.

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*
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
ssyr
subroutine ssyr(uplo, n, alpha, x, incx, a, lda)
SSYR
Definition ssyr.f:132
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slansy
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
strmv
subroutine strmv(uplo, trans, diag, n, a, lda, x, incx)
STRMV
Definition strmv.f:147
Here is the call graph for this function:
Here is the caller graph for this function: