subroutine bl3dfac(n, k, E, D, F, pivot)

c	This version of bl3dfac dated March 15, 2000.

c************************************************************************

c	bl3dfac performs an LU factorization of the block tridiagonal
c	matrix given by:

c	D(1) E(1)  O   O    O ............  O
c	F(1) D(2) E(2) O    O ............  O
c	 O   F(2) D(3) E(3) O ............  O
c	 ....................................
c
c        O    O   ........F(n-1) D(n-1) E(n-1)
c        O    O   .............. F(n-1)   D(n)

c	where each block, E(j), D(j), F(j) is k by k.
c
c	Pivoting is done only within the diagonal blocks.

c	The L matrix is NOT unit lower triangular, but the U matrix
c	is unit UPPER triangular, with the identity matrix on the
c	diagonal. 

c************************************************************************

c	Uses: Lapack routines dgetrf and dgetrs for Gauss LU factorization
c	and solution, and the BLAS routine dgemm for matrix-matrix 
c       multiplication.

c************************************************************************

c	Input variables:

	integer n, k
	double precision E(k,k,n-1), D(k,k,n), F(k,k,n-1)
	integer pivot(k,n)

c	Input:

c	k           Integer, the size of the blocks.
c	n           Integer, the number of blocks.
c	E(k,k,n-1)  Double precision, the sub-diagonal blocks.
c	D(k,k,n)    Double precision, the diagonal blocks.
c	F(k,k,n-1)  Double precision, the super-diagonal blocks.

c	Output:

c	E(k,k,n-1)  The lower triangular blocks of L in the LU factorization.
c	F(k,k,n-1)  The upper triangular blocks of U in the LU factorization.
c	D(k,k,n)    The diagonal blocks of L in the LU factorization.
c                   The diagonal blocks of U are unit upper triangular.
c	pivot(k,n)  Integer pivot vector used in factoring the diagonal blocks.
c		    pivot(1:k,p) records the pivoting done in diagonal block p.

c	Local variables:

	integer j, info
	character *1 trans
	double precision one
	parameter ( one = 1.0d0 )

c************************************************************************

c	Executable statements.

c************************************************************************

c	First, in main blocks 1 to n-1:

	trans = 'N'    
        
	do 100 j = 1,n-1

c          First, factor D(j).
           
	   call dgetrf(k, k, D(1,1,j), k, pivot(1,j), info)

	   if ( info .ne. 0 ) then
	      print *,'At block ',j,','
	      print *,'problem, info = ', info   ! make this better later.
              return
	   endif

c          Now, compute E(j) from D(j) * E(j) = E(j).
	   
	   call dgetrs(trans, k, k, D(1,1,j), k, pivot(1,j),
     *                 E(1,1,j), k, info)
	   if ( info .lt. 0 ) then
	      print *,'Illegal parameter number ',-info
	      return
	   endif

c	   Finally, compute D(j+1) = D(j+1) - F(j) * E(j).

	   call dgemm(trans, trans, k, k, k, -one, F(1,1,j), k,
     *                E(1,1,j), k, one, D(1,1,j+1), k)

  100 	continue

c	Finally, obtain the LU factorization of D(n).

	call dgetrf(k, k, D(1,1,n), k, pivot(1,n), info)
	if ( info .ne. 0 ) then
	   print *,'At block ',j,','
	   print *,'problem, info = ', info   ! make this better later.
           return
	endif

c************************************************************************

	return
	end

c************************ END OF bl3dfac ********************************