/*  
  You will need something like this within your main program

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define min(a, b)  (((a) < (b)) ? (a) : (b))  

#define dim_a 11
#define dim_b 3

  */

void SLVDJI(int IOPT, double B[dim_b+1][dim_a+1], double R[], int NEQ, int 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
*/
   int NRS,IHBP,I,J,K,KK,KC,NN,KI,LIM,IBAND,JC,JI,IBACK,JP,KR,MR;
   double PIVOT, PD, C, SUM;

   NRS=NEQ-1;
   IHBP=IHALFB+1;
   if (IOPT != 2) {
      //
      //  TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
      //
      for (K = 1; K <= NRS; K++) { 
         PIVOT=B[IHBP][K];
         PD = 1./PIVOT;
         KK=K+1;
         KC=IHBP;

         for (I = KK; I <= NEQ; I++) {
         	KC=KC-1;
         	if (KC <= 0) break;
         	C=-B[KC][I]*PD;
         	if (C != 0.0) {
         		B[KC][I]=C;
         		KI=KC+1;
         		LIM=KC+IHALFB;
         		for (J=KI; J<=LIM; J++) {
         			JC=IHBP+J-KC;
         			B[J][I]=B[J][I]+C*B[JC][K];
                        
         		}
         	}
         }
      }
   }
   
   if (IOPT != 1) {
      //
      //  MODIFY LOAD VECTOR R
      //
      NN=NEQ+1;
      IBAND=2*IHALFB+1;
      for (I=2; I<=NEQ; I++) {
         JC=IHBP-I+1;
         JI=1;
         if (JC<=0) {
         	JC=1;
         	JI=I-IHBP+1;
         }
         SUM=0.0;
         for(J=JC; J<=IHALFB; J++) {
         	SUM=SUM+B[J][I]*R[JI];
         	JI=JI+1;
         }
         R[I]=R[I]+SUM;
      }
      //
      //   BACK SOLUTION
      //
      R[NEQ]=R[NEQ]/B[IHBP][NEQ];
      for (IBACK=2; IBACK<=NEQ; IBACK++) {
         I=NN-IBACK;
         JP=I;
         KR=IHBP+1;
         MR=min(IBAND,IHALFB+IBACK);
         SUM=0.0;
         for(J=KR; J<=MR; J++) {
         	JP=JP+1;
         	SUM=SUM+B[J][I]*R[JP];
         }
         R[I]=(R[I]-SUM)/B[IHBP][I];
      }
   }
}


/*
     And here is an example main program that would call the routine
        for testing.
   
void main()
{
   int i;
   
   // Assign the values
   // What you should do is in this section
   //*************************************************************************************
   /*
   k = flag 1, 2 or 3    integer
   B = banded matrix  double reals
   rhs = right hand side vector   double reals
   nn = number of nodes (or equations)   integer
   hbw = half bandwidth    integer
   dim_a = array size for b(__, dim_a) and rhs(dim_a)   integer
   dim_b = array size for b(dim_b, ___)   integer
   */
/*
   int k=3;
   int nn=11;
   int hbw=1;

   // Be sure that the array index begin from 1, NOT 0
   double B[dim_b+1][dim_a+1]={
      {0.0,  0.0,    0.0,    0.0,   0.0,   	 0.0,	  0.0,		 0.0,		0.0,	0.0,    0.0,		0.0,},
      {0.0,  0.0,   2.353,  8.838,  6.508,   7.508,	  8.508,	 5.517,		6.517,	4.533,		5.533, 11.517},
      {0.0,  1.0,  -9.117, -13.307, -11.967,   -13.967,  -11.950,	-9.933,		-8.9,	-7.867,		-14.9, -11.467},
      {0.0,  0.0,  6.838,  4.508,   5.508,   6.508,    3.517,	4.517,		2.533,	3.533,		9.517, 0.0}};

   double rhs[12]={0.0, 5.0, 0.709, 1.058, 2.063, 3.188, 8.156, 16.5, 47.25, 102.0, 110.25, 23.25};
   //*************************************************************************************

   SLVDJI(k,B,rhs,nn,hbw);

   // Output the results
   for (i=1;i<=nn;i++) printf("%13.7f\n",rhs[i]);
}


*/