%  
function [B, R] = slvdji(IOPT, B, R, NEQ, IHALFB)
%    
%   ASYMMETRIC BAND MATRIX EQUATION SOLVER
%   DOCTORED TO IGNORE ZEROS IN LU DECOMP. STEP
%
%   IOPT=1 TRIANGULARIZES THE BAND MATRIX B
%   IOPT=2 SOLVES FOR RIGHT SIDE R, SOLUTION RETURNS IN R
%   IOPT=3 Performs IOPT = 1 THEN IOPT = 2 before returning
%
%   MODIFICATIONS - DOUBLE PRECISION ON ALL REAL VARIABLES
%                  - REVERSE ROW COLUMN ORDERING - I.E. J,I
%
%   Assumes B(MaxBandWidth, MaxUnknowns) is the Banded Matrix
%           R(MaxUnknowns) is the Right Hand Side Vector
%               (The solution, if IOPT=2 or 3, is returned in R)
%           NEQ is the number of equations to solve
%           IHALFB = Is the half-bandwidth for the specific problem.
%
NRS=NEQ-1;
IHBP=IHALFB+1;
if IOPT ~= 2
    %
    %     TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
    %      
    for K = 1:NRS
        PIVOT=B(IHBP,K);
        PD = 1.0/PIVOT;
        KK=K+1;
        KC=IHBP;
        for I = KK:NEQ
            KC=KC-1;
            if KC <= 0
                break
            end
            C=-B(KC,I)*PD;
            if C ~= 0.0
                B(KC,I)=C;
                KI=KC+1;
                LIM=KC+IHALFB;
                for J=KI:LIM
                    JC=IHBP+J-KC;
                    B(J,I)=B(J,I)+C*B(JC,K);
                end   %  J loop
            end  % if C ~= 0  loop
        end   % I loop
    end	% K loop
end % IOPT ~= 2 condition
if IOPT ~= 1
    %
    %      MODIFY LOAD VECTOR R
    %      
    NN=NEQ+1;
    IBAND=2*IHALFB+1;
    for I=2:NEQ
        JC=IHBP-I+1;
        JI=1;
        if JC <= 0
            JC=1;
            JI=I-IHBP+1;
        end  % if loop
        SUM=0.0;
        for J=JC:IHALFB
            SUM=SUM+B(J,I)*R(JI);
            JI=JI+1;
        end   % J loop
        R(I)=R(I)+SUM;
    end   % I loop
    %
    %      BACK SOLUTION
    %      
    R(NEQ)=R(NEQ)/B(IHBP,NEQ);
    for IBACK=2:NEQ
        I=NN-IBACK;
        JP=I;
        KR=IHBP+1;
        MR=min(IBAND,IHALFB+IBACK);
        SUM=0.0;
        for J=KR:MR
            JP=JP+1;
            SUM=SUM+B(J,I)*R(JP);
        end  % J loop
        R(I)=(R(I)-SUM)/B(IHBP,I);
    end   % IBACK loop
end  %  if IOPT ~= 1 loop