diff --git a/PETE/QDPOps.in b/PETE/QDPOps.in
index c5e09af18..5a07a8e22 100644
--- a/PETE/QDPOps.in
+++ b/PETE/QDPOps.in
@@ -373,6 +373,16 @@ binaryOps
   @return (real(trace(adj(l)*r))
   @ingroup group1
   @relates QDPType */"
+-----
+  TAG = "FnLocalColorInnerProduct"
+  FUNCTION = "localColorInnerProduct"
+  EXPR = "return (localColorInnerProduct(a,b));"
+  COMMENTA = "//! InnerProduct on only color fiber indices"
+  COMMENTB = "/*! L2 norm only color fiber indices
+  @sa adj(), trace()
+  @return (traceColor(adj(l)*r)
+  @ingroup group1
+  @relates QDPType */"
 -----
   TAG = "FnQuarkContract13"
   FUNCTION = "quarkContract13"
diff --git a/examples/t_basic.cc b/examples/t_basic.cc
index ee216dff0..392c3c179 100644
--- a/examples/t_basic.cc
+++ b/examples/t_basic.cc
@@ -46,6 +46,8 @@ int main(int argc, char *argv[])
   LatticeFermion    lftmp1;
   LatticeFermion    lftmp2;
   LatticeFermion    lftmp3;
+  LatticeColorVectorSpinMatrix lcvsp1;
+  LatticeColorVectorSpinMatrix lcvsp2;
   LatticeReal R1;
   LatticeReal Rtmp;
   LatticeComplex c;
@@ -147,7 +149,7 @@ int main(int argc, char *argv[])
   pop(xml_out);
 
 
-#if 0
+#if 1
 
   mu = 0;
   nu = 1;
@@ -181,7 +183,7 @@ int main(int argc, char *argv[])
   pop(xml_out);
 
   /* test 5 */
-  rtmp = real(trace(lctmp1 * u));
+  LatticeReal rtmp = real(trace(lctmp1 * u));
   push(xml_out,"TRACE_MULTIPLY_MATRIX_realpart");
   write(xml_out, "rtmp", rtmp);
   pop(xml_out);
@@ -391,7 +393,7 @@ int main(int argc, char *argv[])
 
   /* test 27 */
   gaussian(lqtmp1);
-  r = real(trace(adj(q) * lqtmp1));
+  LatticeReal r = real(trace(adj(q) * lqtmp1));
   push(xml_out,"TRACE_MULT_PROP");
   write(xml_out, "r", r);
   pop(xml_out);
@@ -406,6 +408,33 @@ int main(int argc, char *argv[])
   write(xml_out, "lqtmp24", quarkContract24(lqtmp1, lqtmp2));
   pop(xml_out);
 
+#if 1
+  /* test 29 */
+  gaussian(lcvsp1);
+  gaussian(lcvsp2);
+  //LatticeSpinMatrix spm = localColorInnerProduct(lcvsp1, lcvsp2);
+  auto spm = localColorInnerProduct(lcvsp1, lcvsp2);
+  push(xml_out,"LOCAL_COLOR_INNER_PRODUCT");
+  write(xml_out, "spm", spm);
+  pop(xml_out);
+
+  /* test 30 */
+  gaussian(lcvsp1);
+  gaussian(lcvsp2);
+  multi1d<SpinMatrixD> fred1 = sumMulti(spm, QDP::rb);
+  push(xml_out,"SUMMULTI_LOCAL_COLOR_INNER_PRODUCT1");
+  write(xml_out, "fred1", fred1);
+  pop(xml_out);
+
+  /* test 31 */
+  gaussian(lcvsp1);
+  gaussian(lcvsp2);
+  multi1d<SpinMatrixD> fred2 = sumMulti(localColorInnerProduct(lcvsp1, lcvsp2), QDP::rb);
+  push(xml_out,"SUMMULTI_LOCAL_COLOR_INNER_PRODUCT2");
+  write(xml_out, "fred2", fred2);
+  pop(xml_out);
+#endif
+
 #endif
 
 
diff --git a/include/QDPOperators.h b/include/QDPOperators.h
index 766d438b7..7f21dcb38 100644
--- a/include/QDPOperators.h
+++ b/include/QDPOperators.h
@@ -1,4 +1,5 @@
 // -*- C++ -*-
+// $Id: Header.h,v 1.3 2002-10-14 02:06:56 edwards Exp $
 
 /*! @file
  * @brief Bulk of QDP operators produced by PETE
@@ -461,6 +462,17 @@ struct FnLocalInnerProductReal
   }
 };
 
+struct FnLocalColorInnerProduct
+{
+  PETE_EMPTY_CONSTRUCTORS(FnLocalColorInnerProduct)
+  template<class T1, class T2>
+  inline typename BinaryReturn<T1, T2, FnLocalColorInnerProduct >::Type_t
+  operator()(const T1 &a, const T2 &b) const
+  {
+    return (localColorInnerProduct(a,b));
+  }
+};
+
 struct FnQuarkContract13
 {
   PETE_EMPTY_CONSTRUCTORS(FnQuarkContract13)
@@ -1635,6 +1647,28 @@ localInnerProductReal(const QDPType<T1,C1> & l,const QDPType<T2,C2> & r)
     CreateLeaf<QDPType<T2,C2> >::make(r)));
 }
 
+//! InnerProduct on only color fiber indices
+/*! L2 norm only color fiber indices
+  @sa adj(), trace()
+  @return (traceColor(adj(l)*r)
+  @ingroup group1
+  @relates QDPType */
+template<class T1,class C1,class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+  typename CreateLeaf<QDPType<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPType<T1,C1> & l,const QDPType<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+    typename CreateLeaf<QDPType<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPType<T1,C1> >::make(l),
+    CreateLeaf<QDPType<T2,C2> >::make(r)));
+}
+
 //! Contraction for quark propagators
 /*! 
    epsilon contract 2 quark propagators and return a quark propagator.
@@ -2401,6 +2435,22 @@ localInnerProductReal(const QDPType<T1,C1> & l,const QDPExpr<T2,C2> & r)
     CreateLeaf<QDPExpr<T2,C2> >::make(r)));
 }
 
+template<class T1,class C1,class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+  typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPType<T1,C1> & l,const QDPExpr<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+    typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPType<T1,C1> >::make(l),
+    CreateLeaf<QDPExpr<T2,C2> >::make(r)));
+}
+
 template<class T1,class C1,class T2,class C2>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
@@ -2945,6 +2995,22 @@ localInnerProductReal(const QDPExpr<T1,C1> & l,const QDPType<T2,C2> & r)
     CreateLeaf<QDPType<T2,C2> >::make(r)));
 }
 
+template<class T1,class C1,class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+  typename CreateLeaf<QDPType<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPExpr<T1,C1> & l,const QDPType<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+    typename CreateLeaf<QDPType<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPExpr<T1,C1> >::make(l),
+    CreateLeaf<QDPType<T2,C2> >::make(r)));
+}
+
 template<class T1,class C1,class T2,class C2>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
@@ -3489,6 +3555,22 @@ localInnerProductReal(const QDPType<T1,C1> & l,const typename WordType<C1>::Type
     CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t(r))));
 }
 
+template<class T1,class C1>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+  typename CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::Leaf_t>,
+  typename BinaryReturn<C1,typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPType<T1,C1> & l,const typename WordType<C1>::Type_t & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
+    typename CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPType<T1,C1> >::make(l),
+    CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t(r))));
+}
+
 template<class T1,class C1>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<QDPType<T1,C1> >::Leaf_t,
@@ -4033,6 +4115,22 @@ localInnerProductReal(const typename WordType<C2>::Type_t & l,const QDPType<T2,C
     CreateLeaf<QDPType<T2,C2> >::make(r)));
 }
 
+template<class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
+  typename CreateLeaf<QDPType<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const typename WordType<C2>::Type_t & l,const QDPType<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
+    typename CreateLeaf<QDPType<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t(l)),
+    CreateLeaf<QDPType<T2,C2> >::make(r)));
+}
+
 template<class T2,class C2>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
@@ -5547,6 +5645,22 @@ localInnerProductReal(const QDPExpr<T1,C1> & l,const QDPExpr<T2,C2> & r)
     CreateLeaf<QDPExpr<T2,C2> >::make(r)));
 }
 
+template<class T1,class C1,class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+  typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPExpr<T1,C1> & l,const QDPExpr<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+    typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPExpr<T1,C1> >::make(l),
+    CreateLeaf<QDPExpr<T2,C2> >::make(r)));
+}
+
 template<class T1,class C1,class T2,class C2>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
@@ -6129,6 +6243,22 @@ localInnerProductReal(const QDPExpr<T1,C1> & l,const typename WordType<C1>::Type
     CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t(r))));
 }
 
+template<class T1,class C1>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+  typename CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::Leaf_t>,
+  typename BinaryReturn<C1,typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const QDPExpr<T1,C1> & l,const typename WordType<C1>::Type_t & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
+    typename CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<C1,typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<QDPExpr<T1,C1> >::make(l),
+    CreateLeaf<typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C1>::Type_t>::Type_t(r))));
+}
+
 template<class T1,class C1>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<QDPExpr<T1,C1> >::Leaf_t,
@@ -6673,6 +6803,22 @@ localInnerProductReal(const typename WordType<C2>::Type_t & l,const QDPExpr<T2,C
     CreateLeaf<QDPExpr<T2,C2> >::make(r)));
 }
 
+template<class T2,class C2>
+inline typename MakeReturn<BinaryNode<FnLocalColorInnerProduct,
+  typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
+  typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t>,
+  typename BinaryReturn<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t,C2,FnLocalColorInnerProduct>::Type_t >::Expression_t
+localColorInnerProduct(const typename WordType<C2>::Type_t & l,const QDPExpr<T2,C2> & r)
+{
+  typedef BinaryNode<FnLocalColorInnerProduct,
+    typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
+    typename CreateLeaf<QDPExpr<T2,C2> >::Leaf_t> Tree_t;
+  typedef typename BinaryReturn<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t,C2,FnLocalColorInnerProduct>::Type_t Container_t;
+  return MakeReturn<Tree_t,Container_t>::make(Tree_t(
+    CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::make(typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t(l)),
+    CreateLeaf<QDPExpr<T2,C2> >::make(r)));
+}
+
 template<class T2,class C2>
 inline typename MakeReturn<BinaryNode<FnQuarkContract13,
   typename CreateLeaf<typename SimpleScalar<typename WordType<C2>::Type_t>::Type_t >::Leaf_t,
diff --git a/include/qdp_globalfuncs_subtype.h b/include/qdp_globalfuncs_subtype.h
index d1accaa0a..e106ed7b3 100644
--- a/include/qdp_globalfuncs_subtype.h
+++ b/include/qdp_globalfuncs_subtype.h
@@ -99,6 +99,51 @@ namespace QDP
 #endif
 
 
+  template<class T1,class C1,class T2,class C2>
+  typename QDPSubTypeTrait< typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Type_t
+  localColorInnerProduct(const QDPSubType<T1,C1> & l,const QDPType<T2,C2> & r)
+  {
+    if (!l.getOwnsMemory())
+      QDP_error_exit("localColorInnerProduct with subtype view called");
+
+    typename QDPSubTypeTrait< typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Type_t ret;
+    ret.setSubset( l.subset() );
+
+    //QDP_info("localColorInnerProduct %d sites",l.subset().numSiteTable());
+
+    const int *tab = l.subset().siteTable().slice();
+    for(int j=0; j < l.subset().numSiteTable(); ++j)
+      {
+	int i = tab[j];
+	FnLocalColorInnerProduct op;
+	ret.getF()[j] = op( l.getF()[j] , r.elem(i) );
+      }
+
+    return ret;
+  }
+
+  template<class T1,class C1,class T2,class C2>
+  typename QDPSubTypeTrait< typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Type_t
+  localColorInnerProduct(const QDPType<T1,C1> & l,const QDPSubType<T2,C2> & r)
+  {
+    if (!r.getOwnsMemory())
+      QDP_error_exit("localColorInnerProduct with subtype view called");
+
+    typename QDPSubTypeTrait< typename BinaryReturn<C1,C2,FnLocalColorInnerProduct>::Type_t >::Type_t ret;
+    ret.setSubset( r.subset() );
+
+    const int *tab = r.subset().siteTable().slice();
+    for(int j=0; j < r.subset().numSiteTable(); ++j)
+      {
+	int i = tab[j];
+	FnLocalColorInnerProduct op;
+	ret.getF()[j] = op( l.elem(i) , r.getF()[j] );
+      }
+
+    return ret;
+  }
+
+
 
 template<class T>
 typename UnaryReturn<OLattice<T>, FnSum>::Type_t
diff --git a/include/qdp_inner.h b/include/qdp_inner.h
index ac4c5a714..dabd585f7 100644
--- a/include/qdp_inner.h
+++ b/include/qdp_inner.h
@@ -833,11 +833,6 @@ struct UnaryReturn<IScalar<T>, FnPeekSite> {
   typedef IScalar<typename UnaryReturn<T, FnPeekSite>::Type_t>  Type_t;
 };
 
-template<class T>
-struct UnaryReturn<IScalar<T>, FnSumMulti> {
-  typedef IScalar<typename UnaryReturn<T, FnSumMulti>::Type_t>  Type_t;
-};
-
 
 // Gamma algebra
 template<int N, int m, class T2, class OpGammaConstMultiply>
@@ -1759,6 +1754,12 @@ globalMin(const IScalar<T>& s1)
 }
 
 
+// sumMulti
+template<class T>
+struct UnaryReturn<IScalar<T>, FnSumMulti> {
+  typedef IScalar<typename UnaryReturn<T, FnSumMulti>::Type_t>  Type_t;
+};
+
 // InnerProduct (norm-seq) global sum = sum(tr(adj(s1)*s1))
 template<class T>
 struct UnaryReturn<IScalar<T>, FnNorm2> {
@@ -1798,6 +1799,21 @@ localInnerProduct(const IScalar<T1>& s1, const IScalar<T2>& s2)
 }
 
 
+//! IScalar = localColorInnerProduct(adj(IScalar)*IScalar)
+template<class T1, class T2>
+struct BinaryReturn<IScalar<T1>, IScalar<T2>, FnLocalColorInnerProduct> {
+  typedef IScalar<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
+template<class T1, class T2>
+inline typename BinaryReturn<IScalar<T1>, IScalar<T2>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const IScalar<T1>& s1, const IScalar<T2>& s2)
+{
+  return localColorInnerProduct(s1.elem(), s2.elem());
+}
+
+
+
 //! IScalar = where(IScalar, IScalar, IScalar)
 /*!
  * Where is the ? operation
diff --git a/include/qdp_outer.h b/include/qdp_outer.h
index e9ad9fdea..2bae99e88 100644
--- a/include/qdp_outer.h
+++ b/include/qdp_outer.h
@@ -770,6 +770,11 @@ struct BinaryReturn<OScalar<T1>, OScalar<T2>, FnLocalInnerProductReal > {
   typedef OScalar<typename BinaryReturn<T1, T2, FnLocalInnerProductReal>::Type_t>  Type_t;
 };
 
+template<class T1, class T2>
+struct BinaryReturn<OScalar<T1>, OScalar<T2>, FnLocalColorInnerProduct > {
+  typedef OScalar<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
 
 // Gamma algebra
 template<int N, int m, class T2, class OpGammaConstMultiply>
@@ -998,6 +1003,11 @@ struct BinaryReturn<OLattice<T1>, OLattice<T2>, FnLocalInnerProductReal > {
   typedef OLattice<typename BinaryReturn<T1, T2, FnLocalInnerProductReal>::Type_t>  Type_t;
 };
 
+template<class T1, class T2>
+struct BinaryReturn<OLattice<T1>, OLattice<T2>, FnLocalColorInnerProduct > {
+  typedef OLattice<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
 
 // Gamma algebra
 template<int N, int m, class T2, class OpGammaConstMultiply>
@@ -1199,6 +1209,11 @@ struct BinaryReturn<OLattice<T1>, OScalar<T2>, FnLocalInnerProductReal > {
   typedef OLattice<typename BinaryReturn<T1, T2, FnLocalInnerProductReal>::Type_t>  Type_t;
 };
 
+template<class T1, class T2>
+struct BinaryReturn<OLattice<T1>, OScalar<T2>, FnLocalColorInnerProduct > {
+  typedef OLattice<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
 template<class T1, class T2>
 struct BinaryReturn<OScalar<T1>, OLattice<T2>, FnLocalInnerProduct > {
   typedef OLattice<typename BinaryReturn<T1, T2, FnLocalInnerProduct>::Type_t>  Type_t;
@@ -1209,6 +1224,11 @@ struct BinaryReturn<OScalar<T1>, OLattice<T2>, FnLocalInnerProductReal > {
   typedef OLattice<typename BinaryReturn<T1, T2, FnLocalInnerProductReal>::Type_t>  Type_t;
 };
 
+template<class T1, class T2>
+struct BinaryReturn<OScalar<T1>, OLattice<T2>, FnLocalColorInnerProduct > {
+  typedef OLattice<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
 
 // Local operations
 template<class T1, class T2 >
diff --git a/include/qdp_primcolormat.h b/include/qdp_primcolormat.h
index b4855b1e9..cf6e40e87 100644
--- a/include/qdp_primcolormat.h
+++ b/include/qdp_primcolormat.h
@@ -187,6 +187,11 @@ struct UnaryReturn<PColorMatrix<T,N>, FnImagTrace > {
   typedef PScalar<typename UnaryReturn<T, FnImagTrace>::Type_t>  Type_t;
 };
 
+template<class T, int N>
+struct UnaryReturn<PColorMatrix<T,N>, FnSumMulti > {
+  typedef PColorMatrix<typename UnaryReturn<T, FnSumMulti>::Type_t, N>  Type_t;
+};
+
 template<class T, int N>
 struct UnaryReturn<PColorMatrix<T,N>, FnNorm2 > {
   typedef PScalar<typename UnaryReturn<T, FnNorm2>::Type_t>  Type_t;
diff --git a/include/qdp_primcolorvec.h b/include/qdp_primcolorvec.h
index 86807bbaa..67f829375 100644
--- a/include/qdp_primcolorvec.h
+++ b/include/qdp_primcolorvec.h
@@ -153,6 +153,11 @@ struct BinaryReturn<PColorVector<T1,N>, PScalar<T2>, OpDivideAssign > {
  
 
 // ColorVector
+template<class T, int N>
+struct UnaryReturn<PColorVector<T,N>, FnSumMulti > {
+  typedef PColorVector<typename UnaryReturn<T, FnSumMulti>::Type_t, N>  Type_t;
+};
+
 template<class T, int N>
 struct UnaryReturn<PColorVector<T,N>, FnNorm2 > {
   typedef PScalar<typename UnaryReturn<T, FnNorm2>::Type_t>  Type_t;
@@ -190,6 +195,28 @@ struct BinaryReturn<PColorVector<T1,N>, PColorVector<T2,N>, FnLocalInnerProductR
 // Operators
 //-----------------------------------------------------------------------------
 
+//! PScalar = localColorInnerProduct(PColorVector,PColorVector)
+template<class T1, class T2, int N>
+struct BinaryReturn<PColorVector<T1,N>, PColorVector<T2,N>, FnLocalColorInnerProduct> {
+  typedef PScalar<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
+template<class T1, class T2, int N>
+inline typename BinaryReturn<PColorVector<T1,N>, PColorVector<T2,N>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const PColorVector<T1,N>& l, const PColorVector<T2,N>& r)
+{
+  typename BinaryReturn<PColorVector<T1,N>, PColorVector<T2,N>, FnLocalColorInnerProduct>::Type_t  d;
+
+  // The traceColor is eaten here
+  d.elem() = adj(l.elem(0)) * r.elem(0);
+  for(int k=1; k < N; ++k)
+    d.elem() += adj(l.elem(k)) * r.elem(k);
+
+  return d;
+}
+
+
+//-----------------------------------------------
 // Peeking and poking
 //! Extract color vector components 
 template<class T, int N>
diff --git a/include/qdp_primscalar.h b/include/qdp_primscalar.h
index 9448e3c10..af68e119e 100644
--- a/include/qdp_primscalar.h
+++ b/include/qdp_primscalar.h
@@ -1415,6 +1415,12 @@ sum(const PScalar<T>& s1)
 #endif
 
 
+// ! PScalar
+template<class T>
+struct UnaryReturn<PScalar<T>, FnSumMulti > {
+  typedef PScalar<typename UnaryReturn<T, FnSumMulti>::Type_t>  Type_t;
+};
+
 // InnerProduct (norm-seq) global sum = sum(tr(adj(s1)*s1))
 template<class T>
 struct UnaryReturn<PScalar<T>, FnNorm2 > {
@@ -1499,6 +1505,20 @@ localInnerProductReal(const PScalar<T1>& s1, const PScalar<T2>& s2)
 }
 
 
+//! PScalar<T> = localColorInnerProduct(adj(PScalar<T1>)*PScalar<T2>)
+template<class T1, class T2>
+struct BinaryReturn<PScalar<T1>, PScalar<T2>, FnLocalColorInnerProduct > {
+  typedef PScalar<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
+template<class T1, class T2>
+inline typename BinaryReturn<PScalar<T1>, PScalar<T2>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const PScalar<T1>& s1, const PScalar<T2>& s2)
+{
+  return localColorInnerProduct(s1.elem(), s2.elem());
+}
+
+
 //! PScalar<T> = where(PScalar, PScalar, PScalar)
 /*!
  * Where is the ? operation
diff --git a/include/qdp_primspinmat.h b/include/qdp_primspinmat.h
index c2408b6df..1ad8cfe87 100644
--- a/include/qdp_primspinmat.h
+++ b/include/qdp_primspinmat.h
@@ -212,6 +212,11 @@ struct UnaryReturn<PSpinMatrix<T,N>, FnImagTrace > {
   typedef PScalar<typename UnaryReturn<T, FnImagTrace>::Type_t>  Type_t;
 };
 
+template<class T, int N>
+struct UnaryReturn<PSpinMatrix<T,N>, FnSumMulti > {
+  typedef PSpinMatrix<typename UnaryReturn<T, FnSumMulti>::Type_t, N>  Type_t;
+};
+
 template<class T, int N>
 struct UnaryReturn<PSpinMatrix<T,N>, FnNorm2 > {
   typedef PScalar<typename UnaryReturn<T, FnNorm2>::Type_t>  Type_t;
@@ -465,6 +470,31 @@ transposeSpin(const PSpinMatrix<T,N>& s1)
 }
 
 
+//-----------------------------------------------
+//! PScalar = localColorInnerProduct(PSpinMatrix,PSpinMatrix)
+template<class T1, class T2, int N>
+struct BinaryReturn<PSpinMatrix<T1,N>, PSpinMatrix<T2,N>, FnLocalColorInnerProduct> {
+  typedef PSpinMatrix<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t, N>  Type_t;
+};
+
+template<class T1, class T2, int N>
+inline typename BinaryReturn<PSpinMatrix<T1,N>, PSpinMatrix<T2,N>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const PSpinMatrix<T1,N>& l, const PSpinMatrix<T2,N>& r)
+{
+  typename BinaryReturn<PSpinMatrix<T1,N>, PSpinMatrix<T2,N>, FnLocalColorInnerProduct>::Type_t  d;
+
+  for(int i=0; i < N; ++i)
+    for(int j=0; j < N; ++j)
+    {
+      d.elem(i,j) = localColorInnerProduct(l.elem(0,i), r.elem(0,j));
+      for(int k=1; k < N; ++k)
+	d.elem(i,j) += localColorInnerProduct(l.elem(k,i), r.elem(k,j));
+    }
+
+  return d;
+}
+
+
 //-----------------------------------------------
 // OuterProduct must be handled specially for each color and spin
 // The problem is the traits class - I have no way to say to PVector's
diff --git a/include/qdp_primspinvec.h b/include/qdp_primspinvec.h
index b6e3fbae4..fbf458c30 100644
--- a/include/qdp_primspinvec.h
+++ b/include/qdp_primspinvec.h
@@ -198,7 +198,6 @@ XMLWriter& operator<<(XMLWriter& xml, const PSpinVector<T,N>& d)
 
 
 // Primitive Vectors
-
 template<class T1, int N>
 inline typename UnaryReturn<PSpinVector<T1,N>, OpUnaryPlus>::Type_t
 operator+(const PSpinVector<T1,N>& l)
@@ -903,6 +902,11 @@ struct BinaryReturn<PSpinVector<T1,N>, PScalar<T2>, OpDivideAssign > {
 
 
 // SpinVector
+template<class T, int N>
+struct UnaryReturn<PSpinVector<T,N>, FnSumMulti > {
+  typedef PSpinVector<typename UnaryReturn<T, FnSumMulti>::Type_t, N>  Type_t;
+};
+
 template<class T, int N>
 struct UnaryReturn<PSpinVector<T,N>, FnNorm2 > {
   typedef PScalar<typename UnaryReturn<T, FnNorm2>::Type_t>  Type_t;
diff --git a/include/qdp_reality.h b/include/qdp_reality.h
index 29989d942..18da0a12d 100644
--- a/include/qdp_reality.h
+++ b/include/qdp_reality.h
@@ -1614,6 +1614,11 @@ globalMin(const RScalar<T>& s1)
 
 
 //------------------------------------------
+template<class T>
+struct UnaryReturn<RScalar<T>, FnSumMulti > {
+  typedef RScalar<typename UnaryReturn<T, FnSumMulti>::Type_t>  Type_t;
+};
+
 // InnerProduct (norm-seq) global sum = sum(tr(adj(s1)*s1))
 template<class T>
 struct UnaryReturn<RScalar<T>, FnNorm2 > {
@@ -1673,6 +1678,20 @@ localInnerProductReal(const RScalar<T1>& s1, const RScalar<T2>& s2)
 }
 
 
+//! RScalar<T> = localInnerProduct(adj(RScalar<T1>)*RScalar<T2>)
+template<class T1, class T2>
+struct BinaryReturn<RScalar<T1>, RScalar<T2>, FnLocalColorInnerProduct > {
+  typedef RScalar<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
+template<class T1, class T2>
+inline typename BinaryReturn<RScalar<T1>, RScalar<T2>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const RScalar<T1>& s1, const RScalar<T2>& s2)
+{
+  return localColorInnerProduct(s1.elem(), s2.elem());
+}
+
+
 //! RScalar<T> = where(RScalar, RScalar, RScalar)
 /*!
  * Where is the ? operation
@@ -2344,6 +2363,12 @@ sum(const RComplex<T>& s1)
 #endif
 
 
+// Sum
+template<class T>
+struct UnaryReturn<RComplex<T>, FnSumMulti > {
+  typedef RComplex<typename UnaryReturn<T, FnSumMulti>::Type_t>  Type_t;
+};
+
 // InnerProduct (norm-seq) global sum = sum(tr(adj(s1)*s1))
 template<class T>
 struct UnaryReturn<RComplex<T>, FnNorm2 > {
@@ -2406,6 +2431,24 @@ localInnerProductReal(const RComplex<T1>& l, const RComplex<T2>& r)
 }
 
 
+//! RComplex<T> = localColorInnerProduct(adj(RComplex<T1>)*RComplex<T2>)
+template<class T1, class T2>
+struct BinaryReturn<RComplex<T1>, RComplex<T2>, FnLocalColorInnerProduct > {
+  typedef RComplex<typename BinaryReturn<T1, T2, FnLocalColorInnerProduct>::Type_t>  Type_t;
+};
+
+template<class T1, class T2>
+inline typename BinaryReturn<RComplex<T1>, RComplex<T2>, FnLocalColorInnerProduct>::Type_t
+localColorInnerProduct(const RComplex<T1>& l, const RComplex<T2>& r)
+{
+  typedef typename BinaryReturn<RComplex<T1>, RComplex<T2>, FnLocalColorInnerProduct>::Type_t  Ret_t;
+
+  return Ret_t(localColorInnerProduct(l.real(),r.real()) + localColorInnerProduct(l.imag(),r.imag()),
+	       localColorInnerProduct(l.real(),r.imag()) - localColorInnerProduct(l.imag(),r.real()));
+}
+
+
+
 //! RComplex<T> = where(RScalar, RComplex, RComplex)
 /*!
  * Where is the ? operation
diff --git a/include/qdp_scalar_specific.h b/include/qdp_scalar_specific.h
index 372a90308..fe4002699 100644
--- a/include/qdp_scalar_specific.h
+++ b/include/qdp_scalar_specific.h
@@ -678,7 +678,7 @@ sum(const QDPExpr<RHS,OLattice<T> >& s1)
  * slices
  */
 template<class RHS, class T>
-typename UnaryReturn<OScalar<T>, FnSum>::Type_t
+typename UnaryReturn<OScalar<T>, FnSumMulti>::Type_t
 sumMulti(const QDPExpr<RHS,OScalar<T> >& s1, const Set& ss)
 {
   typename UnaryReturn<OScalar<T>, FnSumMulti>::Type_t  dest(ss.numSubsets());
@@ -838,7 +838,7 @@ sumMulti(const QDPExpr<RHS,OLattice<T> >& s1, const Set& ss)
  * slices
  */
 template<class T>
-multi2d<typename UnaryReturn<OScalar<T>, FnSum>::Type_t>
+multi2d<typename UnaryReturn<OScalar<T>, FnSumMulti>::Type_t>
 sumMulti(const multi1d< OScalar<T> >& s1, const Set& ss)
 {
   multi2d<typename UnaryReturn<OScalar<T>, FnSum>::Type_t>  dest(s1.size(),ss.numSubsets());
@@ -873,7 +873,7 @@ sumMulti(const multi1d< OScalar<T> >& s1, const Set& ss)
  * version is fine.
  */
 template<class T>
-multi2d<typename UnaryReturn<OLattice<T>, FnSum>::Type_t>
+multi2d<typename UnaryReturn<OLattice<T>, FnSumMulti>::Type_t>
 sumMulti(const multi1d< OLattice<T> >& s1, const Set& ss)
 {
   multi2d<typename UnaryReturn<OLattice<T>, FnSum>::Type_t>  dest(s1.size(),ss.numSubsets());
diff --git a/include/qdp_scalarsite_defs.h b/include/qdp_scalarsite_defs.h
index 4c368d05e..497874c68 100644
--- a/include/qdp_scalarsite_defs.h
+++ b/include/qdp_scalarsite_defs.h
@@ -36,6 +36,7 @@ typedef OLattice< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 4> > LatticeDi
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 1> > LatticeStaggeredFermion;
 typedef OLattice< PSpinMatrix< PColorMatrix< RComplex<REAL>, Nc>, 4> > LatticeDiracPropagator;
 typedef OLattice< PSpinMatrix< PColorMatrix< RComplex<REAL>, Nc>, 1> > LatticeStaggeredPropagator;
+typedef OLattice< PSpinMatrix< PColorVector< RComplex<REAL>, Nc>, 4> > LatticeColorVectorSpinMatrix;
 
 typedef OScalar< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 4> > DiracFermion;
 typedef OScalar< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 1> > StaggeredFermion;
@@ -196,6 +197,8 @@ typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, 3>, 4> > LatticeD
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, 2>, 4> > LatticeDiracFermionF2;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, 1>, 4> > LatticeDiracFermionF1;
 
+typedef OLattice< PSpinMatrix< PColorVector< RComplex<REAL32>, Nc>, 4> > LatticeColorVectorSpinMatrixF;
+
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, Nc>, 1> > LatticeStaggeredFermionF;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, 3>, 1> > LatticeStaggeredFermionF3;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL32>, 2>, 1> > LatticeStaggeredFermionF2;
@@ -315,6 +318,8 @@ typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, 3>, 4> > LatticeD
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, 2>, 4> > LatticeDiracFermionD2;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, 1>, 4> > LatticeDiracFermionD1;
 
+typedef OLattice< PSpinMatrix< PColorVector< RComplex<REAL64>, Nc>, 4> > LatticeColorVectorSpinMatrixD;
+
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, Nc>, 1> > LatticeStaggeredFermionD;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, 3>, 1> > LatticeStaggeredFermionD3;
 typedef OLattice< PSpinVector< PColorVector< RComplex<REAL64>, 2>, 1> > LatticeStaggeredFermionD2;
@@ -445,6 +450,7 @@ typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 4> > SubLat
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL>, Nc>, 1> > SubLatticeStaggeredFermion;
 typedef OSubLattice< PSpinMatrix< PColorMatrix< RComplex<REAL>, Nc>, 4> > SubLatticeDiracPropagator;
 typedef OSubLattice< PSpinMatrix< PColorMatrix< RComplex<REAL>, Nc>, 1> > SubLatticeStaggeredPropagator;
+typedef OSubLattice< PSpinMatrix< PColorVector< RComplex<REAL>, Nc>, 4> > SubLatticeColorVectorSpinMatrix;
 
 // Floating aliases
 typedef OSubLattice< PScalar< PColorVector< RComplex<REAL>, Nc> > > SubLatticeColorVector;
@@ -469,6 +475,7 @@ typedef OSubLattice< PSpinVector< PColorVector< RComplex<DOUBLE>, Nc>, 4> > SubL
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<DOUBLE>, Nc>, 1> > SubLatticeDStaggeredFermion;
 typedef OSubLattice< PSpinMatrix< PColorMatrix< RComplex<DOUBLE>, Nc>, 4> > SubLatticeDDiracPropagator;
 typedef OSubLattice< PSpinMatrix< PColorMatrix< RComplex<DOUBLE>, Nc>, 1> > SubLatticeDStaggeredPropagator;
+typedef OSubLattice< PSpinMatrix< PColorVector< RComplex<DOUBLE>, Nc>, 4> > SubLatticeDColorVectorSpinMatrix;
 
 // Floating aliases
 typedef OSubLattice< PScalar< PColorVector< RComplex<DOUBLE>, Nc> > > SubLatticeDColorVector;
@@ -526,6 +533,8 @@ typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, 3>, 4> > SubLa
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, 2>, 4> > SubLatticeDiracFermionF2;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, 1>, 4> > SubLatticeDiracFermionF1;
 
+typedef OSubLattice< PSpinMatrix< PColorVector< RComplex<REAL32>, Nc>, 4> > SubLatticeColorVectorSpinMatrixF;
+
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, Nc>, 1> > SubLatticeStaggeredFermionF;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, 3>, 1> > SubLatticeStaggeredFermionF3;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL32>, 2>, 1> > SubLatticeStaggeredFermionF2;
@@ -578,6 +587,8 @@ typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, 3>, 4> > SubLa
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, 2>, 4> > SubLatticeDiracFermionD2;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, 1>, 4> > SubLatticeDiracFermionD1;
 
+typedef OSubLattice< PSpinMatrix< PColorVector< RComplex<REAL64>, Nc>, 4> > SubLatticeColorVectorSpinMatrixD;
+
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, Nc>, 1> > SubLatticeStaggeredFermionD;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, 3>, 1> > SubLatticeStaggeredFermionD3;
 typedef OSubLattice< PSpinVector< PColorVector< RComplex<REAL64>, 2>, 1> > SubLatticeStaggeredFermionD2;
diff --git a/include/qdp_simpleword.h b/include/qdp_simpleword.h
index a74e09728..1e1158ba6 100644
--- a/include/qdp_simpleword.h
+++ b/include/qdp_simpleword.h
@@ -376,6 +376,12 @@ int localInnerProduct(int s1, int s2)
   return s1*s2;
 }
 
+inline
+int localColorInnerProduct(int s1, int s2)
+{
+  return s1*s2;
+}
+
 inline
 unsigned int sum(unsigned int s1)
 {
@@ -393,6 +399,11 @@ unsigned int localInnerProduct(unsigned int s1, unsigned int s2)
 {
   return s1*s2;
 }
+inline
+unsigned int localColorInnerProduct(unsigned int s1, unsigned int s2)
+{
+  return s1*s2;
+}
 
 inline
 double sum(float s1)
@@ -411,6 +422,11 @@ double localInnerProduct(float s1, float s2)
 {
   return double(s1*s2);
 }
+inline
+double localColorInnerProduct(float s1, float s2)
+{
+  return double(s1*s2);
+}
 
 inline
 double sum(double s1)
@@ -430,18 +446,36 @@ double localInnerProduct(float s1, double s2)
   return double(s1)*s2;
 }
 
+inline
+double localColorInnerProduct(float s1, double s2)
+{
+  return double(s1)*s2;
+}
+
 inline
 double localInnerProduct(double s1, float s2)
 {
   return s1*double(s2);
 }
 
+inline
+double localColorInnerProduct(double s1, float s2)
+{
+  return s1*double(s2);
+}
+
 inline
 double localInnerProduct(double s1, double s2)
 {
   return s1*s2;
 }
 
+inline
+double localColorInnerProduct(double s1, double s2)
+{
+  return s1*s2;
+}
+
 
 
 /*! @} */  // end of group simpleword
@@ -496,6 +530,11 @@ struct BinaryReturn<int, int, FnLocalInnerProduct> {
   typedef int  Type_t;
 };
 
+template<>
+struct BinaryReturn<int, int, FnLocalColorInnerProduct> {
+  typedef int  Type_t;
+};
+
 template<>
 struct TrinaryReturn<bool, int, int, FnWhere> {
   typedef int  Type_t;
@@ -542,6 +581,11 @@ struct BinaryReturn<float, float, FnLocalInnerProduct> {
   typedef double  Type_t;
 };
 
+template<>
+struct BinaryReturn<float, float, FnLocalColorInnerProduct> {
+  typedef double  Type_t;
+};
+
 template<>
 struct TrinaryReturn<bool, float, float, FnWhere> {
   typedef float  Type_t;
@@ -588,6 +632,11 @@ struct BinaryReturn<double, double, FnLocalInnerProduct> {
   typedef double  Type_t;
 };
 
+template<>
+struct BinaryReturn<double, double, FnLocalColorInnerProduct> {
+  typedef double  Type_t;
+};
+
 template<>
 struct TrinaryReturn<bool, double, double, FnWhere> {
   typedef double  Type_t;
diff --git a/other_libs/filedb b/other_libs/filedb
index 7ddc3e5d4..6937e71d6 160000
--- a/other_libs/filedb
+++ b/other_libs/filedb
@@ -1 +1 @@
-Subproject commit 7ddc3e5d41d5e0eaad72d1064b331bdb37bc8bf1
+Subproject commit 6937e71d693fdc7da8b696e2077a2a6413a09a01