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

Susan Blackford
Tue May 13 09:21:01 EDT 1997