Precision problem

["Matt Hauer" <matthew dot hauer at gmail dot com>]

I've been beating my head against my implementation of a Kalman
filter, and tracked it down to the residual being ill-conditioned.  I
have the same thing working in Matlab, and it's doing fine.  I've
tracked the problem down to calculation inaccuracies perturbing the
residual and making it worse conditioned in Fortran than it is in
Matlab, to the tune of from 10^8 to 10^19 or so.

I'm reading some values out of a file written by Matlab with a mess of
8-byte doubles.  Pretty standard binary file.  Those values seem to be
okay, base on variables like Pm.

(caveat: Matlab prints 14 digits and Fortran 16.  log_10(2^52) = 15.6.
 Matlab appears to round the last two off for clarity.)

Ie, Pm (straight from files) =
M - {1.06570664344110,   1.63000000000000,   0.850000000000000}
F - {1.0657066434411027, 1.6299999999999988, 0.84999999999999998}

"Close enough" to be called equally precise.  Now, doing some simple
operations doesn't throw the numbers off very far, as I found with
calculating E.

(eprime is read from the file)
E = DBLE(ZABS(eprime(:,1)))

E =
M - {1.08329039771713,   1.10643027162989,   1.07690305534670}
F - {1.0832903977171333, 1.1064302716298939, 1.0769030553467029}

Trying to calculate V, where delta is hardcoded (with a d0 suffix),
RecV_t is from the file, and E is from above...
E_gen = DCMPLX(E * DCOS(delta), E * DSIN(delta)) theta_a =
MATMUL(RecV_t, E_gen) !9x3 * 3x1, both complex*16
DO i=1, nB
 V(i) = DBLE(ZABS(theta_a(i)))
END DO

V  =
M - {1.04000000000000,   1.02533000000000,   1.02536000000000}
F - {1.0400000064889618, 1.0253299314257127, 1.0253599693910453}
Dif {6.489e-9, 6.857e-8, 3.061e-8}

At this point, I can only count up to six? multiplications and one
sqrt between calculating E and calculating V, which shouldn't lose
more than two decimals of precision.  My best guess is that DCOS or
DSIN are introducing an unexpected degree of error.  I've been going
over this for a couple days now, and become increasingly exasperated.
Does anyone have some insight on what sort of precision errors are
introduced by using the Fortran trig functions?

For reference, I've dropped this code through both g95 and gfortran,
both of which I've updated to make sure that there aren't any known
bugs lurking in them.