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

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