function F = btrix_3x3(IL,IU,KM,B,C,D,R) % % 3X3 BLOCK Tri-DIAGONAL SOLVER % %USE HERE: IL=2 IU=Jmax-1 KM=KM B=b_x C=c_x D=d_x R=R 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.*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 M=1:3 % 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 % % for I=IS:IU L = I - 1; LL= I - 2; 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('3: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); while ~(I == IU) for M=1:3 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 if I == IU-1, break, 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); end end F = R;