LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 136 of file spst01.f.

136 *
137 * -- LAPACK test routine (version 3.4.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 * November 2011
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 *
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:53
strmv
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
Definition: strmv.f:149
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
ssyr
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:134
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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, or the element of largest absolute value of a real symmetric matrix.
Definition: slansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function: