diff --git a/bitvector.cc b/bitvector.cc
index 85361fc..ad3b12e 100644
--- a/bitvector.cc
+++ b/bitvector.cc
@@ -38,6 +38,31 @@ for(unsigned int i=0; i<x.size(); ++i) {x.assign(i,str[i]!='0');}
 return s;
 }
 
+std::ostream & operator<<(std::ostream &s, const bitmatrix &x)
+{
+s<<x.nrows()<<" "<<x.ncols()<<std::endl;
+for(unsigned int i=0; i<x.nrows(); ++i) 
+	{
+	for(unsigned int j=0; j<x.ncols(); ++j) s<< x(i,j);
+	s<<std::endl;
+	}
+return s;
+}
+
+std::istream & operator>>(std::istream  &s, bitmatrix &x)
+{
+unsigned int n,m;
+s >> n >> m;
+x.resize(n,m);
+for(unsigned int i=0; i<n; ++i) 
+        {
+        for(unsigned int j=0; j<m; ++j) {int z; s>>z; if(z) x.set(i,j); else x.reset(i,j);}
+        }
+return s; 
+}
+
+
+
 
 void bitvector::zero_padding() const
 {
diff --git a/bitvector.h b/bitvector.h
index cc838c9..6c759d6 100644
--- a/bitvector.h
+++ b/bitvector.h
@@ -163,6 +163,23 @@ public:
 
 extern bitvector find_irreducible(int deg, int pop= -1, int nth=1); //degree and requested Hamming weight or -1 for random trial
 
+class bitmatrix : public bitvector
+        {
+private:
+        unsigned int nn,mm;
+public:
+	unsigned int nrows() const {return nn;}
+	unsigned int ncols() const {return mm;}
+	bitmatrix() : nn(0),mm(0) {};
+	bitmatrix(unsigned int nn0,unsigned int mm0) : nn(nn0),mm(mm0){bitvector::resize(nn*mm,false);};
+	void set(const unsigned int i, const unsigned int j) {bitvector::set(i*mm+j);};
+	void reset(const unsigned int i, const unsigned int j) {bitvector::reset(i*mm+j);};
+	void flip(const unsigned int i, const unsigned int j) {bitvector::flip(i*mm+j);};
+	const bool assign(const unsigned int i, const unsigned int j, const bool r) {if(r) bitvector::set(i*mm+j); else bitvector::reset(i*mm+j); return r;};
+	const bool operator()(const unsigned int i,  const unsigned int j) const {return bitvector::operator[](i*mm+j);};
+	void resize(const unsigned int n, const unsigned int m) {nn=n; mm=m; bitvector::resize(n*m,false);};
+};
+
 //expand to separate bytes or ints
 template <typename T>
 void bitvector_expand(const bitvector &v, NRVec<T> &r)
@@ -215,6 +232,11 @@ for(int i=0; i<n; ++i) if(v[i]) r.set(i);
 extern std::ostream & operator<<(std::ostream &s, const bitvector &x);
 extern std::istream & operator>>(std::istream  &s, bitvector &x);
 
+extern std::ostream & operator<<(std::ostream &s, const bitmatrix &x);
+extern std::istream & operator>>(std::istream  &s, bitmatrix &x);
+
+
+
 class bitvector_from1 : public bitvector
 	{
 public:
diff --git a/t.cc b/t.cc
index 46fe072..31e4978 100644
--- a/t.cc
+++ b/t.cc
@@ -3412,7 +3412,7 @@ cout <<e;
 cout <<eu;
 }
 
-if(1)
+if(0)
 {
 NRVec<double> a({1.,10.,100.});
 NRVec<double> b({1.,2.,3.});
@@ -3421,4 +3421,23 @@ Tensor<double> bb(b);
 cout << aa*bb;
 }
 
+if(1)
+{
+int seed;
+int f=open("/dev/random",O_RDONLY);
+if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
+close(f);
+srand(seed);
+
+bitmatrix a(4,4);
+a.randomize();
+a.reset(0,0);
+a.set(3,3);
+bitmatrix b(a);
+bitmatrix c=a;
+cout <<a;
+cout <<b;
+cout <<c;
+}
+
 }