Next: Parallel LU Factorization Up: Translating LAPACK-based programs to Previous: Translating LAPACK-based programs to

## 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