- % Unfortunately, since we don't know the matrix ``C``, it isn't
- % easy to compute alpha_k with an existing step size function.
- alpha_k = (rk' * zk)/(dk' * A * dk);
- x_next = x + alpha_k*dk;
- r_next = rk + alpha_k*A*dk;
+ % Used twice, avoid recomputation.
+ rkzk = rk' * zk;
+
+ % The term alpha_k*dk appears twice, but so does Q*dk. We can't
+ % do them both, so we precompute the more expensive operation.
+ Qdk = Q * dk;
+
+ % After substituting the two previously-created variables, the
+ % following algorithm occurs verbatim in the reference.
+ alpha_k = rkzk/(dk' * Qdk);
+ x_next = xk + (alpha_k * dk);
+ r_next = rk + (alpha_k * Qdk);