diff --git a/R/linkage.R b/R/linkage.R
new file mode 100644
index 0000000..5c127fd
--- /dev/null
+++ b/R/linkage.R
@@ -0,0 +1,76 @@
+#' @title Linkage
+#' @description Create hierarchical cluster tree.
+#' @details Z = LINKAGE(Y) creates a hierarchical cluster tree, using the single
+#' linkage algorithm.  The input Y is a distance matrix such as is generated by
+#' PDIST. Y may also be a more general dissimilarity matrix conforming to the
+#' output format of PDIST.
+#' @param Y data
+#' @param method either 'si', 'av', 'co' 'ce' or 'wa'
+#' @export
+linkage <- function(Y, method) {
+	k <- size(Y)[1]
+	n <- size(Y)[2]
+	m <- (1 + sqrt(1 + 8 * n)) / 2
+	if ((k != 1) | (m != trunc(m))) {
+		stop(
+			'The first input has to match the output',
+			'of the PDIST function in size.'
+		)
+	}
+	if (nargin == 1) { # set default switch to be 'co'
+		method <- 'co'
+	}
+	method <- lower(method[1:2]) # simplify the switch string.
+	monotonic <- 1
+	Z <- zeros(m - 1, 3) # allocate the output matrix.
+	N <- zeros(1, 2 * m - 1)
+	N[1:m] <- 1
+	n <- m; # since m is changing, we need to save m in n.
+	R <- 1:n
+	for (s in 1:(n-1)) {
+		X <- Y
+		v <- min(X)[1]
+		k <- min(X)[2]
+		i <- floor(m + 1 / 2 - sqrt(m ^ 2 - m + 1 / 4 - 2 * (k - 1)))
+		j <- k - (i - 1) * (m - i / 2) + i
+		Z[s, ] <- c(R[i], R[j], v) # update one more row to the output matrix A
+		# Temp variables
+		I1 <- 1:(i - 1)
+		I2 <- (i + 1):(j - 1)
+		I3 <- (j + 1):m
+
+		U <- c(I1, I2, I3)
+		I <- c(
+			I1 * (m - (I1 + 1) / 2) - m + i,
+			i * (m - (i + 1) / 2) - m + I2,
+			i * (m - (i + 1) / 2) - m + I3
+		)
+		J <- c(
+			I1 * (m - (I1 + 1) / 2) - m + j,
+			I2 * (m - (I2 + 1) / 2) - m + j,
+			j * (m - (j + 1) / 2) - m + I3
+		)
+
+		switch(method,
+			'si' = Y[I] <- min(Y[I], Y[J]), # single linkage
+			'av' = Y[I] <- Y[I] + Y[J], # average linkage
+			'co' = Y[I] <- max(Y[I], Y[J]), #complete linkage
+			'ce' = {
+				K <- N[R[i]] + N[R[j]] # centroid linkage
+				Y[I] <- (N[R[i]] * Y[I] + N[R[j]] * Y[J] -
+					(N[R[i]] * N[R[j]] * v ^ 2) / K) / K
+			},
+			'wa' = Y[I] <- ((N[R[U]] + N[R[i]]) * Y[I] + (N[R[U]] + N[R[j]]) *
+				Y[J] - N[R[U]] * v) / (N[R[i]] + N[R[j]] + N[R[U]])
+		)
+		J <- c(J, i * (m - (i + 1) / 2) - m + j)
+		Y[J] <- vector() # no need for the cluster information about j
+
+		# update m, N, R
+		m <- m - 1
+		N[n + s] <- N[R[i]] + N[R[j]]
+		R[i] <- n + s
+		R[j:(n - 1)] <- R[(j + 1):n]
+	}
+	return(Z)
+}
\ No newline at end of file
diff --git a/man/linkage.Rd b/man/linkage.Rd
new file mode 100644
index 0000000..57d5113
--- /dev/null
+++ b/man/linkage.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/linkage.R
+\name{linkage}
+\alias{linkage}
+\title{Linkage}
+\usage{
+linkage(Y, method)
+}
+\arguments{
+\item{Y}{data}
+
+\item{method}{either 'si', 'av', 'co' 'ce' or 'wa'}
+}
+\description{
+Create hierarchical cluster tree.
+}
+\details{
+Z = LINKAGE(Y) creates a hierarchical cluster tree, using the single
+linkage algorithm.  The input Y is a distance matrix such as is generated by
+PDIST. Y may also be a more general dissimilarity matrix conforming to the
+output format of PDIST.
+}