diff --git a/matrix/dw_matrix.c b/matrix/dw_matrix.c
index fb978baaad3403511e92fb7b955e83010f04b369..128fd7e34028c55bde55406f7931ecd274ea76af 100644
--- a/matrix/dw_matrix.c
+++ b/matrix/dw_matrix.c
@@ -1558,6 +1558,54 @@ TMatrix UpdateProductMM(PRECISION a, TMatrix X, PRECISION b, TMatrix Y, TMatrix
     bProductMM_Update(pElementM(X),pElementM(Y),pElementM(Z),a,b,RowM(X),ColM(X),ColM(Y),MajorForm(X),MajorForm(Y),MajorForm(Z));
   return X;
 }
+
+/*
+   Assumes
+     a : scalar
+     X : m x n matrix
+     b : scalar
+     Y : m x k matrix
+     Z : k x n matrix
+
+   Results
+     X = a*X + b*Y'*Z
+
+   Returns
+     Returns X upon success and null on failure.  Call GetError() to
+     determine the cause of failure.
+
+   Notes
+     X, Y, and Z do not have to be distinct matrices.
+*/
+TMatrix UpdateTransposeProductMM(PRECISION a, TMatrix X, PRECISION b, TMatrix Y, TMatrix Z)
+{
+  PRECISION *ptr;
+  if (!X || !Y || ! Z)
+    {
+      dw_Error(NULL_ERR);
+      return (TMatrix)NULL;
+    }
+  if ((RowM(X) != RowM(Y)) || (ColM(X) != ColM(Z)) || (ColM(Y) != RowM(Z)))
+    {
+      dw_Error(SIZE_ERR);
+      return (TMatrix)NULL;
+    }
+  if ((X == Y) || (X == Z))
+    {
+      if (!(ptr=(PRECISION*)dw_malloc(RowM(X)*ColM(X)*sizeof(PRECISION))))
+	{
+	  dw_Error(MEM_ERR); 
+	  return (TMatrix)NULL; 
+	}
+      memcpy(ptr,pElementM(X),RowM(X)*ColM(X)*sizeof(PRECISION));
+      bProductMM_Update(pElementM(X),(X == Y) ? ptr : pElementM(Y),(X == Z) ? ptr : pElementM(Z),a,b,RowM(X),
+			    ColM(X),ColM(Y),MajorForm(X),1^MajorForm(Y),MajorForm(Z));
+      dw_free(ptr);
+    }
+  else
+    bProductMM_Update(pElementM(X),pElementM(Y),pElementM(Z),a,b,RowM(X),ColM(X),ColM(Y),MajorForm(X),1^MajorForm(Y),MajorForm(Z));
+  return X;
+}
 /*******************************************************************************/
 /*******************************************************************************/
 /*******************************************************************************/