diff --git a/lib/meas/glue/mesplq.cc b/lib/meas/glue/mesplq.cc
index a88f5ef83..1f6c06391 100644
--- a/lib/meas/glue/mesplq.cc
+++ b/lib/meas/glue/mesplq.cc
@@ -24,7 +24,7 @@ namespace Chroma
    
 
   template<typename Q>
-  void MesPlq_t(const multi1d<Q>& u, 
+  void MesPlq(const multi1d<Q>& u, 
 	      multi2d<Double>& plane_plaq, Double& link)
   {
     START_CODE();
@@ -41,10 +41,10 @@ namespace Chroma
 	// This is the the longer way to write the 1-liner in the else clause
 
 	/* tmp_0 = u(x+mu,nu)*u_dag(x+nu,mu) */
-	LatticeColorMatrix tmp_0 = shift(u[nu],FORWARD,mu) * adj(shift(u[mu],FORWARD,nu));
+	Q tmp_0 = shift(u[nu],FORWARD,mu) * adj(shift(u[mu],FORWARD,nu));
 
 	/* tmp_1 = tmp_0*u_dag(x,nu)=u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu) */
-	LatticeColorMatrix tmp_1 = tmp_0 * adj(u[nu]);
+	Q tmp_1 = tmp_0 * adj(u[nu]);
 
 	/* tmp = sum(tr(u(x,mu)*tmp_1=u(x,mu)*u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu))) */
 	Double tmp = sum(real(trace(u[mu]*tmp_1)));
@@ -80,17 +80,6 @@ namespace Chroma
     END_CODE();
   }
 
-  void MesPlq(const multi1d<LatticeColorMatrixF3>& u, 
-	      multi2d<Double>& plane_plaq, Double& link) 
-  {
-      MesPlq_t(u,plane_plaq, link);
-  }
-
-  void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
-	      multi2d<Double>& plane_plaq, Double& link)
-  {
-      MesPlq_t(u,plane_plaq, link);
-  }
 
   //! Return the value of the average plaquette normalized to 1
   /*!
@@ -104,7 +93,7 @@ namespace Chroma
    * \param link        space-time average link (Write)
    */
   template<typename Q>
-  void MesPlq_t(const multi1d<Q>& u, 
+  void MesPlq(const multi1d<Q>& u, 
 	      Double& w_plaq, Double& s_plaq, Double& t_plaq, 
 	      multi2d<Double>& plane_plaq,
 	      Double& link)
@@ -143,21 +132,6 @@ namespace Chroma
     END_CODE();
   }
 
-  void MesPlq(const multi1d<LatticeColorMatrixF3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, 
-	      multi2d<Double>& plane_plaq,
-	      Double& link) 
-  {
-     MesPlq_t(u,w_plaq,s_plaq,t_plaq, plane_plaq, link);
-  }
-
-  void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, 
-	      multi2d<Double>& plane_plaq,
-	      Double& link)
-  {
-     MesPlq_t(u,w_plaq,s_plaq,t_plaq, plane_plaq, link);
-  }
 
   //! Return the value of the average plaquette normalized to 1
   /*!
@@ -170,29 +144,6 @@ namespace Chroma
    * \param link        space-time average link (Write)
    */
 
-  void MesPlq(const multi1d<LatticeColorMatrixF3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link)
-  {
-    START_CODE();
-
-    multi2d<Double> plane_plaq;
-
-    MesPlq(u, w_plaq, s_plaq, t_plaq, plane_plaq, link);
-
-    END_CODE();
-  }
- 
-  void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link)
-  {
-    START_CODE();
-
-    multi2d<Double> plane_plaq;
-
-    MesPlq(u, w_plaq, s_plaq, t_plaq, plane_plaq, link);
-
-    END_CODE();
-  }
 
   //! Print the value of the average plaquette normalized to 1
   /*!
@@ -203,7 +154,7 @@ namespace Chroma
    * \param u          gauge field (Read)
    */
   template<typename Q>
-  void MesPlq_t(XMLWriter& xml, 
+  void MesPlq(XMLWriter& xml, 
 	        const std::string& xml_group,
 	        const multi1d<Q>& u)
   {
@@ -264,18 +215,49 @@ namespace Chroma
     END_CODE();
   }
 
- void MesPlq(XMLWriter& xml, 
-	     const std::string& xml_group,
-	     const multi1d<LatticeColorMatrixF3>& u)
- {
-   MesPlq_t(xml, xml_group, u);
- }
+  // ***** Instantiate templates for needed types *******
 
- void MesPlq(XMLWriter& xml, 
+  template void MesPlq(const multi1d<LatticeColorMatrix>& u, 
+		       multi2d<Double>& plane_plaq, Double& link);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
+		       multi2d<Double>& plane_plaq, Double& link);
+#endif
+
+  template void MesPlq(const multi1d<LatticeColorMatrix>& u,
+		       Double& w_plaq, Double& s_plaq, Double& t_plaq, 
+		       multi2d<Double>& plane_plaq,
+		       Double& link);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
+		       Double& w_plaq, Double& s_plaq, Double& t_plaq, 
+		       multi2d<Double>& plane_plaq,
+		       Double& link);
+#endif
+
+  template<typename Q>
+  void MesPlq(const multi1d<Q>& u, 
+	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link)
+  {
+    START_CODE();
+    multi2d<Double> plane_plaq;
+    MesPlq(u, w_plaq, s_plaq, t_plaq, plane_plaq, link);
+    END_CODE();
+  }
+  template void MesPlq(const multi1d<LatticeColorMatrix>& u, 
+	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
+	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link);
+#endif
+
+  template void MesPlq(XMLWriter& xml, 
+		       const std::string& xml_group,
+		       const multi1d<LatticeColorMatrix>& u);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void MesPlq(XMLWriter& xml, 
 	        const std::string& xml_group,
-	        const multi1d<LatticeColorMatrixD3>& u)
- {
-   MesPlq_t(xml, xml_group, u);
- }
+		       const multi1d<LatticeColorMatrixD3>& u);
+#endif
 
 }  // end namespace Chroma
diff --git a/lib/meas/glue/mesplq.h b/lib/meas/glue/mesplq.h
index d63e239bf..0b6df6db7 100644
--- a/lib/meas/glue/mesplq.h
+++ b/lib/meas/glue/mesplq.h
@@ -6,21 +6,9 @@
 namespace Chroma 
 {
 
-  //! Return the value of the average plaquette normalized to 1
-  /*!
-   * \ingroup glue
-   *
-   * \param u         gauge field (Read)
-   * \param w_plaq    plaquette average (Write)
-   * \param s_plaq    space-like plaquette average (Write)
-   * \param t_plaq    time-like plaquette average (Write)
-   * \param link      space-time average link (Write)
-   */
-  void MesPlq(const multi1d<LatticeColorMatrixF3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link);
-
-  void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link);
+  template<typename Q>
+  void MesPlq(const multi1d<Q>& u, 
+	      multi2d<Double>& plane_plaq, Double& link);
 
   //! Return the value of the average plaquette normalized to 1
   /*!
@@ -34,15 +22,14 @@ namespace Chroma
    * \param link        space-time average link (Write)
    */
 
-  void MesPlq(const multi1d<LatticeColorMatrixF3>& u, 
-	      Double& w_plaq, Double& s_plaq, Double& t_plaq, 
-	      multi2d<Double>& plane_plaq,
-	      Double& link);
-
-  void MesPlq(const multi1d<LatticeColorMatrixD3>& u, 
+  template<typename Q>
+  void MesPlq(const multi1d<Q>& u, 
 	      Double& w_plaq, Double& s_plaq, Double& t_plaq, 
 	      multi2d<Double>& plane_plaq,
 	      Double& link);
+  template<typename Q>
+  void MesPlq(const multi1d<Q>& u, 
+	      Double& w_plaq, Double& s_plaq, Double& t_plaq, Double& link);
 
   //! Print the value of the average plaquette normalized to 1
   /*!
@@ -51,13 +38,10 @@ namespace Chroma
    * \param xml    plaquette average (Write)
    * \param u      gauge field (Read)
    */
-  void MesPlq(XMLWriter& xml,
-	      const std::string& xml_group,
-	      const multi1d<LatticeColorMatrixF3>& u);
-
-  void MesPlq(XMLWriter& xml,
+  template<typename Q>
+  void MesPlq(XMLWriter& xml, 
 	      const std::string& xml_group,
-	      const multi1d<LatticeColorMatrixD3>& u);
+	      const multi1d<Q>& u);
 
 }  // end namespace Chroma
 
diff --git a/lib/meas/glue/polylp.cc b/lib/meas/glue/polylp.cc
index 2539ee7b6..cd2da37fe 100644
--- a/lib/meas/glue/polylp.cc
+++ b/lib/meas/glue/polylp.cc
@@ -17,16 +17,16 @@ namespace Chroma
    * \param mu         direction of Polyakov loop (Read)
    */
   template<typename Q> 
-  void polylp_t(const multi1d<Q>& u, DComplex& poly_loop, int mu)
+  void polylp(const multi1d<Q>& u, DComplex& poly_loop, int mu)
   {
     START_CODE();
         
     // Initial Polyakov loop
-    LatticeColorMatrix poly = u[mu];
+    Q poly = u[mu];
 
     for(int n = 1; n < Layout::lattSize()[mu]; ++n)  // run over all links in mu dir
     {    
-      LatticeColorMatrix tmp = shift(poly, FORWARD, mu);
+      Q tmp = shift(poly, FORWARD, mu);
       poly = u[mu] * tmp;
     }
 
@@ -36,17 +36,6 @@ namespace Chroma
     END_CODE();
   }
 
-
-  void polylp(const multi1d<LatticeColorMatrixF3>& u, DComplex& poly_loop, int mu) 
-  {
-     polylp_t( u, poly_loop, mu);
-  }
-
-  void polylp(const multi1d<LatticeColorMatrixD3>& u, DComplex& poly_loop, int mu)
-  {
-     polylp_t(u, poly_loop, mu);
-  }
-
   //! Compute Polyakov loop
   /*!
    * \ingroup glue
@@ -55,7 +44,7 @@ namespace Chroma
    * \param poly_loop  Polyakov loop average (Write) 
    */
   template<typename Q>
-  void polylp_t(const multi1d<Q>& u, multi1d<DComplex>& poly_loop)
+  void polylp(const multi1d<Q>& u, multi1d<DComplex>& poly_loop)
   {
     START_CODE();
  
@@ -68,13 +57,16 @@ namespace Chroma
   }
 
 
-  void polylp(const multi1d<LatticeColorMatrixF3>& u, multi1d<DComplex>& poly_loop)
-  { 
-      polylp_t(u,poly_loop);
-  }
+  // ***** Instantiate templates for needed types *******
+
+  template void polylp(const multi1d<LatticeColorMatrix>& u, DComplex& poly_loop, int mu);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void polylp(const multi1d<LatticeColorMatrixD3>& u, DComplex& poly_loop, int mu);
+#endif
+
+  template void polylp(const multi1d<LatticeColorMatrix>& u, multi1d<DComplex>& poly_loop);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void polylp(const multi1d<LatticeColorMatrixD3>& u, multi1d<DComplex>& poly_loop);
+#endif
 
-  void polylp(const multi1d<LatticeColorMatrixD3>& u, multi1d<DComplex>& poly_loop)
-  {
-      polylp_t(u,poly_loop);
-  }
 }  // end namespace Chroma
diff --git a/lib/meas/glue/polylp.h b/lib/meas/glue/polylp.h
index 8c586a5a8..cce62f9ae 100644
--- a/lib/meas/glue/polylp.h
+++ b/lib/meas/glue/polylp.h
@@ -18,9 +18,8 @@ namespace Chroma
    * \param mu         direction of Polyakov loop (Read)
    */
 
-  void polylp(const multi1d<LatticeColorMatrixF3>& u, DComplex& poly_loop, int mu);
-
-  void polylp(const multi1d<LatticeColorMatrixD3>& u, DComplex& poly_loop, int mu);
+  template<typename Q> 
+  void polylp(const multi1d<Q>& u, DComplex& poly_loop, int mu);
 
   //! Compute Polyakov loop
   /*!
@@ -30,8 +29,8 @@ namespace Chroma
    * \param poly_loop  Polyakov loop average (Write) 
    */
 
-  void polylp(const multi1d<LatticeColorMatrixF3>& u, multi1d<DComplex>& poly_loop);
-  void polylp(const multi1d<LatticeColorMatrixD3>& u, multi1d<DComplex>& poly_loop);
+  template<typename Q>
+  void polylp(const multi1d<Q>& u, multi1d<DComplex>& poly_loop);
 
 }  // end namespace Chroma
 
diff --git a/lib/util/gauge/reunit.cc b/lib/util/gauge/reunit.cc
index 9bedf749e..c241186a9 100644
--- a/lib/util/gauge/reunit.cc
+++ b/lib/util/gauge/reunit.cc
@@ -24,7 +24,6 @@
 
 namespace Chroma {
 
-
   namespace ReunitEnv {
     static double time_spent =0;
     double getTime() { return time_spent; }
@@ -447,135 +446,116 @@ namespace Chroma {
     ReunitEnv::time_spent += swatch.getTimeInSeconds();
     END_CODE();
   }
-  
+
+
+  // ***** Variable argument lists *******
+
   // Overloaded definitions
-  // SINGLE
-  void reunit(LatticeColorMatrixF3& xa)
+  // Generic
+  void reunit(LatticeColorMatrix& xa)
   {
     START_CODE();
-
     LatticeBoolean bad;
     int numbad;
-    
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, REUNITARIZE, all);
-    
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, REUNITARIZE, all);
     END_CODE();
   }
-
-   // Overloaded definitions
+  // Overloaded definitions
   // DOUBLE
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
   void reunit(LatticeColorMatrixD3& xa)
   {
     START_CODE();
-
     LatticeBoolean bad;
     int numbad;
-    
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, REUNITARIZE, all);
-    
     END_CODE();
   }
+#endif
 
-
-  // SINGLE
-  void reunit(LatticeColorMatrixF3& xa,
+  // Generic
+  void reunit(LatticeColorMatrix& xa,
 	      const Subset& mstag)
   {
     START_CODE();
-
     LatticeBoolean bad;
     int numbad;
-    
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
-    
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
     END_CODE();
   }
-
   // DOUBLE
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
   void reunit(LatticeColorMatrixD3& xa,
 	      const Subset& mstag)
   {
     START_CODE();
-
     LatticeBoolean bad;
     int numbad;
-    
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
-    
     END_CODE();
   }
+#endif
 
   // Overloaded definitions, with numbad and ruflag
-  // Single
-
-  void reunit(LatticeColorMatrixF3& xa,
+  // Generic
+  void reunit(LatticeColorMatrix& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag)
   {
     START_CODE();
-
     LatticeBoolean bad;
-    
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, REUNITARIZE, all);
-    
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, REUNITARIZE, all);
     END_CODE();
   }
-
   // DOUBLE
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
   void reunit(LatticeColorMatrixD3& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag)
   {
     START_CODE();
-
     LatticeBoolean bad;
-    
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, REUNITARIZE, all);
-    
     END_CODE();
   }
-  
+#endif
 
-  // SINGLE
-  void reunit(LatticeColorMatrixF3& xa,
+  // Generic
+  void reunit(LatticeColorMatrix& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag)
   {
     START_CODE();
-
     LatticeBoolean bad;
-    
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
-    
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
     END_CODE();
   }
-  
-   // DOUBLE
-   void reunit(LatticeColorMatrixD3& xa,
+  // DOUBLE
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  void reunit(LatticeColorMatrixD3& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag)
   {
     START_CODE();
-
     LatticeBoolean bad;
-    
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, REUNITARIZE, mstag);
-    
     END_CODE();
   }
+#endif
+
   // Overloaded definitions, with bad, numbad and ruflag
-  // SINGLE
-  void reunit(LatticeColorMatrixF3& xa, 
+  // Generic
+  void reunit(LatticeColorMatrix& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
 	      enum Reunitarize ruflag)
   {
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, ruflag, all);
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, ruflag, all);
   }
-
   // DOUBLE
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
   void reunit(LatticeColorMatrixD3& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
@@ -583,19 +563,19 @@ namespace Chroma {
   {
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, ruflag, all);
   }
+#endif
 
-  // Single
-  void reunit(LatticeColorMatrixF3& xa, 
+  // Generic
+  void reunit(LatticeColorMatrix& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag)
   {
-    reunit_t<LatticeColorMatrixF3, LatticeComplexF, LatticeRealF, Subset>(xa, bad, numbad, ruflag, mstag);
+    reunit_t<LatticeColorMatrix, LatticeComplex, LatticeReal, Subset>(xa, bad, numbad, ruflag, mstag);
   }
-  
-  
  // Double
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
   void reunit(LatticeColorMatrixD3& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
@@ -604,5 +584,6 @@ namespace Chroma {
   {
     reunit_t<LatticeColorMatrixD3, LatticeComplexD, LatticeRealD, Subset>(xa, bad, numbad, ruflag, mstag);
   }
+#endif
 
 } // End namespace
diff --git a/lib/util/gauge/reunit.h b/lib/util/gauge/reunit.h
index c70d5bf38..39ea31c77 100644
--- a/lib/util/gauge/reunit.h
+++ b/lib/util/gauge/reunit.h
@@ -27,53 +27,48 @@ namespace Chroma {
   };
 
   enum Reunitarize {REUNITARIZE, REUNITARIZE_ERROR, REUNITARIZE_LABEL};
-  
-  void reunit(LatticeColorMatrixF3& xa);
+
+
+  void reunit(LatticeColorMatrix& xa);
   void reunit(LatticeColorMatrixD3& xa);
   
-  void reunit(LatticeColorMatrixF3& xa,
+  void reunit(LatticeColorMatrix& xa,
 	      const Subset& mstag);
-
   void reunit(LatticeColorMatrixD3& xa,
 	      const Subset& mstag);
   
   // With ruflag
-  void reunit(LatticeColorMatrixF3& xa,
+  void reunit(LatticeColorMatrix& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag);
-
   void reunit(LatticeColorMatrixD3& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag);
   
-  
-  void reunit(LatticeColorMatrixF3& xa,
+  void reunit(LatticeColorMatrix& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag);
-  
   void reunit(LatticeColorMatrixD3& xa,
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag);
   
   // With numbad and ruflag
-  void reunit(LatticeColorMatrixF3& xa, 
+  void reunit(LatticeColorMatrix& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
 	      enum Reunitarize ruflag);
-  
   void reunit(LatticeColorMatrixD3& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
 	      enum Reunitarize ruflag);
   
-  void reunit(LatticeColorMatrixF3& xa, 
+  void reunit(LatticeColorMatrix& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
 	      enum Reunitarize ruflag,
 	      const Subset& mstag);
-
   void reunit(LatticeColorMatrixD3& xa, 
 	      LatticeBoolean& bad, 
 	      int& numbad, 
diff --git a/lib/util/gauge/unit_check.cc b/lib/util/gauge/unit_check.cc
index 1be830b42..60f45ae65 100644
--- a/lib/util/gauge/unit_check.cc
+++ b/lib/util/gauge/unit_check.cc
@@ -17,7 +17,8 @@ namespace Chroma
    *
    * \param  u  The LatticeColorMatrix to be tested
    */
-  void unitarityCheck(const multi1d<LatticeColorMatrixF3>& u)
+  template<typename Q>
+  void unitarityCheck(const multi1d<Q>& u)
   {
     START_CODE();
 
@@ -25,27 +26,18 @@ namespace Chroma
 
     for (int mu=0; mu < Nd; ++mu)
     {
-      LatticeColorMatrixF3 u_tmp = u[mu];
+      Q u_tmp = u[mu];
       reunit(u_tmp, numbad, REUNITARIZE_ERROR);
     }
     
     END_CODE();
   }
 
- void unitarityCheck(const multi1d<LatticeColorMatrixD3>& u)
-  {
-    START_CODE();
-
-    int numbad;
-
-    for (int mu=0; mu < Nd; ++mu)
-    {
-      LatticeColorMatrixD3 u_tmp = u[mu];
-      reunit(u_tmp, numbad, REUNITARIZE_ERROR);
-    }
-    
-    END_CODE();
-  }
+  // Instantiate template for needed types.
+  template void unitarityCheck(const multi1d<LatticeColorMatrix>& u);
+#if DEFAULT_PRECISION!=64 && QDP_NC==3
+  template void unitarityCheck(const multi1d<LatticeColorMatrixD3>& u);
+#endif
 
 }
 
diff --git a/lib/util/gauge/unit_check.h b/lib/util/gauge/unit_check.h
index 1749b86bb..e42b50cf1 100644
--- a/lib/util/gauge/unit_check.h
+++ b/lib/util/gauge/unit_check.h
@@ -11,8 +11,8 @@ namespace Chroma
    *
    * \param  u  The multi1d<LatticeColorMatrix> to be tested
    */
-  void unitarityCheck(const multi1d<LatticeColorMatrixF3>& u);
-  void unitarityCheck(const multi1d<LatticeColorMatrixD3>& u);
+  template <typename Q>
+    void unitarityCheck(const multi1d<Q>& u);
 
 }