      subroutine kji_block( m, n, a, lda, ipiv, info, nb )
*
*  -- Blocked KJI-SAXPY method to illustrate the use of Level 3 BLAS
*     Northeast Parallel Architecture Center, Syracuse University.
*     Dec. 26, 1990
*     A. Gaber Mohamed, NPAC, Syracuse University, Syracuse, New York.
*     .. scalar arguments ..
      integer            info, lda, m, n
*     ..
*     .. array arguments ..
      integer            ipiv( * )
      double precision   a( lda, * )
*     ..
*
*  purpose
*  =======
*
*  kji_block computes an lu factorization of a general m-by-n matrix a
*  using partial pivoting with row interchanges.
*
*  At the kth step of the elemination, a block column of both matrices U and L
*  is comuted.
*
*  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).
*
*
*  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.
*
*  a       (input/output) double precision array, dimension (lda,n)
*          on entry, the m x n matrix to be factored.
*          on exit, the factors l and u; the unit diagonal elements of l
*          are not stored.
*
*  lda     (input) integer
*          the leading dimension of the array a.  lda >= max(1,m).
*
*  ipiv    (output) integer array, dimension (min(m,n))
*          the pivot indices.  row i of the matrix was interchanged with
*          row ipiv(i).
*
*  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.
*
*  =====================================================================
*
*     .. parameters ..
      double precision   one
      parameter          ( one = 1.0d+0 )
*     ..
*     .. local scalars ..
      integer            i, k, kb, nb
*     ..
*     .. external subroutines ..
      external           dgemm, dgetf2, dlaswp, dtrsm, xerbla
*     ..
*     .. intrinsic functions ..
      intrinsic          max, min
*     ..
*     .. executable statements ..
*
*     test the input parameters.
*
      info = 0
      if( m.lt.0 ) then
         info = -1
      else if( n.lt.0 ) then
         info = -2
      else if( lda.lt.max( 1, m ) ) then
         info = -4
      end if
      if( info.ne.0 ) then
         call xerbla( 'kji_block', -info )
         return
      end if
*
*     quick return if possible.
*
      if( m.eq.0 .or. n.eq.0 )
     $   return
*
*
      if( nb.eq.1 ) then
         call jik_noblock( m, n, a, lda, ipiv, info )
      else
*
         do 20 k = 1, min( m, n ), nb
            kb = min( min( m, n )-k+1, nb )
*
*           factorize diagonal and subdiagonal blocks 
*
         call jik_noblock( m-k+1, kb, a( k, k ), lda, ipiv( k ), info )
*
*           update pivot indices and apply the interchanges to 
*           previous and following blocks of the current block.
*
            do 15 i = k, min( m, k+kb-1 )
               ipiv( i ) = k - 1 + ipiv( i )
   15       continue
         call dlaswp( k-1, a, lda, k, k+kb-1, ipiv, 1 )
         call dlaswp( n-k-kb+1, a( 1, k+kb ), lda, k, k+kb-1, ipiv, 1 )
*
            if( info.ne.0 ) go to 30

*           Compute block row of U

               call dtrsm( 'left', 'lower', 'no transpose', 'unit', kb,
     $                     n-k-kb+1, one, a( k, k ), lda, a(k, k+kb), lda )

*
*       Update right-hand bottom, (n-k-kb+1) X (n-k-kb+1) block

               call dgemm( 'no transpose', 'no transpose', n-k-kb+1,
     $                     n-k-kb+1, kb, -one, a(k+kb,k), lda, 
     $                     a(k,k+kb), lda, one, a(k+kb,k+kb), lda)
   20    continue
      end if
      return
*
   30 continue
      info = info + k - 1
      return
*
*     end of kji_block
*
      end
