subroutine bl3dfacp(n, k, E, D, F, work, pivot, ifail) c This version of bl3dfacp dated March 29, 2000. c************************************************************************ c bl3dfacp 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 among the rows of D(j), F(j), resulting in c fill in. c************************************************************************ c************************************************************************ c Input variables: integer n, k double precision E(k,k,n-1), D(k,k,n), F(k,k,n-1), work(k,k,n-2) integer pivot(*), ifail 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 work(k,k,n-2) Double precision, the fill in diagonal. Need not be c initialized, and may be declared as work(M) where c M >= k*k*(n-2). c Output: c The factors L,U are in the lower and upper triangles, respectively. c pivot(k*n) Integer pivot vector used in the elimination. c ifail = 0, if all is well. c j, 1 <= j <= n, if a zero pivot is encountered at the c jth stage. c Local variables: c Uses: c idamax, daxpy, dscal from the Blas. integer col, nb, j, index1, index2 double precision one, scale, zero parameter ( one = 1.0d0, zero = 0.0d0 ) intrinsic abs integer idamax c************************************************************************ c Executable statements. c************************************************************************ ifail = 0 c Initialize work: do 500 nb = 1,n-2 do 500 col = 1,k do 500 j = 1,k 500 work(col,j,nb) = zero c First, in main blocks 1 to n-2: c************************************************************************ do 1000 nb = 1,n-2 c In each block of k rows: do 900 col = 1,k c Find index of maximum absolute value in D(col:k,col,nb) index1 = idamax(k+1-col, D(col,col, nb), 1) + col - 1 c Find index of maximum absolute value in F(1:k,col,nb) index2 = idamax(k, F(1,col, nb), 1) if (abs(D(index1,col,nb)).ge. abs(F(index2,col,nb)) ) then if (D(index1,col,nb) .eq. zero ) then ifail = (nb-1)*k+col return endif pivot((nb-1)*k+col) = (nb-1)*k+index1 c Swap row "col", and row "index1" in D(nb),E(nb), work(nb). c Do not swap the stored multipliers in the lower triangle c of D(nb). call dswap(k+1-col, D(col,col,nb), k, * D(index1,col,nb), k) call dswap(k, E(col,1,nb), k, E(index1,1,nb), k) call dswap(k, work(col,1,nb), k, work(index1,1,nb), k) else if (E(index2,col,nb) .eq. zero ) then ifail = (nb-1)*k+col return endif pivot((nb-1)*k+col) = nb*k+index2 c Swap row "col" in D(nb), and row "index2" in F(nb), c and also in E(nb), D(nb+1). and in work(nb) and E(nb+1). c Do not swap the multipliers. call dswap(k+1-col, D(col,col,nb), k, * F(index2,col,nb), k) call dswap(k, E(col,1,nb), k, D(index2,1,nb+1), k) call dswap(k, work(col,1,nb), k, E(index2,1,nb+1), k) endif scale = -one/D(col,col,nb) call dscal(k-col, scale, D(col+1,col,nb), 1) call dscal(k, scale, F(1,col,nb), 1) c Subtract the appropriate multiples of row col in D(nb), c from the diagonal across, from rows col+1 to k of D(nb), c and from rows 1 to k of F(nb). Subtract the same multiples c of the full row col in E(nb) from rows col+1 to k of E(nb) c and from rows 1 to k of D(nb+1). Do the same for work(nb) and c E(nb+1). do 700 j = col+1,k call daxpy(k-col, D(j,col,nb), D(col,col+1,nb), k, * D(j,col+1,nb),k) call daxpy(k, D(j,col,nb), E(col,1,nb), k, * E(j,1,nb),k) call daxpy(k, D(j,col,nb), work(col,1,nb), k, * work(j,1,nb),k) 700 continue do 800 j = 1,k call daxpy(k-col, F(j,col,nb), D(col,col+1,nb), k, * F(j,col+1,nb),k) call daxpy(k, F(j,col,nb), E(col,1,nb), k, * D(j,1,nb+1),k) call daxpy(k, F(j,col,nb), work(col,1,nb), k, * E(j,1,nb+1),k) 800 continue 900 continue 1000 continue c Now, the operations in block n-1. do 2000 col = 1,k c Find index of maximum absolute value in D(col:k,col,n-1) index1 = idamax(k+1-col, D(col,col, n-1), 1) + col - 1 c Find index of maximum absolute value in F(1:k,col,n-1) index2 = idamax(k, F(1,col, n-1), 1) if (abs(D(index1,col,n-1)).ge. abs(F(index2,col,n-1)) ) then if (D(index1,col,n-1) .eq. zero ) then ifail = (n-2)*k+col return endif pivot((n-2)*k+col) = (n-2)*k+index1 c Swap row "col", and row "index1" in D(n-1),E(n-1). c Do not swap the stored multipliers in the lower triangle c of D(n-1). call dswap(k+1-col, D(col,col,n-1), k, * D(index1,col,n-1), k) call dswap(k, E(col,1,n-1), k, E(index1,1,n-1), k) else if (E(index2,col,n-1) .eq. zero ) then ifail = (n-2)*k+col return endif pivot((n-2)*k+col) = (n-1)*k+index2 c Swap row "col" in D(n-1), and row "index2" in F(n-1), c and also in E(n-1), D(n). c Do not swap the multipliers. call dswap(k+1-col, D(col,col,n-1), k, * F(index2,col,n-1), k) call dswap(k, E(col,1,n-1), k, D(index2,1,n-1+1), k) endif scale = -one/D(col,col,n-1) call dscal(k-col, scale, D(col+1,col,n-1), 1) call dscal(k, scale, F(1,col,n-1), 1) c Subtract the appropriate multiples of row col in D(n-1), c from the diagonal across, from rows col+1 to k of D(n-1), c and from rows 1 to k of F(n-1). Subtract the same multiples c of the full row col in E(n-1) from rows col+1 to k of E(n-1) c and from rows 1 to k of D(n). do 1700 j = col+1,k call daxpy(k-col, D(j,col,n-1), D(col,col+1,n-1), k, * D(j,col+1,n-1),k) call daxpy(k, D(j,col,n-1), E(col,1,n-1), k, * E(j,1,n-1),k) 1700 continue do 1800 j = 1,k call daxpy(k-col, F(j,col,n-1), D(col,col+1,n-1), k, * F(j,col+1,n-1),k) call daxpy(k, F(j,col,n-1), E(col,1,n-1), k, * D(j,1,n),k) 1800 continue 2000 continue c Finally, reduce the last block, D(n), to upper triangular form. do 3000 col = 1,k-1 c Find index of maximum absolute value in D(col:k,col,n) index1 = idamax(k+1-col, D(col,col, n), 1) + col - 1 if (D(index1,col,n) .eq. zero ) then ifail = (n-1)*k+col return endif pivot((n-1)*k+col) = (n-1)*k+index1 c Swap row "col", and row "index1" in D(n). c Do not swap the stored multipliers in the lower triangle c of D(n). call dswap(k+1-col, D(col,col,n), k, * D(index1,col,n), k) scale = -one/D(col,col,n) call dscal(k-col, scale, D(col+1,col,n), 1) c Subtract the appropriate multiples of row col in D(n), c from the diagonal across, from rows col+1 to k of D(n), do 2700 j = col+1,k call daxpy(k-col, D(j,col,n), D(col,col+1,n), k, * D(j,col+1,n),k) 2700 continue 3000 continue return end c************************ END OF bl3dfacp *******************************