LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
Definition: strmv.f:147
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:132
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: