diff --git a/sklearn/gaussian_process/gaussian_process.py b/sklearn/gaussian_process/gaussian_process.py
index 8d2c3adb107a3af85bd3a24b830daae4baa482b1..7f7b1dc1e6ef8b108fd8401b64e82cdc4b0c9b47 100644
--- a/sklearn/gaussian_process/gaussian_process.py
+++ b/sklearn/gaussian_process/gaussian_process.py
@@ -414,7 +414,6 @@ class GaussianProcess(BaseEstimator, RegressorMixin):
 
             # Get pairwise componentwise L1-distances to the input training set
             dx = manhattan_distances(X, Y=self.X, sum_over_features=False)
-
             # Get regression function and correlation
             f = self.regr(X)
             r = self.corr(self.theta, dx).reshape(n_eval, n_samples)
diff --git a/sklearn/metrics/pairwise.py b/sklearn/metrics/pairwise.py
index d9500781e5cd9b0415a64fbf5f6a063680968c6b..c9c3301a58ebef2ecb0b49af0df91b25089eb40a 100644
--- a/sklearn/metrics/pairwise.py
+++ b/sklearn/metrics/pairwise.py
@@ -196,12 +196,14 @@ def manhattan_distances(X, Y=None, sum_over_features=True):
     array([[1]])
     >>> manhattan_distances(2, 3)
     array([[1]])
+    >>> manhattan_distances([[1, 2], [3, 4]], [[1, 2], [0, 3]])
+    array([[0, 2], [4, 4]])
     >>> import numpy as np
     >>> X = np.ones((1, 2))
     >>> y = 2 * np.ones((2, 2))
     >>> manhattan_distances(X, y, sum_over_features=False)
-    array([[[ 1.,  1.],
-            [ 1.,  1.]]])
+    array([[ 1.,  1.],
+           [ 1.,  1.]])
     """
     X, Y = check_pairwise_arrays(X, Y)
     n_samples_X, n_features_X = X.shape
@@ -211,6 +213,8 @@ def manhattan_distances(X, Y=None, sum_over_features=True):
     D = np.abs(X[:, np.newaxis, :] - Y[np.newaxis, :, :])
     if sum_over_features:
         D = np.sum(D, axis=2)
+    else:
+        D = D.reshape((n_samples_X * n_samples_Y, n_features_X))
     return D