// a class for band matrices

// a band matrix A with a half bandwidth bw >= 0 is a square
// ncol x ncol matrix with all entries A(i,j) = 0 if |i-j| > bw.

// A diagonal matrix is a band matrix with bw = 0.
// A tridiagonal matrix is a band matrix with bw = 1.

#ifndef BANDH
#define BANDH
#include <vector>
using namespace std;

class bandMatrix{
  vector<double> d; // an array for diagonal entries
  vector<double> a; // an array for superdiagonal entries
  vector<double> b; // an array for subdiagonal entries
  int ncol, bw;
 public:
  bandMatrix( int n, int bandwidth );
  double& operator() (int i, int j);
  double operator() (int i, int j) const;
  int cols() const {return ncol;}
  int hbw() const {return bw;}
};
// print small matrices
std::ostream& operator <<( std::ostream& os, const bandMatrix& A );
std::ostream& operator <<( std::ostream& os, const vector<double>& x );

// return A X
vector<double> operator*(const bandMatrix& A, const vector<double>& X );

// return the solution X of A X = Y
vector<double> solve(const bandMatrix& A, const vector<double> & Y );

#endif