      SUBROUTINE PZLUBR( M, N, NB, A, LDA, WORK, LDWORK, IPIV, INFO )
*
*     Parallel blocked LU factorization
*
*
*
*     .. Scalar Arguments ..
      INTEGER            LDA, M, N, NB, LDWORK, INFO
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      COMPLEX*16         A( LDA, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PZLUBR computes an LU factorization of a general m-by-n matrix A
*  using partial pivoting with row interchanges.  It assumes this
*  node's portion of the matrix is in array A.  A block wrapped 
*  storage scheme with blocksize nb is assumed.
*
*  The factorization has the form
*     A = P * L * U
*  where P is a permutation matrix, L is lower triangular with unit
*  diagonal elements (lower trapezoidal if m > n), and U is upper
*  triangular (upper trapezoidal if m < n).
*
*  This is based on the right-looking Level 3 BLAS LAPACK version
*  of the algorithm.
*
*  NOTE: not currently fully tested for N .ne. M
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrix A.  N >= 0.
*
*  NB      (input) INTEGER
*          The block size for the blocked algorithm.  NB > 0.
*
*  A       (input/output) COMPLEX*16 array, dimension (LDA,LOCALN)
*          On entry, this node's portion of m by n matrix to be factored.
*          On exit, this node's factors L and U; the unit diagonal 
*          elements of L are not stored.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A that holds the local
*          portion of the matrix to be factored.  LDA >= max(1,LOCALM).
*
*  WORK    COMPLEX*16 array, dimension LDWORK
*          Work array.
*
*  LDWORK  (input) INTEGER
*          Size of work array WORK.  LDWORK >= (LOCALM + LOCALN + NB )*NB
*
*  IPIV    (output) INTEGER array, dimension (min(M,N))
*          The pivot indices.  Row i of the matrix was interchanged with
*          row IPIV(i) if IPIV(i) is a row that resides on this node.
*          otherwise -IPIV(i) equals the row index of the node that
*          holds the pivot row.
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     VERSION 2.0 DATED Aug. 1992
*     R. VAN DE GEIJN
*
*
*        This code resulted from a joint collaboration between  
*              Robert van de Geijn 
*     and the group for 'Research in Parallel Computational
*     Electromagnetics' at the Jet Propulsion Laboratory.  Futher 
*     information can be found in the paper,
*
*         T. Cwik, R. van de Geijn, and J. Patterson, 
*         "Application of Massively Parallel Computation to 
*          Integral Equation Models of Electromagnetic Scattering 
*          (Invited Paper)," J. Opt. Soc. Am. A, Vol. 11, No. 4, 
*          April 1994.
*
*     Please email your comments to the authors at
*
*        rvdg@cs.utexas.edu (Robert van de Geijn)
*        cwik@yosemite.jpl.nasa.gov (Tom Cwik)
*        jep@yosemite.jpl.nasa.gov  (Jean Patterson)
*
*     Use of this code must properly reference that paper.
*
*     The code is free, and comes with no implicit or explicit 
*     guarantees.  
*
*     All rights reserved
*
*     .. Parameters ..
      COMPLEX*16         ONE
      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
*     ..
*     .. Local Scalars ..
*
*     ..
*     .. Misc. ..
      INTEGER            I, IINFO, J, JB
*
*     nprow          row dimension of node grid
*     npcol          column dimension of node grid
*     myrow          my row index
*     mycol          my column index
*
      INTEGER           NPROW, NPCOL, MYROW, MYCOL
      INTEGER           LOCALM, LOCALN
      INTEGER           ICURROW, ICURCOL, IWORK2
      INTEGER           II, JJ, IDUMMY
*     ..
*     .. External Subroutines ..
      EXTERNAL           ZGEMM, ZLASWP, PZLUR, ZTRSM, XERBLA
*
      INTEGER            ITYPE_ROW, ITYPE_COL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0

*
*     get machine parameters
*
      CALL PLAMCH2( NPROW, NPCOL, MYROW, MYCOL )
*
*     determine local number of rows and columns
*
      CALL IMYPART( 1, M, NB, IDUMMY, LOCALM, MYROW, NPROW )
      CALL IMYPART( 1, N, NB, IDUMMY, LOCALN, MYCOL, NPCOL )
      
*
*     check input parameters
*
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, LOCALM ) ) THEN
         INFO = -5
      ELSE IF( NB.LE.0 ) THEN
         INFO = -3
      ELSE IF (LDWORK .LT. (LOCALM+NB+LOCALN)*NB ) THEN
         INFO =-7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'PZLUBR ', -INFO )
         RETURN
      END IF

      IWORK2 = NB*(LOCALM+NB)+1
*     
      ICURROW = 0
      ICURCOL = 0

      II = 1
      JJ = 1
      DO 20 J = 1, MIN( M, N ), NB
         JB = MIN( MIN( M, N )-J+1, NB )
*
*           Factor diagonal and subdiagonal blocks and test for exact
*           singularity.
*
         IF (MYCOL .EQ. ICURCOL) THEN
            CALL PZLUR( M-J+1, JB, A( II, JJ ), LDA, IPIV( J ), 
     $           ICURROW, WORK, NB, NB, 
     $           IINFO )
         ENDIF
*
*           Adjust INFO and the pivot indices.
*
         IF( INFO.EQ.0 .AND. IINFO.GT.0 )
     $        INFO = IINFO + J - 1
         DO 10 I = J, J + JB - 1
            IF (IPIV(I) .GT. 0) 
     $           IPIV( I ) = II - 1 + IPIV( I )
 10      CONTINUE
*
*     broadcast pivot info
*
         CALL IGEBC2D( 'RING', 'ROW', JB, 1, IPIV(J), JB,            
     $        MYROW, ICURCOL, ITYPE_ROW())
*           Apply interchanges to columns outside this block.
*
         IF (MYCOL .NE. ICURCOL) THEN
            CALL PZLASWP( LOCALN, A, LDA, J, J+JB-1, II, 
     $           ICURROW, IPIV, 1 )
         ELSE
            IF (JJ .GT. 1) THEN
               CALL PZLASWP( JJ-1, A, LDA, J, J+JB-1, II, 
     $              ICURROW, IPIV, 1 )
            ENDIF
            IF (JJ+JB .LE. LOCALN) THEN
               CALL PZLASWP( LOCALN-JJ-JB+1, A(1,JJ+JB), LDA, 
     $              J, J+JB-1, II, ICURROW, IPIV, 1)
            ENDIF
         ENDIF
*
         IF (J+JB .LE. N) THEN
*
*     broadcast factored panel
*
            IF (MYROW .EQ. ICURROW) II = II + JB
            IF (MYCOL .EQ. ICURCOL) 
*
*     copy this node's part of factored panel into work array
*
     $           CALL ZLACPY( 'G', LOCALM-II+1, JB, A( II, JJ ),
     $           LDA, WORK( NB*NB+1 ), LOCALM-II+1 )
*     
*     NOTICE: for some rows, sending too much
*
            CALL ZGEBCS ('CTREE', 'ROW', LOCALM-II+1+NB, JB, 
     $           WORK, LOCALM-II+1+NB, MYROW, ICURCOL, 
     $           ITYPE_ROW() )
            IF (MYROW .EQ. ICURROW) II = II-JB
         ENDIF
*            
         IF (MYCOL .EQ. ICURCOL) JJ = JJ+JB
*
         IF (JJ .LE. LOCALN) THEN
*     
*     update horizontal panel
*
            IF ( MYROW .EQ. ICURROW) THEN
               CALL ZTRSM( 'LEFT', 'LOWER', 'NO TRANSPOSE', 
     $              'UNIT', JB, LOCALN-JJ+1, ONE, WORK, NB,
     $              A( II, JJ ), LDA )
               CALL ZLACPY( 'G', JB, LOCALN-JJ+1, A( II, JJ ),
     $              LDA, WORK(IWORK2), JB )
               CALL ZGEBCS ( 'CTREE', 'COLUMN', JB, LOCALN-JJ+1, 
     $              WORK(IWORK2), JB, ICURROW, MYCOL,
     $              ITYPE_COL() )
               II = II+JB
            ELSE
               CALL ZGEBCS ( 'CTREE', 'COLUMN', JB, LOCALN-JJ+1, 
     $              WORK(IWORK2), JB, ICURROW, MYCOL,
     $              ITYPE_COL() )
            ENDIF
*
            IF( II .LE. LOCALM ) THEN
*
*                 Update trailing submatrix.
*
               CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', 
     $              LOCALM-II+1, LOCALN-JJ+1, JB, -ONE, 
     $              WORK(NB*NB+1), 
     $              LOCALM-II+1, WORK(IWORK2), JB, 
     $              ONE, A(II, JJ),
     $              LDA )
            END IF
         ELSEIF (MYROW .EQ. ICURROW) THEN
            II = II+JB
         END IF
         ICURROW = MOD( ICURROW+1, NPROW )
         ICURCOL = MOD( ICURCOL+1, NPCOL )
 20      CONTINUE
      RETURN
*
*     End of PZLUBR
*
      END
