ENH: Solve toeplitz linear systems 

This PR implements Levinson recursion to solve toeplitz linear systems. Much of code is based on @argriffing's PR earlier PR #2754, but I used cython, and the performance is substantially improved. This is also related to #2432. The underlying cython code returns all the workspace variables, which I think is something that @josef-pkt mentioned are of interest for some timeseries applications.

```
$ python scipy/linalg/benchmarks/bench_solve_toeplitz.py

                    solve_toeplitz vs. generic solver
==============================================================
      shape      |  solver     |        dtype       |   time
                 |             |                    | (seconds)
--------------------------------------------------------------
      (100, 100) |   generic   |            float64 |   0.04
      (100, 100) |   toeplitz  |            float64 |   0.00
      (100, 100) |   generic   |         complex128 |   0.00
      (100, 100) |   toeplitz  |         complex128 |   0.00
      (300, 300) |   generic   |            float64 |   0.00
      (300, 300) |   toeplitz  |            float64 |   0.00
      (300, 300) |   generic   |         complex128 |   0.01
      (300, 300) |   toeplitz  |         complex128 |   0.00
    (1000, 1000) |   generic   |            float64 |   0.10
    (1000, 1000) |   toeplitz  |            float64 |   0.00
    (1000, 1000) |   generic   |         complex128 |   0.34
    (1000, 1000) |   toeplitz  |         complex128 |   0.01
```
Yeah, I'm not sure how to deal with that. The problem is that `r=None` is the same default as in `scipy.linalg.toeplitz`, and I wanted this to have the same signature. And you can't have a required argument after an optional argument. But I agree it's not pretty...