SUBROUTINE DGEFAM (A,LDA,N,IPVT,INFO)
      INTEGER LDA,N,IPVT(N),INFO
      DOUBLE PRECISION A(LDA,N)
C
C     DGEFAM FACTORS A DOUBLE PRECISION REAL MATRIX BY
C     GAUSSIAN ELIMINATION.
C
C     DGEFAM IS USUALLY CALLED BY DGECOM, BUT IT CAN BE CALLED
C     DIRECTLY WITH A MODEST SAVING IN TIME IF RCOND IS NOT
C     NEEDED.
C
C     ON ENTRY:
C
C        A       DOUBLE PRECISION(LDA,N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY A AS DECLARED
C                IN THE MAIN CALLING PROGRAM.
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX A.
C
C     ON RETURN:
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L * U WHERE
C                L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C                   IPVT(K) = THE INDEX OF THE K-TH PIVOT ROW
C                THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY
C                   DET(A) = S*A(1,1)*A(2,2)* ... *A(N,N)
C                WHERE S = (-1)**(NUMBER OF TIMES IPVT(K) .NE. K)
C                BUT THIS RESULT SHOULD BE USED WITH CAUTION AS IT
C                MAY EASILY UNDERFLOW OR OVERFLOW.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0D0.  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT DGESLM WILL DIVIDE BY ZERO IF
C                     CALLED.
C                     USE RCOND IN DGECOM FOR A RELIABLE INDICATION OF
C                     SINGULARITY.
C
C     THIS VERSION IS ADAPTED FROM THE LINPACK SUBROUTINE DGEFA BY
C     REPLACEMENT OF CALLS TO THE BLAS WITH IN-LINE CODE.  MODIFICATION
C     DONE BY ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA,
C     APRIL 1980.
C
C     FORTRAN FUNCTIONS CALLED:
C
      DOUBLE PRECISION DABS
C
C     INTERNAL VARIABLES:
C
      INTEGER I,J,K,KP1,L,NM1
      DOUBLE PRECISION T
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      IPVT(N) = 1
      INFO = 0
      NM1 = N-1
      IF (NM1 .LT. 1) GO TO 100
      DO 90 K=1,NM1
         KP1 = K+1
C
C        FIND L = PIVOT INDEX
C
         L = K
         DO 10 I=KP1,N
            IF (DABS(A(I,K)) .GT. DABS(A(L,K))) L = I
   10    CONTINUE
         IPVT(K) = L
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (A(L,K) .EQ. 0.0D0) GO TO 70
C
C        INTERCHANGE IF NECESSARY
C
         IF (L .EQ. K) GO TO 20
         T = A(L,K)
         A(L,K) = A(K,K)
         A(K,K) = T
   20    CONTINUE
C
C        COMPUTE MULTIPLIERS
C
         T = 1.0D0/A(K,K)
         DO 30 I=KP1,N
            A(I,K) = T*A(I,K)
   30    CONTINUE
C
C        ROW ELIMINATION WITH COLUMN INDEXING
C
         DO 60 J=KP1,N
            T = A(L,J)
            IF (L .EQ. K) GO TO 40
            A(L,J) = A(K,J)
            A(K,J) = T
   40       CONTINUE
            DO 50 I=KP1,N
               A(I,J) = A(I,J)-T*A(I,K)
   50       CONTINUE
   60    CONTINUE
         GO TO 80
   70    CONTINUE
         INFO = K
   80    CONTINUE
   90 CONTINUE
  100 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D0) INFO = N
      RETURN
      END