Sequential LU Factorization

      SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
*
*  LU factorization of a M-by-N matrix A using partial pivoting with 
*  row interchanges.
*
      INTEGER            INFO, LDA, M, N, IPIV( * )
      DOUBLE PRECISION   A( LDA, * )
*
      INTEGER            I, IINFO, J, JB, NB
      PARAMETER          ( NB = 64 )
      EXTERNAL           DGEMM, DGETF2, DLASWP, DTRSM
      INTRINSIC          MIN
*
      DO 20 J = 1, MIN(M,N), NB
         JB = MIN( MIN(M,N)-J+1, NB )
*
*        Factor diagonal block and test for exact singularity.
*
         CALL DGETF2( M-J+1, JB, A(J,J), LDA, IPIV(J), IINFO )
*
*        Adjust INFO and the pivot indices.
*
         IF( INFO.EQ.0 .AND. IINFO.GT.0 ) INFO = IINFO + J - 1
         DO 10 I = J, MIN(M,J+JB-1)
            IPIV(I) = J - 1 + IPIV(I)
   10    CONTINUE
*
*        Apply interchanges to columns 1:J-1 and J+JB:N.
*
         CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
         IF( J+JB.LE.N ) THEN
            CALL DLASWP( N-J-JB+1, A(1,J+JB), LDA, J, J+JB-1, IPIV, 1 )
*
*           Compute block row of U and update trailing submatrix.
*
            CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
     $                  N-J-JB+1, 1.0D+0, A(J,J), LDA, A(J,J+JB), LDA )
            IF( J+JB.LE.M )
     $         CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1,
     $                     N-J-JB+1, JB, -1.0D+0, A(J+JB,J), LDA,
     $                     A(J,J+JB), LDA, 1.0D+0, A(J+JB,J+JB), LDA )
         END IF
   20 CONTINUE
      RETURN
*
      END