Returns in w the LUfactorization (without pivoting) of the banded matrix a of order nrow with (nbandl + 1 + nbandu) bands or diagonals in the work array w .
gauss elimination without pivoting is used. the routine is intended for use with matrices a which do not require row inter changes during factorization, especially for the totally positive matrices which occur in spline calculations. the routine should not be used for an arbitrary banded matrix.
Input
w array of size (nroww,nrow) contains the interesting
part of a banded matrix a , with the diagonals or bands of a
stored in the rows of w , while columns of a correspond to
columns of w . this is the storage mode used in linpack and
results in efficient innermost loops.
explicitly, a has nbandl bands below the diagonal
+ 1 (main) diagonal
+ nbandu bands above the diagonal
and thus, with middle = nbandu + 1,
a(i+j,j) is in w(i+middle,j) for i=nbandu,...,nbandl
j=1,...,nrow .
for example, the interesting entries of a (1,2)banded matrix
of order 9 would appear in the first 1+1+2 = 4 rows of w
as follows.
13 24 35 46 57 68 79
12 23 34 45 56 67 78 89
11 22 33 44 55 66 77 88 99
21 32 43 54 65 76 87 98
all other entries of w not identified in this way with an en
try of a are never referenced .
Output
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(inout),  dimension(nroww,nrow)  ::  w 
work array. See header for details. 

integer(kind=ip),  intent(in)  ::  nroww 
row dimension of the work array w. must be >= nbandl + 1 + nbandu. 

integer(kind=ip),  intent(in)  ::  nrow 
matrix order 

integer(kind=ip),  intent(in)  ::  nbandl 
number of bands of a below the main diagonal 

integer(kind=ip),  intent(in)  ::  nbandu 
number of bands of a above the main diagonal 

integer(kind=ip),  intent(out)  ::  iflag 
indicating success(=1) or failure (=2) 
pure subroutine dbnfac(w,nroww,nrow,nbandl,nbandu,iflag) integer(ip),intent(in) :: nroww !! row dimension of the work array w. must be >= nbandl + 1 + nbandu. integer(ip),intent(in) :: nrow !! matrix order integer(ip),intent(in) :: nbandl !! number of bands of a below the main diagonal integer(ip),intent(in) :: nbandu !! number of bands of a above the main diagonal integer(ip),intent(out) :: iflag !! indicating success(=1) or failure (=2) real(wp),dimension(nroww,nrow),intent(inout) :: w !! work array. See header for details. integer(ip) :: i, ipk, j, jmax, k, kmax, middle, midmk, nrowm1 real(wp) :: factor, pivot iflag = 1_ip middle = nbandu + 1_ip ! w(middle,.) contains the main diagonal of a. nrowm1 = nrow  1_ip if (nrowm1 < 0_ip) then iflag = 2_ip return else if (nrowm1 == 0_ip) then if (w(middle,nrow)==0.0_wp) iflag = 2_ip return end if if (nbandl<=0_ip) then ! a is upper triangular. check that diagonal is nonzero . do i=1_ip,nrowm1 if (w(middle,i)==0.0_wp) then iflag = 2_ip return end if end do if (w(middle,nrow)==0.0_wp) iflag = 2_ip return end if if (nbandu<=0_ip) then ! a is lower triangular. check that diagonal is nonzero and ! divide each column by its diagonal. do i=1_ip,nrowm1 pivot = w(middle,i) if (pivot==0.0_wp) then iflag = 2_ip return end if jmax = min(nbandl,nrowi) do j=1_ip,jmax w(middle+j,i) = w(middle+j,i)/pivot end do end do return end if ! a is not just a triangular matrix. construct lu factorization do i=1_ip,nrowm1 ! w(middle,i) is pivot for ith step . pivot = w(middle,i) if (pivot==0.0_wp) then iflag = 2_ip return end if ! jmax is the number of (nonzero) entries in column i ! below the diagonal. jmax = min(nbandl,nrowi) ! divide each entry in column i below diagonal by pivot. do j=1_ip,jmax w(middle+j,i) = w(middle+j,i)/pivot end do ! kmax is the number of (nonzero) entries in row i to ! the right of the diagonal. kmax = min(nbandu,nrowi) ! subtract a(i,i+k)*(ith column) from (i+k)th column ! (below row i). do k=1_ip,kmax ipk = i + k midmk = middle  k factor = w(midmk,ipk) do j=1_ip,jmax w(midmk+j,ipk) = w(midmk+j,ipk)  w(middle+j,i)*factor end do end do end do ! check the last diagonal entry. if (w(middle,nrow)==0.0_wp) iflag = 2_ip end subroutine dbnfac