      SUBROUTINE SGET1( N, AFAC, LDAFAC, IPIV, LDA, TMP, RESID )
      implicit none
*
*  -- LAPACK test routine --
*     Argonne National Lab, Courant Institute, and N.A.G. Ltd.
*     April 1, 1989
*
*  -- NPAC, Syracuse University --
*     Manish Parashar
*     April 1992
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDAFAC, N
      REAL               RESID, sum, slarnd
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), iseed(4)
      REAL               A, AFAC( LDAFAC, * ), TMP( * )
*     ..
*
*  Purpose
*  =======
*
*     SGET1 reconstructs the original matrix  A  from its
*     triangular factors  L  and  U  and computes the residual
*          norm( L*U - A ) / ( sqrt( N ) * norm( A ) * EPS ),
*     where EPS is the machine epsilon.
*
*  Arguments
*  ==========
*
*  N      - INTEGER
*           On entry, N specifies the order of the matrix A.
*           N should be at least 0.
*           Unchanged on exit.
*
*  AFAC   - REAL array of dimension( LDA, N )
*           On entry, AFAC contains the information about the LU
*           factorization of the original matrix.
*           Overwritten with the reconstructed matrix, and then
*           with the difference LU - A.
*
*  LDAFAC - INTEGER
*           On entry, LDAFAC specifies the first dimension of AFAC as
*           declared in the calling program.
*           LDAFAC should be at least N.
*           Unchanged on exit.
*
*  IPIV   - INTEGER array of dimension( N )
*           On entry, IPIV contains the pivot indices as determined
*           by SGETRF.
*           Unchanged on exit.
*
*  A      - REAL array of dimension( LDA, N )
*           On entry, A contains the original unfactored matrix.
*           Unchanged on exit.
*
*  LDA    - INTEGER
*           On entry, LDA specifies the first dimension of A as
*           declared in the calling program.
*           LDA should be at least N.
*           Unchanged on exit.
*
*  TMP    - REAL array of dimension( N )
*           TMP is an N-vector of work space required by the norm
*           routines.
*
*  RESID  - REAL
*           On return, RESID contains the ratio described above,
*           which measures the difference between the product LU
*           and the original matrix A.
*
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               ANORM
*     ..
*     .. External Functions ..
      REAL               R1MACH, SLANGE, SLAMCH
      EXTERNAL           R1MACH, SLANGE, SLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLAPGE
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, SQRT
*     ..
*     .. Executable Statements ..
*
*     Quick exit if N = 0.
*
      IF( N.LE.0 ) THEN
         RESID = ZERO
         RETURN
      END IF

c     print *,"LU Factorized Matrix"
c     call printmat(afac,n,n,"LU Factorized Matrix\0")

*
*     Call SLAPGE to reconstruct A from the triangular factors L and U .
*     The reconstructed matrix overwrites the array variable AFAC.
*
      CALL SLAPGE( N, AFAC, LDAFAC, IPIV )
*
*     Compute the difference  L*U - A  and store in AFAC.
*
      iseed(1) = 617
      iseed(2) = 1553
      iseed(3) = 2191
      iseed(4) = 3771
      anorm = 0.0
      DO 20 J = 1, N
	    SUM = 0.0
         DO 10 I = 1, N
            A = slarnd(2,iseed) 
c           write(*,*) 'a(i,j) = ',A
	    SUM = SUM + abs(A)
            AFAC( I, J ) = AFAC( I, J ) - A 
   10    CONTINUE
      ANORM = amax1(sum,anorm)
   20 CONTINUE
*
c     print *,"LU - A Matrix"
c     call printmat(afac,n,n,"LU-A matrix\0")

*     Compute norm( L*U - A ) / ( SQRT( N ) * NORM( A ) * EPS ) .
*
C     ANORM = SLANGE( '1', N, N, A, LDA, TMP )
      RESID = SLANGE( '1', N, N, AFAC, LDAFAC, TMP )
*
      RESID = ((RESID / SQRT(REAL(N))) / ANORM) / SLAMCH('Epsilon')
*
      RETURN
*
*     End of SGET1
*
      END

