Parallel LU Factorization

 

      SUBROUTINE PDGETRF( M, N, A, IA, JA, DESCA, IPIV, INFO )
*
      INTEGER            IA, INFO, JA, M, N, DESCA( * ), IPIV( * )
      DOUBLE PRECISION   A( * )
*
*  LU factorization of a M-by-N distributed matrix A(IA:IA+M-1,JA:JA+N-1)
*  using partial pivoting with row interchanges.
*
      INTEGER            I, IINFO, J, JB
      EXTERNAL           IGAMN2D, PDGEMM, PDGETF2, PDLASWP, PDTRSM
      INTRINSIC          MIN
*
      DO 10 J = JA, JA+MIN(M,N)-1, DESCA( NB_ )
         JB = MIN( MIN(M,N)-J+JA, DESCA( NB_ ) )
         I = IA + J - JA
*
*        Factor diagonal block and test for exact singularity.
*
         CALL PDGETF2( M-J+JA, JB, A, I, J, DESCA, IPIV, IINFO )
         IF( INFO.EQ.0 .AND. IINFO.GT.0 ) INFO = IINFO + J - JA
*
*        Apply interchanges to columns JA:J-JA and J+JB:JA+N-1.
*
         CALL PDLASWP( 'Forward', 'Rows', J-JA, A, IA, JA, DESCA,
     $                 I, I+JB-1, IPIV )
         IF( J-JA+JB+1.LE.N ) THEN
            CALL PDLASWP( 'Forward', 'Rows', N-J-JB+JA, A, IA, J+JB,
     $                    DESCA, I, I+JB-1, IPIV )
*
*           Compute block row of U and update trailing submatrix.
*
            CALL PDTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
     $                   N-J-JB+JA, 1.0D+0, A, I, J, DESCA, A, I, J+JB,
     $                   DESCA )
            IF( J-JA+JB+1.LE.M ) THEN
     $         CALL PDGEMM( 'No transpose', 'No transpose', M-J-JB+JA,
     $                      N-J-JB+JA, JB, -1.0D+0, A, I+JB, J, DESCA, A,
     $                      I, J+JB, DESCA, 1.0D+0, A, I+JB, J+JB, DESCA )
         END IF
   10 CONTINUE
      IF( INFO.EQ.0 ) INFO = MIN(M,N) + 1
      CALL IGAMN2D( ICTXT, 'Row', ' ', 1, 1, INFO, 1, I, J, -1, -1, MYCOL )
      IF( INFO.EQ.MIN(M,N)+1 ) INFO = 0
*
      RETURN
*
      END

The required steps to call a ScaLAPACK routine from a parallel program are demonstrated in Example Program #1 in Chapter 2 and explained in section 2.3.