POV-Ray : Newsgroups : povray.programming : Speeding up MInverse Server Time
29 Mar 2024 07:05:45 EDT (-0400)
  Speeding up MInverse (Message 1 to 4 of 4)  
From: Algo
Subject: Speeding up MInverse
Date: 27 Jun 2008 15:25:01
Message: <web.48653da6c2d04118d9b2aafa0@news.povray.org>
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

From: Warp
Subject: Re: Speeding up MInverse
Date: 27 Jun 2008 15:49:39
Message: <48654453@news.povray.org>
It would be interesting to see some actual rendertime comparisons.

-- 
                                                          - Warp


Post a reply to this message

From: Algo
Subject: Re: Speeding up MInverse
Date: 29 Jun 2008 19:05:00
Message: <web.486814be53fb1507aeb59a6c0@news.povray.org>
I would be glad to compile into a working POV render.

I need access to a development source code.

Algo


Post a reply to this message

From: Warp
Subject: Re: Speeding up MInverse
Date: 30 Jun 2008 10:24:26
Message: <4868ec9a@news.povray.org>
Algo <Gem### [at] hotmailcom> wrote:
> I would be glad to compile into a working POV render.

> I need access to a development source code.

  The source codes downloadable from povray.org are not sufficient?

-- 
                                                          - Warp


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.