Why does k-means clustering algorithm use only Euclidean distance metric? Is there a specific purpose in terms of efficiency or functionality why the k-means algorithm does not use for example cosine (dis)similarity as a distance metric, but can only use the Euclidean norm? In general, will K-means method comply and be correct when other distances than Euclidean are considered or used?
[Addition by @ttnphns. The question is two-fold. "(Non)Euclidean distance" may concern distance between two data points or distance between a data point and a cluster centre. Both ways have been attempted to address in the answers so far.]
 A: K-Means procedure - which is a vector quantization method often used as a clustering method - does not explicitly use pairwise distances between data points at all (in contrast to hierarchical and some other clusterings which allow for arbitrary proximity measure). It amounts to repeatedly assigning points to the closest centroid thereby using Euclidean distance from data points to a centroid. However, K-Means is implicitly based on pairwise Euclidean distances between data points, because the sum of squared deviations from centroid is equal to the sum of pairwise squared Euclidean distances divided by the number of points. The term "centroid" is itself from Euclidean geometry. It is multivariate mean in euclidean space. Euclidean space is about euclidean distances. Non-Euclidean distances will generally not span Euclidean space. That's why K-Means is for Euclidean distances only.
But a Euclidean distance between two data points can be represented in a number of alternative ways. For example, it is closely tied with cosine or scalar product between the points. If you have cosine, or covariance, or correlation, you can always (1) transform it to (squared) Euclidean distance, and then (2) create data for that matrix of Euclidean distances (by means of Principal Coordinates or other forms of metric Multidimensional Scaling) to (3) input those data to K-Means clustering. Therefore, it is possible to make K-Means "work with" pairwise cosines or such; in fact, such implementations of K-Means clustering exist. See also about "K-means for distance matrix" implementation.
It is possible to program K-means in a way that it directly calculate on the square matrix of pairwise Euclidean distances, of course. But it will work slowly, and so the more efficient way is to create data for that distance matrix (converting the distances into scalar products and so on - the pass that is outlined in the previous paragraph) - and then apply standard K-means procedure to that dataset.
Please note I was discussing the topic whether euclidean or noneuclidean dissimilarity between data points is compatible with K-means. It is related to but not quite the same question as whether noneuclidean deviations from centroid (in wide sense, centre or quasicentroid) can be incorporated in K-means or modified "K-means".
See related question K-means: Why minimizing WCSS is maximizing Distance between clusters?.
A: Since this is apparently now a canonical question, and it hasn't been mentioned here yet:
One natural extension of k-means to use distance metrics other than the standard Euclidean distance on $\mathbb R^d$ is to use the kernel trick. This refers to the idea of implicitly mapping the inputs to a high-, or infinite-, dimensional Hilbert space, where distances correspond to the distance function we want to use, and run the algorithm there. That is, letting $\varphi : \mathbb R^p \to \mathcal H$ be some feature map such that the desired metric $d$ can be written $d(x, y) = \lVert \varphi(x) - \varphi(y) \rVert_{\mathcal H}$, we run k-means on the points $\{ \varphi(x_i) \}$. In many cases, we can't compute the map $\varphi$ explicitly, but we can compute the kernel $k(x, y) = \langle \varphi(x), \varphi(y) \rangle_{\mathcal H}$. Not all distance metrics fit this model, but many do, and there are such functions defined on strings, graphs, images, probability distributions, and more....
In this situation, in the standard (Lloyd's) k-means algorithm, we can assign easily points to their clusters, but we represent the cluster centers implicitly (as linear combinations of the input points in Hilbert space). Finding the best representation in the input space would require finding a Fréchet mean, which is quite expensive. So it's easy to get cluster assignments with a kernel, harder to get the means.
The following paper discusses this algorithm, and relates it to spectral clustering:

I. Dhillon, Y. Guan, and B. Kulis. Kernel k-means, Spectral Clustering and Normalized Cuts. KDD 2005.

A: See also @ttnphns answer for an interpretation of k-means that actually involves pointwise Euclidean distances.
The way k-means is constructed is not based on distances.
K-means minimizes within-cluster variance. Now if you look at the definition of variance, it is identical to the sum of squared Euclidean distances from the center. (@ttnphns answer refers to pairwise Euclidean distances!)
The basic idea of k-means is to minimize squared errors. There is no "distance" involved here.
Why it is not correct to use arbitary distances: because k-means may stop converging with other distance functions. The common proof of convergence is like this: the assignment step and the mean update step both optimize the same criterion. There is a finite number of assignments possible. Therefore, it must converge after a finite number of improvements.
To use this proof for other distance functions, you must show that the mean (note: k-means) minimizes your distances, too.
If you are looking for an Manhattan-distance variant of k-means, there is k-medians.
Because the median is a known best L1 estimator.
If you want arbitrary distance functions, have a look at k-medoids (aka: PAM, partitioning around medoids). The medoid minimizes arbitrary distances (because it is defined as the minimum), and there only exist a finite number of possible medoids, too. It is much more expensive than the mean, though.
A: I've read many interesting comments here, but let me add that Matlab's "personal" implementation of k-means supports 4 non-Euclidean distances [between data points and cluster centres]. The only comment from the documentation I can see about that is:

Distance measure, in p-dimensional space, used for minimization, specified as the comma-separated pair consisting of 'Distance' and a string.
kmeans computes centroid clusters differently for the different, supported distance measures. This table summarizes the available distance measures. In the formulae, x is an observation (that is, a row of X) and c is a centroid (a row vector).

Then a list of functions of c and x follows. Thus, considering that p is the dimensionality of the input data, it seems that no Euclidean embedding is performed beforehand.
BTW in the past I've been using Matlab's k-means with correlation distance and it (unsurprisingly) did what it was supposed to do.
A: From here:


Let us consider two documents A and B represented by the vectors in the above figure. The cosine treats both vectors as unit vectors by normalizing them, giving you a measure of the angle between the two vectors. It does provide an accurate measure of similarity but with no regard to magnitude. But magnitude is an important factor while considering similarity.

A: I might be a little pedantic here,  but K-means is the name given to a particular algorithm that assigns labels to data points such that within cluster variances are minimized, and it is not the name for a "general technique". 
K-means algorithm has been independently proposed from several fields, with strong interpretations applicable to the field. It just turns out, nicely, that it is also euclidean distance to the center. For a brief history of K-means, please read Data Clustering: 50-years beyond K-means
There are a plethora of other clustering algorithms that use metrics other than Euclidean. The most general case I know is of using Bregman Divergences for clustering, of which Euclidean is a special case.
