diff --git a/bfmatrix.cpp b/bfmatrix.cpp
index afcfaa9c31951d743b608a27d719074c7115fdc3..85a8408f74f24598b2d85e9e1d352a015dc5f014 100644
--- a/bfmatrix.cpp
+++ b/bfmatrix.cpp
@@ -42,56 +42,101 @@ void BFMatrix::print(const NEWMAT::Matrix&  m,
 // Member functions for SparseBFMatrix
 //
 
+//
+// Transpose
+//
+
+boost::shared_ptr<BFMatrix> SparseBFMatrix::Transpose()
+const
+{
+  boost::shared_ptr<SparseBFMatrix>   tm(new SparseBFMatrix(mp->t()));
+  return(tm);
+}
+
 //
 // Concatenation of two matrices returning a third
 //
 void SparseBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
 {
+  if (B.Nrows() && Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+
   try {
-    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
-    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
-    if (Nrows() != lB.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
-    *(lAB.mp) = *mp | *(lB.mp);
+    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);      // Means that output is sparse
+    lAB = *this;
+    lAB.HorConcat2MyRight(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error"); 
+    try {
+      FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);        // Means that output is full
+      lAB = FullBFMatrix(this->AsMatrix());
+      lAB.HorConcat2MyRight(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error"); 
+    }
   }
 }
 
 void SparseBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 {
+  if (B.Nrows() && int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+
   try {
-    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
-    if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
-    *(lAB.mp) = *mp | B;
+    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);   // Means that output is sparse
+    lAB = *this;
+    lAB.HorConcat2MyRight(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error"); 
+    try {
+      FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);     // Means that output is full
+      lAB = FullBFMatrix(this->AsMatrix());
+      lAB.HorConcat2MyRight(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error"); 
+    }
   }
 }
 
 void SparseBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
 {
+  if (B.Ncols() && Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+
   try {
-    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
-    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
-    if (Ncols() != lB.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
-    *(lAB.mp) = *mp & *(lB.mp);
+    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);      // Means that output is sparse
+    lAB = *this;
+    lAB.VertConcatBelowMe(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);        // Means that output is full
+      lAB = FullBFMatrix(this->AsMatrix());
+      lAB.VertConcatBelowMe(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    }
   }
 }
 
 void SparseBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 {
+  if (B.Ncols() && int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+
   try {
-    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
-    if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
-    *(lAB.mp) = *mp & B;
+    SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);      // Means that output is sparse
+    lAB = *this;
+    lAB.VertConcatBelowMe(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);        // Means that output is full
+      lAB = FullBFMatrix(this->AsMatrix());
+      lAB.VertConcatBelowMe(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    }
   }
 }
 
@@ -100,36 +145,58 @@ void SparseBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 //
 void SparseBFMatrix::HorConcat2MyRight(const BFMatrix& B)
 {
+  if (!B.Nrows()) return;
+
+  if (Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+
   try {
-    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
-    if (Nrows() != lB.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);   // Means that we want to concatenate a sparse matrix
     *mp |= *(lB.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);     // Means that we want to concatenate a full
+      this->HorConcat2MyRight(lB.AsMatrix());
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: dynamic cast error");
+    } 
   }
 }
 
 void SparseBFMatrix::HorConcat2MyRight(const NEWMAT::Matrix& B)
 {
+  if (!B.Nrows()) return;
+
   if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
   *mp |= B;
 }
 
 void SparseBFMatrix::VertConcatBelowMe(const BFMatrix& B)
 {
+  if (!B.Ncols()) return;
+
+  if (Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+
   try {
-    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
-    if (Ncols() != lB.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+    const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);   // Means that we want to concatenate a sparse matrix
     *mp &= *(lB.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);     // Means that we want to concatenate a full
+      this->VertConcatBelowMe(lB.AsMatrix());
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: dynamic cast error");
+    } 
   }
 }
 
 void SparseBFMatrix::VertConcatBelowMe(const NEWMAT::Matrix& B)
 {
+  if (!B.Ncols()) return;
+
   if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
   *mp &= B;
 }
@@ -146,16 +213,24 @@ NEWMAT::ReturnMatrix SparseBFMatrix::MulByVec(const NEWMAT::ColumnVector& invec)
 // Add another matrix to this one
 void SparseBFMatrix::AddToMe(const BFMatrix& M, double s)
 {
+  if (Ncols() != M.Ncols() || Nrows() != M.Nrows()) {
+    throw BFMatrixException("SparseBFMatrix::AddToMe: Matrix size mismatch");
+  }
+
   try {
-    const SparseBFMatrix& lM = dynamic_cast<const SparseBFMatrix&>(M);
-    if (Ncols() != lM.Ncols() || Nrows() != lM.Nrows()) {
-      throw BFMatrixException("SparseBFMatrix::AddToMe: Matrix size mismatch");
-    }
+    const SparseBFMatrix& lM = dynamic_cast<const SparseBFMatrix&>(M);  // Add sparse matrix to this sparse matrix
     if (s == 1.0) *mp += *(lM.mp);
     else *mp += s * *(lM.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("SparseBFMatrix::AddToMe: dynamic cast error"); 
+    try {
+      const FullBFMatrix& lM = dynamic_cast<const FullBFMatrix&>(M);      // Add full matrix to this sparse matrix
+      if (s == 1.0) *mp += SpMat<double>(lM.ReadAsMatrix());
+      else *mp += s * SpMat<double>(lM.ReadAsMatrix());
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("SparseBFMatrix::AddToMe: dynamic cast error"); 
+    }
   }
 }
 
@@ -177,11 +252,11 @@ NEWMAT::ReturnMatrix SparseBFMatrix::SolveForx(const NEWMAT::ColumnVector& b,
 // Member functions for FullBFMatrix
 //
 
-boost::shared_ptr<BFMatrix> FullBFMatrix::Transpose(boost::shared_ptr<BFMatrix>& pA) const
+boost::shared_ptr<BFMatrix> FullBFMatrix::Transpose() 
+const
 {
-  boost::shared_ptr<FullBFMatrix>  tmp(new FullBFMatrix);
-  *(tmp->mp) = mp->t();
-  return(tmp);
+  boost::shared_ptr<FullBFMatrix>  tm(new FullBFMatrix(mp->t()));
+  return(tm);
 }
 
 //
@@ -190,51 +265,85 @@ boost::shared_ptr<BFMatrix> FullBFMatrix::Transpose(boost::shared_ptr<BFMatrix>&
 
 void FullBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
 {
+  if (B.Nrows() && Nrows() != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
+
   try {
-    const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
-    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
-    if (Nrows() != lB.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
-    *(lAB.mp) = *mp | *(lB.mp);
+    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);  // This means output is a full matrix
+    lAB = *this;
+    lAB.HorConcat2MyRight(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error"); 
+    try {
+      SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+      lAB = SparseBFMatrix(this->AsMatrix());
+      lAB.HorConcat2MyRight(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error"); 
+    }
   }
 }
 
 void FullBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 {
+  if (B.Nrows() && int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
+
   try {
-    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
-    if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
-    *(lAB.mp) = *mp | B;
+    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);  // This means output is a full matrix
+    lAB = *this;
+    lAB.HorConcat2MyRight(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error"); 
+    try {
+      SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+      lAB = SparseBFMatrix(this->AsMatrix());
+      lAB.HorConcat2MyRight(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error"); 
+    }
   }
 }
 
 void FullBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
 {
+  if (B.Ncols() && Ncols() != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
+
   try {
-    const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
-    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
-    if (Ncols() != lB.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
-    *(lAB.mp) = *mp & *(lB.mp);
+    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);  // This means output is a full matrix
+    lAB = *this;
+    lAB.VertConcatBelowMe(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+      lAB = SparseBFMatrix(this->AsMatrix());
+      lAB.VertConcatBelowMe(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error"); 
+    }
   }
 }
 
 void FullBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 {
+  if (B.Ncols() && int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
+
   try {
-    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
-    if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
-    *(lAB.mp) = *mp & B;
+    FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);  // This means output is a full matrix
+    lAB = *this;
+    lAB.VertConcatBelowMe(B);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error"); 
+    try {
+      SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+      lAB = SparseBFMatrix(this->AsMatrix());
+      lAB.VertConcatBelowMe(B);
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error"); 
+    }
   }
 }
 
@@ -243,36 +352,58 @@ void FullBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
 //
 void FullBFMatrix::HorConcat2MyRight(const BFMatrix& B)
 {
+  if (!B.Nrows()) return;
+
+  if (Nrows() != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+
   try {
     const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
-    if (Nrows() != lB.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
     *mp |= *(lB.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: dynamic cast error"); 
+    try {
+      const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+      this->HorConcat2MyRight(lB.AsMatrix());
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: dynamic cast error"); 
+    }
   }
 }
 
 void FullBFMatrix::HorConcat2MyRight(const NEWMAT::Matrix& B)
 {
+  if (!B.Nrows()) return;
+
   if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
   *mp |= B;
 }
 
 void FullBFMatrix::VertConcatBelowMe(const BFMatrix& B)
 {
+  if (!B.Ncols()) return;
+
+  if (Ncols() != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+
   try {
     const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
-    if (Ncols() != lB.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
     *mp &= *(lB.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::HorConcatBelowMe: dynamic cast error"); 
+    try {
+      const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+      this->VertConcatBelowMe(lB.AsMatrix());
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::HorConcatBelowMe: dynamic cast error"); 
+    }
   }
 }
 
 void FullBFMatrix::VertConcatBelowMe(const NEWMAT::Matrix& B)
 {
+  if (!B.Ncols()) return;
+
   if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
   *mp &= B;
 }
@@ -297,15 +428,22 @@ NEWMAT::ReturnMatrix FullBFMatrix::MulByVec(const NEWMAT::ColumnVector& invec) c
 // Add another matrix to this one
 void FullBFMatrix::AddToMe(const BFMatrix& m, double s)
 {
+  if (Ncols() != m.Ncols() || Nrows() != m.Nrows()) {
+    throw BFMatrixException("FullBFMatrix::AddToMe: Matrix size mismatch");
+  }
+
   try {
     const FullBFMatrix& lm = dynamic_cast<const FullBFMatrix&>(m);
-    if (Ncols() != lm.Ncols() || Nrows() != lm.Nrows()) {
-      throw BFMatrixException("FullBFMatrix::AddToMe: Matrix size mismatch");
-    }
     *mp += s*(*lm.mp);
   }
   catch (std::bad_cast) {
-    throw BFMatrixException("FullBFMatrix::AddToMe: dynamic cast error"); 
+    try {
+      const SparseBFMatrix& lm = dynamic_cast<const SparseBFMatrix&>(m);
+      *mp += s*lm.AsMatrix();
+    }
+    catch (std::bad_cast) {
+      throw BFMatrixException("FullBFMatrix::AddToMe: dynamic cast error");
+    }
   }
 }
 
diff --git a/bfmatrix.h b/bfmatrix.h
index 87fe792f743357e035b4ccad9112b3dedde285d0..847701d80d2c59248c238795a24d67c71c579983 100644
--- a/bfmatrix.h
+++ b/bfmatrix.h
@@ -85,7 +85,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) = 0;  
 
   // Transpose
-  virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const = 0;
+  virtual boost::shared_ptr<BFMatrix> Transpose() const = 0;
 
   // Concatenation. Note that desired polymorphism prevents us from using BFMatrix->NEWMAT::Matrix conversion
   // Concatenate two matrices yielding a third
@@ -162,7 +162,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) {mp->AddTo(x,y,val);}
 
   // Transpose. 
-  virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const {throw BFMatrixException("SparseBFMatrix::Transpose: Not yet implemented");}
+  virtual boost::shared_ptr<BFMatrix> Transpose() const;
   
   // Concatenation of two matrices returning a third
   // AB = [*this B] in Matlab lingo
@@ -210,6 +210,7 @@ public:
   }
 
   virtual NEWMAT::ReturnMatrix AsMatrix() const {NEWMAT::Matrix ret; ret = *mp; ret.Release(); return(ret);}
+  virtual const NEWMAT::Matrix& ReadAsMatrix() const {return(*mp);} 
 
   // Basic properties
   virtual unsigned int Nrows() const {return(mp->Nrows());}
@@ -234,7 +235,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) {(*mp)(x,y)+=val;}
 
   // Transpose. 
-  virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const;
+  virtual boost::shared_ptr<BFMatrix> Transpose() const;
   
   // Concatenation of two matrices returning a third
   virtual void HorConcat(const BFMatrix& B, BFMatrix& AB) const;