function F = bpentax_3x3(IL,IU,KM,A,B,C,D,E,R)
%
%            3X3 BLOCK PENTADIAGONAL SOLVER
IS = IL + 2;
I  = IL;
ierr = 0;
ior  = 0;
%
for n = 1:3
  for m = 1:3
    G(KM,n,m) = C(I,KM,n,m);
  end
end
[L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
if ~(ierr==0)
  fprintf('1:Failure of LU decomposition ierr = %g',ierr);
end

D1 = V1.*R(I,KM,1);
D2 = V2.*( R(I,KM,2) - L21.*D1);
D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
R(I,KM,3) = D3;
R(I,KM,2) = D2 - U23.*R(I,KM,3);
R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
%
for M=1:3
  D1 = V1.*E(I,KM,1,M);
  D2 = V2.*(E(I,KM,2,M) - L21.*D1);
  D3 = V3.*(E(I,KM,3,M) - L31.*D1 - L32.*D2);
  V(I,KM,3,M) = D3;
  V(I,KM,2,M) = D2 - U23.*V(I,KM,3,M);
  V(I,KM,1,M) = D1 - U13.*V(I,KM,3,M) - U12.*V(I,KM,2,M);
  D1 = V1.*D(I,KM,1,M);
  D2 = V2.*(D(I,KM,2,M) - L21.*D1);
  D3 = V3.*(D(I,KM,3,M) - L31.*D1 - L32.*D2);
  U(I,KM,3,M) = D3;
  U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
  U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
end
%
%
I = IL + 1;
L = I  - 1;
%
%
for N=1:3
  for M=1:3
    G(KM,N,M) = C(I,KM,N,M) - B(I,KM,N,1).*U(L,KM,1,M) ...
	- B(I,KM,N,2).*U(L,KM,2,M)                  ...
	- B(I,KM,N,3).*U(L,KM,3,M);
  end
end
[L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
if ~(ierr==0)
  fprintf('2:Failure of LU decomposition ierr = %g',ierr);
end
for M=1:3
  R(I,KM,M) = R(I,KM,M) - B(I,KM,M,1).*R(L,KM,1)  ...
      - B(I,KM,M,2).*R(L,KM,2)                 ...
      - B(I,KM,M,3).*R(L,KM,3);
end
D1 = V1.*R(I,KM,1);
D2 = V2.*( R(I,KM,2) - L21.*D1);
D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
R(I,KM,3) = D3;
R(I,KM,2) = D2 - U23.*R(I,KM,3);
R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
for N=1:3
  for M=1:3
    H(N,KM,M) = D(I,KM,N,M) - B(I,KM,N,1).*V(L,KM,1,M) ...
	- B(I,KM,N,2).*V(L,KM,2,M)                  ...
	- B(I,KM,N,3).*V(L,KM,3,M);
  end
end
for M=1:3
  D1 = V1.*E(I,KM,1,M);
  D2 = V2.*(E(I,KM,2,M) - L21.*D1);
  D3 = V3.*(E(I,KM,3,M) - L31.*D1 - L32.*D2);
  V(I,KM,3,M) = D3;
  V(I,KM,2,M) = D2 - U23.*V(I,KM,3,M);
  V(I,KM,1,M) = D1 - U13.*V(I,KM,3,M) - U12.*V(I,KM,2,M);
  %     
  D1 = V1.*H(1,KM,M);
  D2 = V2.*(H(2,KM,M) - L21.*D1);
  D3 = V3.*(H(3,KM,M) - L31.*D1 - L32.*D2);
  U(I,KM,3,M) = D3;
  U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
  U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
end
%
%  
for I=IS:IU
  L = I - 1;
  LL= I - 2;
  for N=1:3
    for M=1:3
      BH(I,KM,N,M) = B(I,KM,N,M)                    ...
	  - A(I,KM,N,1).*U(LL,KM,1,M)  ...
	  - A(I,KM,N,2).*U(LL,KM,2,M)               ...
	  - A(I,KM,N,3).*U(LL,KM,3,M);
    end
  end
  for N=1:3
    for M=1:3
      G(KM,N,M) = C(I,KM,N,M) - A(I,KM,N,1).*V(LL,KM,1,M)  ...
	  - A(I,KM,N,2).*V(LL,KM,2,M)                     ...
	  - A(I,KM,N,3).*V(LL,KM,3,M)                     ...
	  - BH(I,KM,N,1).*U(L,KM,1,M)                     ...
	  - BH(I,KM,N,2).*U(L,KM,2,M)                     ...
	  - BH(I,KM,N,3).*U(L,KM,3,M);
    end
  end
  [L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
  if ~(ierr==0)
    fprintf('3:Failure of LU decomposition ierr = %g',ierr);
  end
  for M=1:3
    R(I,KM,M) = R(I,KM,M) - BH(I,KM,M,1).*R(L,KM,1)   ...
	- BH(I,KM,M,2).*R(L,KM,2)                  ...
	- BH(I,KM,M,3).*R(L,KM,3)                  ...
	- A(I,KM,M,1).*R(LL,KM,1)                  ...
	- A(I,KM,M,2).*R(LL,KM,2)                  ...
	- A(I,KM,M,3).*R(LL,KM,3);
  end
  D1 = V1.*R(I,KM,1);
  D2 = V2.*( R(I,KM,2) - L21.*D1);
  D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
  R(I,KM,3) = D3;
  R(I,KM,2) = D2 - U23.*R(I,KM,3);
  R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
  while ~(I == IU)
    for N=1:3
      for M=1:3
	H(N,KM,M) = D(I,KM,N,M) - BH(I,KM,N,1).*V(L,KM,1,M)  ...
	    - BH(I,KM,N,2).*V(L,KM,2,M)                     ...
	    - BH(I,KM,N,3).*V(L,KM,3,M);
      end
    end
    for M=1:3
      D1 = V1.*H(1,KM,M);
      D2 = V2.*(H(2,KM,M) - L21.*D1);
      D3 = V3.*(H(3,KM,M) - L31.*D1 - L32.*D2);
      U(I,KM,3,M) = D3;
      U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
      U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
    end
    if I == IU-1, break, end
    for M=1:3
      D1 = V1.*E(I,KM,1,M);
      D2 = V2.*(E(I,KM,2,M) - L21.*D1);
      D3 = V3.*(E(I,KM,3,M) - L31.*D1 - L32.*D2);
      V(I,KM,3,M) = D3;
      V(I,KM,2,M) = D2 - U23.*V(I,KM,3,M);
      V(I,KM,1,M) = D1 - U13.*V(I,KM,3,M) - U12.*V(I,KM,2,M);
    end
    I = IU;
  end
end
%
%    BACK SWEEP
%
I = IU - 1;
for M=1:3
  R(I,KM,M) = R(I,KM,M) - U(I,KM,M,1).*R(I+1,KM,1)         ...
      - U(I,KM,M,2).*R(I+1,KM,2) - U(I,KM,M,3).*R(I+1,KM,3);
end
for I = IU-2:-1:IL
  for M=1:3
    R(I,KM,M) = R(I,KM,M) - U(I,KM,M,1).*R(I+1,KM,1)      ...
	- U(I,KM,M,2).*R(I+1,KM,2)                             ...
	- U(I,KM,M,3).*R(I+1,KM,3)                             ...
	- V(I,KM,M,1).*R(I+2,KM,1)                         ...
	- V(I,KM,M,2).*R(I+2,KM,2)                         ...
	- V(I,KM,M,3).*R(I+2,KM,3);
  end
end
F = R;