subroutine bl3dsolmod(n, k, nk, E, D, F, pivot, nrhs, rhs)

c	This version of bl3dsolmod dated April 13, 2000.

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

c	bl3dsol solves the system of linear equations whose coefficient
c	matrix is the block tridiagonal	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, and which has already
c	been factored using bl3dfac.

c	The right hand sides are stored by with each right hand side
c	vector consecutive (see below).

c	The L matrix is NOT unit lower triangular, but the U matrix
c	is unit UPPER triangular.

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

c	Uses: 

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

c	Input variables:

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

c	Input:

c	k             Integer, the size of the blocks.
c	n             Integer, the number of blocks.
c       nk            Integer, the number of equations.
c	E(k,k,n-1)    Double precision, the sub-diagonal blocks.
c	D(k,k,n)      Double precision, the diagonal blocks, after factoring 
c                     by bl3dfac.
c	F(k,k,n-1)    Double precision, the super-diagonal blocks, after 
c                     modification in bl3dfac.
c	pivot(k,n)    Integer pivot vector returned by bl3dfac, used in 
c                     factoring the diagonal blocks.
c		      pivot(1:k,p) records the pivoting done in diagonal 
c                     block p.
c	nrhs          integer, the number of right hand sides.
c	rhs(nk,nrhs) Double precision, right hand sides. These are stored
c	              by blocks - that is, the nrhs right hand sides 
c	              associated with each k by k block are stored 
c	              consecutively. For example, for n = 3, k = 2 and
c                     nrhs = 3, the elements of the right hand sides
c                     are stored in the following order:

c	              [1] [7 ] [13]
c                     [2] [8 ] [14]  First block.

c                     [3] [9 ] [15]
c                     [4] [10] [16]  Second block.

c	              [5] [11] [17]  Third block.
c                     [6] [12] [18]

c	Output:

c	rhs(nk,nrhs) double precision, the solutions.
c	
c	Local variables:

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

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

c	Executable statements.

c************************************************************************
        
	do 1000 nr = 1, nrhs

c	   Modification of the right hand sides, that is, solve the
c	   block bi-diagonal lower triangular system L.x = b, overwriting
c	   b with the solution.

c	   First, solve D(1).x(1..k) = b(1..k).

	   trans = 'N'
	   call dgetrs(trans, k, 1, D(1,1,1), k, pivot(1,1),
     *                 rhs(1,nr), k, info)

	   if ( info .lt. 0 ) then
	      print *,'Illegal parameter number ',-info
	      return
	   endif

c	   Then, for j = 2 to n solve 
c	   D(j)x((j-1)*k+1..j*k) = b((j-1)*k+1..j*k) 
c                                - F(j-1)x((j-2)*k+1..(j-1)*k)

	   do 100 j = 2, n
           
	      call dgemm(trans, trans, k, 1, k, -one, F(1,1,j-1), k,
     *                  rhs(k*(j-2)+1,nr), k, one, rhs(k*(j-1)+1,nr), k)

	      call dgetrs(trans, k, 1, D(1,1,j), k, pivot(1,j),
     *                    rhs(k*(j-1)+1,nr), k, info)
	      if ( info .lt. 0 ) then
	         print *,'Illegal parameter number ',-info
	         return
	      endif
	   

  100 	   continue

	   do 200 j = n-1,1,-1

c	      Form rhs((j-1)*k+1:j*k,nr) = 
c	           rhs((j-1)*k+1:j*k,nr) - E(j)*rhs(j*k+1:(j+1)*k,nr)

	      call dgemm(trans, trans, k, 1, k, -one, E(1,1,j), k,
     *                   rhs(j*k+1,nr), k, one, rhs((j-1)*k+1,nr), k)

  200	   continue

 1000	continue

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

	return
	end

c************************ END OF bl3dsol ********************************