Two fixes + one improvement to lyapunov_symm.m

Fix 1) When method=3 failed, it was using the convergence tolerance
(i.e.. a small number <<1) as the qz_criterium (which should be around
1).
Fix 2) The method failed with sparse matrices. At least the iterative
method ought to work with sparse matrices.
Improvement 1) The iterative algorithm used here for method=3 is
potentially inefficient. A very simple alternative is the doubling
algorithm described here:
http://www.math.cts.nthu.edu.tw/download.php?filename=683_a93bf34d.pdf&dir=publish&title=prep2012-05-005
To see the difference, try:

tmp=sprandn(10,10,0.1);tmp=0.99*tmp/max(abs(eigs(tmp)));tmp2=sprandn(10,10,0.2);tmp2=tmp2*tmp2';
lyapunov_symm( tmp, tmp2, sqrt(eps), 1e-15, 3 )

with line 70/71 commented/uncommented.

In my test it reduced the number of iterations from 1178 to 12, and
convergence was exact (evol=0) rather than just satisfying the bound.

The slight downside is that this requires one extra matrix
multiplication per iteration, which increases the cost per iteration by
around 50%. But 12+50% is still << 1178!