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.

Date

December 2016

Definition at line 136 of file spst01.f.

*
*  -- LAPACK test routine (version 3.7.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     December 2016
*
*     .. Scalar Arguments ..
      REAL               RESID
      INTEGER            LDA, LDAFAC, LDPERM, N, RANK
      CHARACTER          UPLO
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), AFAC( LDAFAC, * ),
     $                   PERM( LDPERM, * ), RWORK( * )
      INTEGER            PIV( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      REAL               ANORM, EPS, T
      INTEGER            I, J, K
*     ..
*     .. External Functions ..
      REAL               SDOT, SLAMCH, SLANSY
      LOGICAL            LSAME
      EXTERNAL           sdot, slamch, slansy, lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           sscal, ssyr, strmv
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          real
*     ..
*     .. Executable Statements ..
*
*     Quick exit if N = 0.
*
      IF( n.LE.0 ) THEN
         resid = zero
         RETURN
      END IF
*
*     Exit with RESID = 1/EPS if ANORM = 0.
*
      eps = slamch( 'Epsilon' )
      anorm = slansy( '1', uplo, n, a, lda, rwork )
      IF( anorm.LE.zero ) THEN
         resid = one / eps
         RETURN
      END IF
*
*     Compute the product U'*U, overwriting U.
*
      IF( lsame( uplo, 'U' ) ) THEN
*
         IF( rank.LT.n ) THEN
            DO 110 j = rank + 1, n
               DO 100 i = rank + 1, j
                  afac( i, j ) = zero
  100          CONTINUE
  110       CONTINUE
         END IF
*
         DO 120 k = n, 1, -1
*
*           Compute the (K,K) element of the result.
*
            t = sdot( k, afac( 1, k ), 1, afac( 1, k ), 1 )
            afac( k, k ) = t
*
*           Compute the rest of column K.
*
            CALL strmv( 'Upper', 'Transpose', 'Non-unit', k-1, afac,
     $                  ldafac, afac( 1, k ), 1 )
*
  120    CONTINUE
*
*     Compute the product L*L', overwriting L.
*
      ELSE
*
         IF( rank.LT.n ) THEN
            DO 140 j = rank + 1, n
               DO 130 i = j, n
                  afac( i, j ) = zero
  130          CONTINUE
  140       CONTINUE
         END IF
*
         DO 150 k = n, 1, -1
*           Add a multiple of column K of the factor L to each of
*           columns K+1 through N.
*
            IF( k+1.LE.n )
     $         CALL ssyr( 'Lower', n-k, one, afac( k+1, k ), 1,
     $                    afac( k+1, k+1 ), ldafac )
*
*           Scale column K by the diagonal element.
*
            t = afac( k, k )
            CALL sscal( n-k+1, t, afac( k, k ), 1 )
  150    CONTINUE
*
      END IF
*
*        Form P*L*L'*P' or P*U'*U*P'
*
      IF( lsame( uplo, 'U' ) ) THEN
*
         DO 170 j = 1, n
            DO 160 i = 1, n
               IF( piv( i ).LE.piv( j ) ) THEN
                  IF( i.LE.j ) THEN
                     perm( piv( i ), piv( j ) ) = afac( i, j )
                  ELSE
                     perm( piv( i ), piv( j ) ) = afac( j, i )
                  END IF
               END IF
  160       CONTINUE
  170    CONTINUE
*
*
      ELSE
*
         DO 190 j = 1, n
            DO 180 i = 1, n
               IF( piv( i ).GE.piv( j ) ) THEN
                  IF( i.GE.j ) THEN
                     perm( piv( i ), piv( j ) ) = afac( i, j )
                  ELSE
                     perm( piv( i ), piv( j ) ) = afac( j, i )
                  END IF
               END IF
  180       CONTINUE
  190    CONTINUE
*
      END IF
*
*     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
*
      IF( lsame( uplo, 'U' ) ) THEN
         DO 210 j = 1, n
            DO 200 i = 1, j
               perm( i, j ) = perm( i, j ) - a( i, j )
  200       CONTINUE
  210    CONTINUE
      ELSE
         DO 230 j = 1, n
            DO 220 i = j, n
               perm( i, j ) = perm( i, j ) - a( i, j )
  220       CONTINUE
  230    CONTINUE
      END IF
*
*     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
*     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
*
      resid = slansy( '1', uplo, n, perm, ldafac, rwork )
*
      resid = ( ( resid / real( n ) ) / anorm ) / eps
*
      RETURN
*
*     End of SPST01
*

Here is the call graph for this function:

Here is the caller graph for this function: