subroutine bl3dsolp(trans, n, k, nk, E, D, F, work, pivot, * nrhs, rhs) c This version of bl3dsolp dated March 31, 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 bl3dfacp. c The right hand sides are stored by with each right hand side c vector consecutive (see below). c************************************************************************ c Uses: c************************************************************************ c Input variables: character *1 trans integer n, k, nk double precision E(k,k,n-1), D(k,k,n), F(k,k,n-1), work(k,k,n-2) integer pivot(nk), nrhs double precision rhs(nk,nrhs) 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 k Integer, the size of the blocks. c n Integer, the number of blocks. c nk Integer, nk = n*k, the size of the system. 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 bl3dpiv. c F(k,k,n-1) Double precision, the super-diagonal blocks, after c modification in bl3dpiv. c work(k,k,n-2) Double precision, work space for the fill-in generated c by bl3dfacp. c pivot(nk) Integer pivot vector returned by bl3dpiv, used in c factoring the diagonal blocks. c pivot((1:k)+(p-1)*k) records the pivoting done in c diagonal block p. c nrhs integer, the number of right hand sides. c rhs(nk,nrhs) Double precision, right hand sides. These are stored c with all elements of one rhs vector consecutive. c For example, for n = 3, k = 2 and nrhs = 3, c the elements of the right hand sides are stored in c 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(kn,nrhs) double precision, the solutions. c c Local variables: integer j, p, nb, col double precision temp c************************************************************************ c Executable statements. c************************************************************************ if ( trans .eq. 'N' .or. trans .eq. 'n') then 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. do 5000 p = 1,nrhs c First, in blocks 1 to n-1. do 1000 nb = 1,n-1 do 800 col = 1,k temp = rhs((nb-1)*k+col,p) rhs((nb-1)*k+col,p) = rhs(pivot((nb-1)*k+col),p) rhs(pivot((nb-1)*k+col),p) = temp do 600 j = col+1,k rhs((nb-1)*k+j,p) = rhs((nb-1)*k+j,p) * + rhs((nb-1)*k+col,p)*D(j,col,nb) 600 continue do 700 j = 1,k rhs(nb*k+j,p) = rhs(nb*k+j,p) * + rhs((nb-1)*k+col,p)*F(j,col,nb) 700 continue 800 continue 1000 continue c Then, in the last block: do 2000 col = 1,k temp = rhs((nb-1)*k+col,p) rhs((nb-1)*k+col,p) = rhs(pivot((nb-1)*k+col),p) rhs(pivot((nb-1)*k+col),p) = temp do 1600 j = col+1,k rhs((nb-1)*k+j,p) = rhs((nb-1)*k+j,p) * + rhs((nb-1)*k+col,p)*D(j,col,nb) 1600 continue 2000 continue c Now, backsubstitution. First, in the last block, involving c only D(n). do 2500 j = k,1,-1 do 2400 col = j+1,k rhs(nk-k+j,p) = rhs(nk-k+j,p) * - D(j,col,n)*rhs(nk-k+col,p) 2400 continue rhs(nk-k+j,p) = rhs(nk-k+j,p)/D(j,j,n) 2500 continue c Now, in the second last block, involving D(n-1) and E(n-1). do 3600 j = k,1,-1 c First, using D(n-1) do 3400 col = j+1,k rhs(nk-2*k+j,p) = rhs(nk-2*k+j,p) * - D(j,col,n-1)*rhs(nk-2*k+col,p) 3400 continue c Then, using E(n-1) do 3500 col = 1,k rhs(nk-2*k+j,p) = rhs(nk-2*k+j,p) * - E(j,col,n-1)*rhs(nk-k+col,p) 3500 continue rhs(nk-2*k+j,p) = rhs(nk-2*k+j,p)/D(j,j,n-1) 3600 continue c Finally, in blocks j = n-2, n-3, ..., 1, which involve c D(j), E(j), G(j). do 4800 nb = n-3,0,-1 do 4700 j = k,1,-1 c First, using D(nb) do 4400 col = j+1,k rhs(nb*k+j,p) = rhs(nb*k+j,p) * - D(j,col,nb+1)*rhs(nb*k+col,p) 4400 continue c Then, using E(nb) do 4500 col = 1,k rhs(nb*k+j,p) = rhs(nb*k+j,p) * - E(j,col,nb+1)*rhs((nb+1)*k+col,p) 4500 continue c Finally, using G(nb) do 4600 col = 1,k rhs(nb*k+j,p) = rhs(nb*k+j,p) * - work(j,col,nb+1)*rhs((nb+2)*k+col,p) 4600 continue rhs(nb*k+j,p) = rhs(nb*k+j,p)/D(j,j,nb+1) 4700 continue 4800 continue 5000 continue else c Solve the transpose problem. do 9000 p = 1,nrhs c First, solve the lower triangular system c D^T(j=1:col,col=1:k,1)rhs(1:k,p) = rhs(1:k,p) do 5100 col = 1,k do 5050 j = 1,col-1 rhs(col,p) = rhs(col,p) - D(j,col,1)*rhs(j,p) 5050 continue rhs(col,p) = rhs(col,p)/D(col,col,1) 5100 continue c Then solve thelower triangular system c D^T(j=1:col,col=1:k,2)rhs(k+1:2k,p) = rhs(k+1:2k,p) c - E^T(1)rhs(1:k,p) do 5300 j = 1,k do 5250 col = 1,k rhs(k+j,p) = rhs(k+j,p) - E(col,j,1)*rhs(col,p) 5250 continue 5300 continue do 5400 col = 1,k do 5350 j = 1,col-1 rhs(k+col,p) = rhs(k+col,p) - D(j,col,2)*rhs(k+j,p) 5350 continue rhs(k+col,p) = rhs(k+col,p)/D(col,col,2) 5400 continue c Then the remaining block rows, solving c D^T(j=1:col,col=1:k,nb)rhs((nb-1)*k+1:nb*k,p) = c rhs((nb-1)*k+1:nb*k,p) - E^T(nb-1)rhs((nb-2)*k+1:(nb-1)*k,p) c - G^T(nb-2)rhs((nb-3)*k+1:(nb-2)*k,p) do 7000 nb = 3,n do 6300 j = 1,k do 6250 col = 1,k rhs((nb-1)*k+j,p) = rhs((nb-1)*k+j,p) * - E(col,j,nb-1)*rhs((nb-2)*k+col,p) * - work(col,j,nb-2)*rhs((nb-3)*k+col,p) 6250 continue 6300 continue do 6400 col = 1,k do 6350 j = 1,col-1 rhs((nb-1)*k+col,p) = rhs((nb-1)*k+col,p) * - D(j,col,nb)*rhs((nb-1)*k+j,p) 6350 continue rhs((nb-1)*k+col,p) = * rhs((nb-1)*k+col,p)/D(col,col,nb) 6400 continue 7000 continue c Now, do the back-substitution. c First, in the bottom block: do 7500 col = k-1,1,-1 do 7400 j = col+1, k rhs((n-1)*k+col,p) = rhs((n-1)*k+col,p) + * D(j,col,n)*rhs((n-1)*k+j,p) 7400 continue if ( pivot((n-1)*k+col) .ne. (n-1)*k+col) then temp = rhs((n-1)*k+col,p) rhs((n-1)*k+col,p) = rhs(pivot((n-1)*k+col),p) rhs(pivot((n-1)*k+col),p) = temp endif 7500 continue c Then, in each main block from block n-1 down to 1: do 8900 nb = n-1,1,-1 do 8500 col = k,1,-1 do 8400 j = col+1,k rhs((nb-1)*k+col,p) = rhs((nb-1)*k+col,p) + * D(j,col,nb)*rhs((nb-1)*k+j,p) 8400 continue do 8450 j = 1,k rhs((nb-1)*k+col,p) = rhs((nb-1)*k+col,p) + * F(j,col,nb)*rhs(nb*k+j,p) 8450 continue if ( pivot((nb-1)*k+col) .ne. (nb-1)*k+col) then temp = rhs((nb-1)*k+col,p) rhs((nb-1)*k+col,p) = rhs(pivot((nb-1)*k+col),p) rhs(pivot((nb-1)*k+col),p) = temp endif 8500 continue 8900 continue 9000 continue endif return end c************************ END OF bl3dsolp *******************************