subroutine bl3dsol(trans, n, k, E, D, F, pivot, nrhs, rhs) c This version of bl3dsol dated April 20, 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 blocks (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: character *1 trans integer n, k double precision E(k,k,n-1), D(k,k,n), F(k,k,n-1) integer pivot(k,n), nrhs double precision rhs(k,nrhs,n) c Input: c trans character*1, 'n' or 'N' is A.x = b is to be solved, c 't' or 'T' is A^T.x = b is to be solved. c n Integer, the number of blocks. c k Integer, the size of the blocks. 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(k,nrhs,n) 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 ] [3 ] [5 ] c [2 ] [4 ] [6 ] First block. c [7 ] [9 ] [11] c [8 ] [10] [12] Second block. c [13] [15] [17] Third block. c [14] [16] [18] c Output: c rhs(k,nrhs,n) double precision, the solutions. c c Local variables: integer j, info double precision one parameter ( one = 1.0d0 ) c************************************************************************ c Executable statements. c************************************************************************ 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. if ( trans .eq. 'n' .or. trans .eq. 'N' ) then c First, solve D(1).x(1..k) = b(1..k). call dgetrs(trans, k, nrhs, D(1,1,1), k, pivot(1,1), * rhs(1,1,1), 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, nrhs, k, -one, F(1,1,j-1), k, * rhs(1,1,j-1), k, one, rhs(1,1,j), k) call dgetrs(trans, k, nrhs, D(1,1,j), k, pivot(1,j), * rhs(1,1,j), 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) = rhs(:,:,j) - E(j)*rhs(:,:,j+1) call dgemm(trans, trans, k, nrhs, k, -one, E(1,1,j), k, * rhs(1,1,j+1), k, one, rhs(1,1,j), k) 200 continue else c Solve the transposed system. c First, the right hand side modification. do 300 j = 2,n c Form rhs(:,:,j) = rhs(:,:,j) - E(j-1)^T*rhs(:,:,j-1) call dgemm(trans, 'N', k, nrhs, k, -one, E(1,1,j-1), k, * rhs(1,1,j-1), k, one, rhs(1,1,j), k) 300 continue c Now, the back substitution. c First, solve D(n)^T.rhs(:,:,n) = rhs(:,:,n) call dgetrs(trans, k, nrhs, D(1,1,n), k, pivot(1,n), * rhs(1,1,n), k, info) if ( info .lt. 0 ) then print *,'Illegal parameter number ',-info return endif c For j = n-1 down to 1, solve c D(j)^T.rhs(:,:,j) = rhs(:,:,j) - F(j)^T.rhs(:,:,j+1) do 500 j = n-1,1,-1 c Form rhs(:,:,j) = rhs(:,:,j) - F(j)^T*rhs(:,:,j+1) call dgemm(trans, 'N', k, nrhs, k, -one, F(1,1,j), k, * rhs(1,1,j+1), k, one, rhs(1,1,j), k) call dgetrs(trans, k, nrhs, D(1,1,j), k, pivot(1,j), * rhs(1,1,j), k, info) if ( info .lt. 0 ) then print *,'Illegal parameter number ',-info return endif 500 continue endif c************************************************************************ return end c************************ END OF bl3dsol ********************************