diff --git a/mat.cc b/mat.cc
index 1c6fb2f..c94c3de 100644
--- a/mat.cc
+++ b/mat.cc
@@ -45,7 +45,7 @@ return r;
 //direct product
 
 template <typename T>
-const NRMat<T> NRMat<T>::otimes(const NRMat<T> &rhs) const
+const NRMat<T> NRMat<T>::otimes(const NRMat<T> &rhs, bool reversecolumns) const
 {
 if(nn==0 && mm == 0) return *this;
 if(rhs.nn==0 && rhs.mm== 0) return rhs;
@@ -54,12 +54,24 @@ NRMat<T> r((T)0,nn*rhs.nn,mm*rhs.mm);
 
 int i,j,k,l;
 
+if(reversecolumns)
+{
+for(i=0;i<nn;i++) for(j=0;j<mm;j++)
+        {
+        T c=(*this)(i,j);
+        for(k=0;k<rhs.mm;k++) for(l=0;l<rhs.mm;l++)
+                r( i*rhs.nn+k , l*nn+j ) = c *rhs(k,l);
+        }
+}
+else
+{
 for(i=0;i<nn;i++) for(j=0;j<mm;j++)
 	{
 	T c=(*this)(i,j);
         for(k=0;k<rhs.mm;k++) for(l=0;l<rhs.mm;l++)
                 r( i*rhs.nn+k , j*rhs.nn+l ) = c *rhs(k,l);
 	}
+}
 
 return r;
 }
diff --git a/mat.h b/mat.h
index 97b132d..5709785 100644
--- a/mat.h
+++ b/mat.h
@@ -59,7 +59,7 @@ public:
 	const T dot(const NRMat &rhs) const; // scalar product of Mat.Mat//@@@for complex do conjugate
 	const NRMat operator*(const NRMat &rhs) const; // Mat * Mat
 	const NRMat oplus(const NRMat &rhs) const; //direct sum
-	const NRMat otimes(const NRMat &rhs) const; //direct product
+	const NRMat otimes(const NRMat &rhs, bool reversecolumns=false) const; //direct product
 	void diagmultl(const NRVec<T> &rhs); //multiply by a diagonal matrix from L
 	void diagmultr(const NRVec<T> &rhs); //multiply by a diagonal matrix from R
 	const NRSMat<T> transposedtimes() const; //A^T . A