Abstract. We consider the numerical solution of the generalized Lyapunov and Stein equations in Rn, arising respectively from stochastic optimal control in continuous- and discrete-time, as well as bilinear model order reduction. Generalizing the Smith method, our algorithms converge quadratically and have an O(n3) computational complexity and memory requirement. For large-scale problems, where the relevant matrix operators are “sparse”, our algorithms are of an O(n) complexity and memory requirement. Applying to Newton’s method for the solution of the rational Riccati equations (of large-scale), algorithms of an O(n3) (or O(n)) complexity and memory requirement are obtained. This contrasts favourably with the naive Newton’s method algorithms of an O(n6) complexity or the slower modified Newton’s methods of an O(n3) complexity. The methods also provide efficient alternatives to the (Bi-)CG methods with ADI preconditions by Benner and Damm (2013) for bilinear model order redu...