Gradient following method implemented.  A diagonal approximation is used for the coef-coef block of Lagrangian hessian.
This prevents the construction and factorization of normal equations, which is the most memory and CPU intensive step.
Other changes:
- Linear search algorithm reworked.  The program now use a more robust algorithm to sample along the step direction and
use linear gradient extrapolation to find better step length.  Each pair of adjacent points are chosen and the pair that
give the best predicted improvements are used.
- Fixed a flag issue which made version subroutine fail to compile with `gfortran`