"cancor.mcd" <- 
function(x, y, covlist, covariance = FALSE, norm = 1)
{
#MCD-based canonical correlation analysis
#
# x ... matrix (data frame) of first variable set
# y ... matrix (data frame) of second variable set
# covlist ... output of covMcd (library(robustbase))
# covarance ... FALSE if the analysis should be based on correlations
# norm ... 1 is canonical vectors should be normed to length 1

	dx <- dim(x)
	dy <- dim(y)
	r <- min(dx[2], dy[2])
	n <- dx[1]
	if(n != dy[1])
		stop("x and y must have same no. of rows")
	z <- cbind(x, y)
	s <- covlist$cov
	if(covariance == FALSE) {
		dd <- diag(diag(s))
		d <- solve(dd)
		corr <- sqrt(d) %*% s %*% sqrt(d)
		s <- corr
	}
	sxx <- s[1:dx[2], 1:dx[2]]
	syy <- s[(dx[2] + 1):(dx[2] + dy[2]), (dx[2] + 1):(dx[2] + dy[2])]
	sxy <- s[1:dx[2], (dx[2] + 1):(dx[2] + dy[2])]
	syx <- s[(dx[2] + 1):(dx[2] + dy[2]), 1:dx[2]]
	require(MASS)
	invsxx <- ginv(sxx)
	invsyy <- ginv(syy)
	sy <- invsyy %*% syx %*% invsxx %*% sxy
	cancorr <- sqrt(eigen(sy)$values)
	canvecty <- eigen(sy)$vectors
	canvecty <- as.matrix(canvecty)
	dimnames(canvecty) <- list(dimnames(y)[[2]],NULL)
	if(norm == 1) {
		normalisation <- t(canvecty) %*% syy %*% canvecty
		canvecty <- t(canvecty)/sqrt(diag(normalisation))
		canvecty <- t(canvecty)
	}
	sx <- invsxx %*% sxy %*% invsyy %*% syx
	cancorr <- sqrt(eigen(sx)$values)
	canvectx <- eigen(sx)$vectors
	canvectx <- as.matrix(canvectx)
	dimnames(canvectx) <- list(dimnames(x)[[2]],NULL)
	if(norm == 1) {
		normalisation <- t(canvectx) %*% sxx %*% canvectx
		canvectx <- t(canvectx)/sqrt(diag(normalisation))
		canvectx <- t(canvectx)
	}
	list(cor=cancorr,xcoef=canvectx,ycoef=canvecty,correlation=s)
}