norm, covar, dlyap cannot handle systems of dimension > 150

```julia
N = 200
dlyap(randn(N,N), Matrix(I(N)))
```
crashes julia (on my machine) due to the naive implementation of the Lyapunov solver that attempts to invert N^2 x N^2 matrices.

`covar` and `norm` rely on `dlyap` so these functions have the same problem.

@mfalt , @baggepinnen we should discuss how to handle this and related problems