      SUBROUTINE PSGEIN3( N, A, DESCA, IPIV, WORK, IWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.4) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*     .. Scalar Arguments ..
      INTEGER            INFO, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), IPIV( * ), IWORK( * )
      REAL   A( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PSGEIN3 computes ...
*
*  This is the right-looking Parallel Level 2 BLAS version of the
*  algorithm.
*
*  Arguments
*  =========
*
*  M       (global input) INTEGER
*          The number of rows to be operated on, i.e. the number of rows
*          of the distributed submatrix sub( A ). M >= 0.
*
*  N       (global input) INTEGER
*          The number of columns to be operated on, i.e. the number of
*          columns of the distributed submatrix sub( A ). N >= 0.
*
*  A       (local input/local output) REAL pointer into the
*          local memory to an array of dimension (LLD_A, LOCq(JA+N-1)).
*          On entry, this array contains the local pieces of the M-by-N
*          distributed matrix sub( A ) to be factored. On exit, this
*          array contains the local pieces of the factors L and U from
*          the factorization sub( A ) = P*L*U; the unit diagonal ele-
*          ments of L are not stored.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  IPIV    (local output) INTEGER array, dimension ( LOCp(M_A)+MB_A )
*          This array contains the pivoting information.
*          IPIV(i) -> The global row local row i was swapped with.
*          This array is tied to the distributed matrix A.
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*          > 0:  If INFO = K, U(IA+K-1,JA+K-1) is exactly zero.
*                The factorization has been completed, but the factor U
*                is exactly singular, and division by zero will occur if
*                it is used to solve a system of equations.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL   ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          COLBTOP, COLCTOP, ROWBTOP
      INTEGER            I, ICOFF, ICTXT, IINFO, IN, IROFF, J, JB, JN,
     $                   MN, MYCOL, MYROW, NPCOL, NPROW,
     $                   ITEMP, IAROW, IACOL, IPB, IPR, MP, K,
     $                   II, JJ, LI, LJ, IDIR
*     ..
*     .. Local Arrays ..
      INTEGER            DESCB( DLEN_ ), DESCR( DLEN_ ),
     $                   IDUM1( 1 ), IDUM2( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, IGAMN2D, PCHK1MAT,
     $                   PTOPGET, PTOPSET, PSGEMM, PSGETF2,
     $                   PSLASWP, PSTRSM, PXERBLA, INFOG1L, INFOG2L,
     $                   DESCINIT
*     ..
*     .. External Functions ..
      INTEGER            ICEIL, INDXG2P, NUMROC
      EXTERNAL           ICEIL, INDXG2P, NUMROC
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(600+CTXT_)
      ELSE
         CALL CHK1MAT( N, 1, N, 2, 1, 1, DESCA, 6, INFO )
         IF( INFO.EQ.0 ) THEN
            IROFF = MOD( 1-1, DESCA( MB_ ) )
            ICOFF = MOD( 1-1, DESCA( NB_ ) )
            IF( IROFF.NE.0 ) THEN
               INFO = -4
            ELSE IF( ICOFF.NE.0 ) THEN
               INFO = -5
            ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(600+NB_)
            END IF
         END IF
         CALL PCHK1MAT( N, 1, N, 2, 1, 1, DESCA, 6, 0, IDUM1,
     $                  IDUM2, INFO )
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PSGEIN3', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( DESCA( M_ ).EQ.1 ) THEN
         IPIV( 1 ) = 1
         RETURN
      ELSE IF( N.EQ.0 ) THEN
         RETURN
      END IF
*
*     Split-ring topology for the communication along process rows
*
      CALL PTOPGET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP )
      CALL PTOPGET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP )
      CALL PTOPGET( ICTXT, 'Combine', 'Columnwise', COLCTOP )
      CALL PTOPSET( ICTXT, 'Broadcast', 'Rowwise', 'S-ring' )
      CALL PTOPSET( ICTXT, 'Broadcast', 'Columnwise', ' ' )
      CALL PTOPSET( ICTXT, 'Combine', 'Columnwise', ' ' )
*
      IPR = 1
      IPB = IPR + DESCA( NB_ )
*
*     Loop over the N columns.
*
      DO J = 1, N, DESCA( NB_ )
        JB = MIN( N-J+1, DESCA( NB_ ) )
        I = J
*
*       Process block A(1:N,J:J+JB).
*       ============================
*
        CALL PSGEI2( N-I+1, JB, A, J, J, DESCA, IPIV, WORK( IPR ), 
     $               IINFO )
*
        IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 
     $    INFO = IINFO + J - 1
*
*       Substracting 1's to the diagonal elements.
*       ==========================================
*
        CALL INFOG1L( J, DESCA( NB_ ), NPCOL, MYCOL, DESCA( CSRC_ ),
     $                LJ, IACOL )
        IF( MYCOL.EQ.IACOL ) THEN
 
          DO K = J, J+JB-1
*
*           Get local indices for I and J, as well as the row number of
*           processors that contain the I-th row.
*
            CALL INFOG1L( K, DESCA( MB_ ), NPROW, MYROW, DESCA( RSRC_ ),
     $                    LI, IAROW )
*
*           Find pivot and test for singularity.
*
            IF( MYROW.EQ.IAROW ) THEN
              IDIR = LI+( LJ+K-J-1 )*DESCA( LLD_ )
              A( IDIR ) = A( IDIR )-ONE
            END IF

          END DO
        END IF
*
*       Apply the interchanges.
*       =======================
*
*       Apply the interchanges to columns J+JB:N.
*
        IF( J+JB.LE.N )
     $    CALL PSLASWP( 'Forward', 'Rowwise', N-J-JB+1, A, 1, J+JB, 
     $                  DESCA, I, I+JB-1, IPIV )
*
*       Apply the interchanges to columns 1:J-1.
*
        IF( J.GT.1 )
     $    CALL PSLASWP( 'Forward', 'Rowwise', J-1, A, 1, 1, DESCA,
     $                  I, I+JB-1, IPIV )
*
*
*       Apply the transformations to columns J+JB:N.
*       ============================================
*
*       Define an array descriptor for workspace block.
*
        CALL INFOG2L( I, J+JB, DESCA, NPROW, NPCOL, MYROW, MYCOL, 
     $                II, JJ, IAROW, IACOL )
*
        CALL DESCINIT( DESCB, DESCA( NB_ ), N-J-JB+1,
     $                 DESCA( NB_ ), DESCA( NB_ ),
     $                 IAROW, IACOL,
     $                 ICTXT, MAX( 1, DESCA( NB_ ) ), IINFO )
*
*       Copy block A(J:J+JB,J+JB:N).
*
        CALL PSLACPY( 'All', JB, N-J-JB+1, A, J, J+JB, DESCA,
     $                  WORK( IPB ), 1, 1, DESCB )
C        CALL PSLAPRN2( ..., ..., WORK( IPB ), 1, 1, DESCB, 0, 0,
C     $               'www', 6, WORK( 10000 ) )
*
*       Multiply blocks...
*
        CALL PSGEMM( 'No tranpose', 'No tranpose', N, N-J-JB+1, JB, 
     $               ONE, A, 1, J, DESCA, 
     $               WORK( IPB ), 1, 1, DESCB,
     $               ONE, A, 1, J+JB, DESCA )
*
*
*       Apply the transformations to columns 1:J-1.
*       ===========================================
*
*       Define an array descriptor for workspace block.
*
        CALL INFOG2L( I, 1, DESCA, NPROW, NPCOL, MYROW, MYCOL, 
     $                II, JJ, IAROW, IACOL )
*
        CALL DESCINIT( DESCB, DESCA( NB_ ), J-1,
     $                 DESCA( NB_ ), DESCA( NB_ ),
     $                 IAROW, IACOL,
     $                 ICTXT, MAX( 1, DESCA( NB_ ) ), IINFO )
*
*       Copy block A(J:J+JB,1:J-1).
*
        CALL PSLACPY( 'All', JB, J-1, A, J, 1, DESCA,
     $                WORK( IPB ), 1, 1, DESCB )
*
*       Multiply blocks...
*
        CALL PSGEMM( 'No tranpose', 'No tranpose', N, J-1, JB, 
     $               ONE, A, 1, J, DESCA, 
     $               WORK( IPB ), 1, 1, DESCB,
     $               ONE, A, 1, 1, DESCA )
*
*
*       Adding 1's to the diagonal elements.
*       ====================================
*
        CALL INFOG1L( J, DESCA( NB_ ), NPCOL, MYCOL, DESCA( CSRC_ ),
     $                LJ, IACOL )
        IF( MYCOL.EQ.IACOL ) THEN
 
          DO K = J, J+JB-1
*
*           Get local indices for I and J, as well as the row number of
*           processors that contain the I-th row.
*
            CALL INFOG1L( K, DESCA( MB_ ), NPROW, MYROW, DESCA( RSRC_ ),
     $                    LI, IAROW )
*
*           Find pivot and test for singularity.
*
            IF( MYROW.EQ.IAROW ) THEN
              IDIR = LI+( LJ+K-J-1 )*DESCA( LLD_ )
              A( IDIR ) = A( IDIR )+ONE
            END IF
 
          END DO
        END IF
*
      END DO
*
*
*     Use the row pivots and apply them to the columns of the global
*     matrix.
*
*
      MP = NUMROC( DESCA( M_ ), DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 
     $             NPROW )
* 
      CALL DESCSET( DESCR, DESCA( M_ ) + DESCA( MB_ )*NPROW, 1,
     $              DESCA( MB_ ), 1, DESCA( RSRC_ ), MYCOL, ICTXT,
     $              MP+DESCA( MB_ ) )
     CALL PSLAPIV( 'Backward', 'Columns', 'Column', N, N, A, 1,
    $              1, DESCA, IPIV, 1, 1, DESCR, IWORK )
*
*
*
      IF( INFO.EQ.0 )
     $   INFO = MN + 1
      CALL IGAMN2D( ICTXT, 'Rowwise', ' ', 1, 1, INFO, 1, IDUM1, IDUM2,
     $              -1, -1, MYCOL )
      IF( INFO.EQ.MN+1 )
     $   INFO = 0
*
      CALL PTOPSET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP )
      CALL PTOPSET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP )
      CALL PTOPSET( ICTXT, 'Combine', 'Columnwise', COLCTOP )
*
      RETURN
*
*     End of PSGEIN3
*
      END
