|
 |
Excessive dependence upon the dividers (floating-point or integer) results in
many lost opportunities in scientific programming. I have a joke about this.
In this case, the POV code, MInverse, calls for 16 FPDIVs when only 1 is
actually necessary. The easily revised code, NEW_MInverse, is below.
I have tested the change for accuracy by comparing the inversion of the
identical test case using both methods. The differences in their answers are
within what would be expected from variations in the least-significant bit of
double precision arithmetic.
RMS difference between NEW and POV = 2.36E-013
To test the absolute accuracy of NEW_MInverse, I simply compute the product of
the computed inverses matrix-multiplied by the original test matrix. The
result should be the identity matrix, of course. I am unable to detect
systematic differences in accuracy between the two methods BECAUSE THE ERROR IN
THIS ALGORITHM IS SO LARGE:
RMS Error POV: 3.13E-009
RMS Error NEW: 3.13E-009
The inverses computed (by either of these methods) are more like
single-precision than double. It is quite possible that these surprisingly
extreme, it is conceivable that POV has actually been a single-precision code
all these years without our knowing it.
It is very easy to write highly accurate matrix inversions, though the result
could never rival MInverse on the basis of elegance and transparency. It would
be a shame to do so unless there were a real reason.
Fortunately, this is a scientific project, so that such worries are not some
vague feeling: they reduce to a measured number. Why we like scientific
projects. What would be the symptoms of inaccuracy in MInverse? And do we see
these symptoms anywhere?
On the other hand, there can be little debate about the advantage of
NEW_MInverse, which has zero downside and is a good deal faster:
POV MInverse: 2.41E-007(sec)
NEW_MInverse: 1.47E-007(sec)
The above was measured with Release version of VC++2008.
--Algo
/////////////////////////////////////////////////////////////////////
////This is ready to plug into POV-3.6///////////////////////////////
/////////////////////////////////////////////////////////////////////
void NEW_MInvers(MATRIX r, MATRIX m)
{
DBL d00, d01, d02, d03;
DBL d10, d11, d12, d13;
DBL d20, d21, d22, d23;
DBL d30, d31, d32, d33;
DBL m00, m01, m02, m03;
DBL m10, m11, m12, m13;
DBL m20, m21, m22, m23;
DBL m30, m31, m32, m33;
DBL D;
m00 = m[0][0]; m01 = m[0][1]; m02 = m[0][2]; m03 = m[0][3];
m10 = m[1][0]; m11 = m[1][1]; m12 = m[1][2]; m13 = m[1][3];
m20 = m[2][0]; m21 = m[2][1]; m22 = m[2][2]; m23 = m[2][3];
m30 = m[3][0]; m31 = m[3][1]; m32 = m[3][2]; m33 = m[3][3];
d00 = m11*m22*m33 + m12*m23*m31 + m13*m21*m32 - m31*m22*m13 - m32*m23*m11 -
m33*m21*m12;
d01 = m10*m22*m33 + m12*m23*m30 + m13*m20*m32 - m30*m22*m13 - m32*m23*m10 -
m33*m20*m12;
d02 = m10*m21*m33 + m11*m23*m30 + m13*m20*m31 - m30*m21*m13 - m31*m23*m10 -
m33*m20*m11;
d03 = m10*m21*m32 + m11*m22*m30 + m12*m20*m31 - m30*m21*m12 - m31*m22*m10 -
m32*m20*m11;
d10 = m01*m22*m33 + m02*m23*m31 + m03*m21*m32 - m31*m22*m03 - m32*m23*m01 -
m33*m21*m02;
d11 = m00*m22*m33 + m02*m23*m30 + m03*m20*m32 - m30*m22*m03 - m32*m23*m00 -
m33*m20*m02;
d12 = m00*m21*m33 + m01*m23*m30 + m03*m20*m31 - m30*m21*m03 - m31*m23*m00 -
m33*m20*m01;
d13 = m00*m21*m32 + m01*m22*m30 + m02*m20*m31 - m30*m21*m02 - m31*m22*m00 -
m32*m20*m01;
d20 = m01*m12*m33 + m02*m13*m31 + m03*m11*m32 - m31*m12*m03 - m32*m13*m01 -
m33*m11*m02;
d21 = m00*m12*m33 + m02*m13*m30 + m03*m10*m32 - m30*m12*m03 - m32*m13*m00 -
m33*m10*m02;
d22 = m00*m11*m33 + m01*m13*m30 + m03*m10*m31 - m30*m11*m03 - m31*m13*m00 -
m33*m10*m01;
d23 = m00*m11*m32 + m01*m12*m30 + m02*m10*m31 - m30*m11*m02 - m31*m12*m00 -
m32*m10*m01;
d30 = m01*m12*m23 + m02*m13*m21 + m03*m11*m22 - m21*m12*m03 - m22*m13*m01 -
m23*m11*m02;
d31 = m00*m12*m23 + m02*m13*m20 + m03*m10*m22 - m20*m12*m03 - m22*m13*m00 -
m23*m10*m02;
d32 = m00*m11*m23 + m01*m13*m20 + m03*m10*m21 - m20*m11*m03 - m21*m13*m00 -
m23*m10*m01;
d33 = m00*m11*m22 + m01*m12*m20 + m02*m10*m21 - m20*m11*m02 - m21*m12*m00 -
m22*m10*m01;
D = m00*d00 - m01*d01 + m02*d02 - m03*d03;
if (D == 0.0)
{
Error("Singular matrix in MInvers.");
}
DBL InvD = 1.0/D;
r[0][0] = d00*InvD; r[0][1] = -d10*InvD; r[0][2] = d20*InvD; r[0][3] =
-d30*InvD;
r[1][0] = -d01*InvD; r[1][1] = d11*InvD; r[1][2] = -d21*InvD; r[1][3] =
d31*InvD;
r[2][0] = d02*InvD; r[2][1] = -d12*InvD; r[2][2] = d22*InvD; r[2][3] =
-d32*InvD;
r[3][0] = -d03*InvD; r[3][1] = d13*InvD; r[3][2] = -d23*InvD; r[3][3] =
d33*InvD;
}
/////////////////////////////////////////////////////////////////////
Post a reply to this message
|
|
com> wrote: