[cleanup] Move the FD approximation of the GFE function gradient from the GFE function into the test

It is a part of the test harness, and should not clutter the implementation of the
GFE functions themselves.

Use 1e-8: The square root of the machine precision is said to be the best choice.

Needed for projection-based finite elements.

This patch also adds a unit test.  So far this test passes.

The program is not really a test, but it doesn't really belong into 'src', either.

[[Imported from SVN: r9997]]

[[Imported from SVN: r9913]]

[[Imported from SVN: r9912]]

[[Imported from SVN: r9877]]

That should be no problem now that we use a 50-decimal digits type
for the FD approximation.

[[Imported from SVN: r9876]]

[[Imported from SVN: r9875]]

It offers a 'double' interface, but we don't really need that.

[[Imported from SVN: r9874]]

[[Imported from SVN: r9873]]

[[Imported from SVN: r9872]]

[[Imported from SVN: r9870]]

This code compares scalar and vector AD, and a finite difference
approximation.

[[Imported from SVN: r9869]]

[[Imported from SVN: r9681]]

[[Imported from SVN: r9680]]

[[Imported from SVN: r9614]]

[[Imported from SVN: r9573]]

[[Imported from SVN: r9436]]

This still fails -- there must still be a bug in the Hessian computation

[[Imported from SVN: r9427]]

[[Imported from SVN: r9420]]

[[Imported from SVN: r9419]]

Hence many input configurations are tested now

[[Imported from SVN: r9418]]

They are easier to debug

[[Imported from SVN: r9417]]

[[Imported from SVN: r9416]]

[[Imported from SVN: r9413]]

[[Imported from SVN: r9411]]

[[Imported from SVN: r9410]]

[[Imported from SVN: r9409]]