Add Schur and Hessenberg functions and a matrix square root.

The matrix square root uses the spectral decomposition for symmetric matrices. For non-symmetric matrices it uses the method of Björck and Hammerling (1983) which was the reference in "A new sqrtm..." that was mentioned on the list.

One peculiar thing with ```sqrtm``` is that
```julia
julia> sqrtm(-1)
0.0 + 1.0im
```
even though the usual square root returns an error. I have chosen it to do so because I find it reasonable that ```sqrtm(-ones(1,1))[1]==sqrtm(-1)```. The former is in line with the behaviour of ```eig``` for real matrices with complex eigenvalues.