diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py
index fc4a2f33..dcafeb88 100644
--- a/GPy/kern/_src/coregionalize.py
+++ b/GPy/kern/_src/coregionalize.py
@@ -6,6 +6,7 @@ import numpy as np
 from scipy import weave
 from ...core.parameterization import Param
 from ...core.parameterization.transformations import Logexp
+from ...util.config import config # for assesing whether to use weave
 
 class Coregionalize(Kern):
     """
@@ -56,6 +57,27 @@ class Coregionalize(Kern):
         self.B = np.dot(self.W, self.W.T) + np.diag(self.kappa)
 
     def K(self, X, X2=None):
+        if config.getboolean('weave', 'working'):
+            try:
+                return self._K_weave(X, X2)
+            except:
+                print "\n Weave compilation failed. Falling back to (slower) numpy implementation\n"
+                config.set('weave', 'working', False)
+                return self._K_numpy(X, X2)
+        else:
+            return self._K_numpy(X, X2)
+
+
+    def _K_numpy(self, X, X2=None):
+        index = np.asarray(X, dtype=np.int)
+        if X2 is None:
+            return self.B[index,index.T]
+        else:
+            index2 = np.asarray(X2, dtype=np.int)
+            return self.B[index,index2.T]
+
+    def _K_weave(self, X, X2=None):
+        """compute the kernel function using scipy.weave"""
         index = np.asarray(X, dtype=np.int)
 
         if X2 is None:
@@ -91,12 +113,33 @@ class Coregionalize(Kern):
 
     def update_gradients_full(self, dL_dK, X, X2=None):
         index = np.asarray(X, dtype=np.int)
-        dL_dK_small = np.zeros_like(self.B)
         if X2 is None:
             index2 = index
         else:
             index2 = np.asarray(X2, dtype=np.int)
 
+        #attempt to use weave for a nasty double indexing loop: fall back to numpy
+        if config.getboolean('weave', 'working'):
+            try:
+                dL_dK_small = self._gradient_reduce_weave(dL_dK, index1, index2)
+            except:
+                print "\n Weave compilation failed. Falling back to (slower) numpy implementation\n"
+                config.set('weave', 'working', False)
+                dL_dK_small = self._gradient_reduce_weave(dL_dK, index1, index2)
+        else:
+            dL_dK_small = self._gradient_reduce_weave(dL_dK, index1, index2)
+
+
+
+        dkappa = np.diag(dL_dK_small)
+        dL_dK_small += dL_dK_small.T
+        dW = (self.W[:, None, :]*dL_dK_small[:, :, None]).sum(0)
+
+        self.W.gradient = dW
+        self.kappa.gradient = dkappa
+
+    def _gradient_reduce_weave(self, dL_dK, index, index2):
+        dL_dK_small = np.zeros_like(self.B)
         code="""
         for(int i=0; i<num_inducing; i++){
           for(int j=0; j<N; j++){
@@ -106,13 +149,19 @@ class Coregionalize(Kern):
         """
         N, num_inducing, output_dim = index.size, index2.size, self.output_dim
         weave.inline(code, ['N', 'num_inducing', 'output_dim', 'dL_dK', 'dL_dK_small', 'index', 'index2'])
+        return dL_dK_small
+1
+    def _gradient_reduce_numpy(self, dL_dK, index, index2):
+        index, index2 = index[:,0], index2[:,0]
+        for i in range(k.output_dim):
+        dL_dK_small = np.zeros_like(self.B)
+            tmp1 = dL_dK[index==i]
+            for j in range(k.output_dim):
+                dL_dK_small[j,i] = tmp1[:,index2==j].sum()
+        return dL_dK_small
+
 
-        dkappa = np.diag(dL_dK_small)
-        dL_dK_small += dL_dK_small.T
-        dW = (self.W[:, None, :]*dL_dK_small[:, :, None]).sum(0)
 
-        self.W.gradient = dW
-        self.kappa.gradient = dkappa
 
     def update_gradients_diag(self, dL_dKdiag, X):
         index = np.asarray(X, dtype=np.int).flatten()
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 95ad7961..cde024f6 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -356,6 +356,50 @@ class KernelTestsNonContinuous(unittest.TestCase):
         X2 = self.X2[self.X2[:,-1]!=2]
         self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
 
+class Coregionalize_weave_test(unittest.TestCase):
+    """
+    Make sure that the coregionalize kernel work with and without weave enabled
+    """
+    k = GPy.kern.coregionalize(1, output_dim=12)
+    N1, N2 = 100, 200
+    X = np.random.randint(0,12,(N1,1))
+    X2 = np.random.randint(0,12,(N2,1))
+
+    #symmetric case
+    dL_dK = np.random.randn(N1, N1)
+    GPy.util.config.config.set('weave', 'working', True)
+    K_weave = k.K(X)
+    k.update_gradients_full(dL_dK, X)
+    grads_weave = k.gradient.copy()
+
+    GPy.util.config.config.set('weave', 'working', False)
+    K_numpy = k.K(X)
+    k.update_gradients_full(dL_dK, X)
+    grads_numpy = k.gradient.copy()
+
+    self.assertTrue(np.allclose(K_numpy, K_weave))
+    self.assertTrue(np.allclose(grads_numpy, grads_weave))
+
+    #non-symmetric case
+    dL_dK = np.random.randn(N1, N2)
+    GPy.util.config.config.set('weave', 'working', True)
+    K_weave = k.K(X, X2)
+    k.update_gradients_full(dL_dK, X, X2)
+    grads_weave = k.gradient.copy()
+
+    GPy.util.config.config.set('weave', 'working', False)
+    K_numpy = k.K(X, X2)
+    k.update_gradients_full(dL_dK, X, X2)
+    grads_numpy = k.gradient.copy()
+
+    self.assertTrue(np.allclose(K_numpy, K_weave))
+    self.assertTrue(np.allclose(grads_numpy, grads_weave))
+
+    #reset the weave state for any other tests
+    GPy.util.config.config.set('weave', 'working', False)
+
+
+
 
 if __name__ == "__main__":
     print "Running unit tests, please be (very) patient..."