%  
function [B, R] = slvdij(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
%                 
%   Assumes B(MaxUnknowns, MaxBandWidth) 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(K,IHBP);
        PD = 1.0/PIVOT; % multiply is somewhat faster than divide
        KK=K+1;
        KC=IHBP;
        for I = KK:NEQ
            KC=KC-1;
            if KC <= 0
                break
            end
            C=-B(I,KC)*PD;
            if C ~= 0.0
                B(I,KC)=C;
                KI=KC+1;
                LIM=KC+IHALFB;
                for J=KI:LIM
                    JC=IHBP+J-KC;
                    B(I,J)=B(I,J)+C*B(K,JC);
                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(I,J)*R(JI);
            JI=JI+1;
        end   % J loop
        R(I)=R(I)+SUM;
    end   % I loop
    %
    %      BACK SOLUTION
    %      
    R(NEQ)=R(NEQ)/B(NEQ,IHBP);
    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(I,J)*R(JP);
        end  % J loop
        R(I)=(R(I)-SUM)/B(I,IHBP);
    end   % IBACK loop
end  %  if IOPT ~= 1 loop