# Bug 6748

Summary: Product: __gluInvertMatrixd:This function has error Mesa shan hao bo GLU mesa-dev RESOLVED FIXED normal medium dcm unspecified x86 (IA32) Windows (All) Test program demonstrating instability of Shan Hao Bo's algorithm

 shan hao bo 2006-04-26 13:00:03 UTC ```\$(mesasrc)/glu/sgi/libutil/project.c:__gluInvertMatrixd:This algorithm is wrong The inverse matrix of identity matrix is itself! but old __gluInvertMatrixd consider identity matrix is bad matrix. modified: /* ** inverse = invert(src) */ static int __gluInvertMatrixd(const GLdouble src[16], GLdouble inverse[16]) { int i, j, k; double t; GLdouble temp[4][4]; for (i=0; i<4; i++) { for (j=0; j<4; j++) { temp[i][j] = src[i*4+j]; } } __gluMakeIdentityd(inverse); for (i = 0; i < 4; i++) { if (temp[i][i] == 0.0f) { /* ** Look for non-zero element in column */ for (j = i + 1; j < 4; j++) { if (temp[j][i] != 0.0f) { break } } if (j != 4) { /* ** Swap rows. */ for (k = 0; k < 4; k++) { t = temp[i][k]; temp[i][k] = temp[j][k]; temp[j][k] = t; t = inverse[i*4+k]; inverse[i*4+k] = inverse[j*4+k]; inverse[j*4+k] = t; } } else { /* ** No non-zero pivot. The matrix is singular, which shouldn't ** happen. This means the user gave us a bad matrix. */ return GL_FALSE; } } t = 1.0f / temp[i][i]; for (k = 0; k < 4; k++) { temp[i][k] *= t; inverse[i*4+k] *= t; } for (j = 0; j < 4; j++) { if (j != i) { t = temp[j][i]; for (k = 0; k < 4; k++) { temp[j][k] -= temp[i][k]*t; inverse[j*4+k] -= inverse[i*4+k]*t; } } } } return GL_TRUE; }``` Michel DÃ¤nzer 2006-04-26 16:17:16 UTC `*** Bug 6749 has been marked as a duplicate of this bug. ***` Brian Paul 2006-04-27 01:36:12 UTC ```This looks like a fairly significant change to the function. Have you thoroughly exercised this change?``` shan hao bo 2006-04-27 13:48:09 UTC ```(In reply to comment #2) > This looks like a fairly significant change to the function. Have you > thoroughly exercised this change? Old function has not wrong. I have made a mistake while optimizing this function. Very sorry. But new function faster than old function; ``` shan hao bo 2006-04-27 13:50:04 UTC ```(In reply to comment #2) > This looks like a fairly significant change to the function. Have you > thoroughly exercised this change? #include #include #include //============================================================================/ / //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FOR TEST >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>// //============================================================================/ / typedef double GLdouble; #undef GL_FALSE #define GL_FALSE 0 #undef GL_TRUE #define GL_TRUE 1 __inline__ unsigned long long int rdtsc() { unsigned long long int counter; __asm__ volatile (".byte 0x0f, 0x31" : "=A" (counter)); return counter; } void randMatrix(GLdouble m[16]) { m[0+4*0] = rand() * 1.0f; m[0+4*1] = rand() * 1.0f; m[0+4*2] = rand() * 1.0f; m[0+4*3] = rand() * 1.0f; m[1+4*0] = rand() * 1.0f; m[1+4*1] = rand() * 1.0f; m[1+4*2] = rand() * 1.0f; m[1+4*3] = rand() * 1.0f; m[2+4*0] = rand() * 1.0f; m[2+4*1] = rand() * 1.0f; m[2+4*2] = rand() * 1.0f; m[2+4*3] = rand() * 1.0f; m[3+4*0] = rand() * 1.0f; m[3+4*1] = rand() * 1.0f; m[3+4*2] = rand() * 1.0f; m[3+4*3] = rand() * 1.0f; } void printMatrix(GLdouble m[16]) { printf("%f, %f, %f, %f\n", m[0+4*0], m[0+4*1], m[0+4*2], m[0+4*3]); printf("%f, %f, %f, %f\n", m[1+4*0], m[1+4*1], m[1+4*2], m[1+4*3]); printf("%f, %f, %f, %f\n", m[2+4*0], m[2+4*1], m[2+4*2], m[2+4*3]); printf("%f, %f, %f, %f\n\n", m[3+4*0], m[3+4*1], m[3+4*2], m[3+4*3]); } //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FOR TEST >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>// /// copy from project.c /* ** Make m an identity matrix */ static void __gluMakeIdentityd(GLdouble m[16]) { m[0+4*0] = 1; m[0+4*1] = 0; m[0+4*2] = 0; m[0+4*3] = 0; m[1+4*0] = 0; m[1+4*1] = 1; m[1+4*2] = 0; m[1+4*3] = 0; m[2+4*0] = 0; m[2+4*1] = 0; m[2+4*2] = 1; m[2+4*3] = 0; m[3+4*0] = 0; m[3+4*1] = 0; m[3+4*2] = 0; m[3+4*3] = 1; } /// copy from project.c and change function name /* ** inverse = invert(src) */ static int old__gluInvertMatrixd(const GLdouble src[16], GLdouble inverse[16]) { int i, j, k, swap; double t; GLdouble temp[4][4]; for (i=0; i<4; i++) { for (j=0; j<4; j++) { temp[i][j] = src[i*4+j]; } } __gluMakeIdentityd(inverse); for (i = 0; i < 4; i++) { /* ** Look for largest element in column */ swap = i; for (j = i + 1; j < 4; j++) { if (fabs(temp[j][i]) > fabs(temp[i][i])) { swap = j; } } if (swap != i) { /* ** Swap rows. */ for (k = 0; k < 4; k++) { t = temp[i][k]; temp[i][k] = temp[swap][k]; temp[swap][k] = t; t = inverse[i*4+k]; inverse[i*4+k] = inverse[swap*4+k]; inverse[swap*4+k] = t; } } if (temp[i][i] == 0) { /* ** No non-zero pivot. The matrix is singular, which shouldn't ** happen. This means the user gave us a bad matrix. */ return GL_FALSE; } t = temp[i][i]; for (k = 0; k < 4; k++) { temp[i][k] /= t; inverse[i*4+k] /= t; } for (j = 0; j < 4; j++) { if (j != i) { t = temp[j][i]; for (k = 0; k < 4; k++) { temp[j][k] -= temp[i][k]*t; inverse[j*4+k] -= inverse[i*4+k]*t; } } } } return GL_TRUE; } /// copy from project.c static void __gluMultMatricesd(const GLdouble a[16], const GLdouble b[16], GLdouble r[16]) { int i, j; for (i = 0; i < 4; i++) { for (j = 0; j < 4; j++) { r[i*4+j] = a[i*4+0]*b[0*4+j] + a[i*4+1]*b[1*4+j] + a[i*4+2]*b[2*4+j] + a[i*4+3]*b[3*4+j]; } } } ///-------------------------------------------------------------------------- /// ///-------------------------------------------------------------------------- /// ///-------------------------------------------------------------------------- /// /// new __gluInvertMatrixd /* ** inverse = invert(src) */ static int __gluInvertMatrixd(const GLdouble src[16], GLdouble inverse[16]) { int i, j, k; double t; GLdouble temp[4][4]; for (i=0; i<4; i++) { for (j=0; j<4; j++) { temp[i][j] = src[i*4+j]; } } __gluMakeIdentityd(inverse); for (i = 0; i < 4; i++) { if (temp[i][i] == 0.0f) { /* ** Look for non-zero element in column */ for (j = i + 1; j < 4; j++) { if (temp[j][i] != 0.0f) { break; } } if (j != 4) { /* ** Swap rows. */ for (k = 0; k < 4; k++) { t = temp[i][k]; temp[i][k] = temp[j][k]; temp[j][k] = t; t = inverse[i*4+k]; inverse[i*4+k] = inverse[j*4+k]; inverse[j*4+k] = t; } } else { /* ** No non-zero pivot. The matrix is singular, which shouldn't ** happen. This means the user gave us a bad matrix. */ return GL_FALSE; } } t = 1.0f / temp[i][i]; for (k = 0; k < 4; k++) { temp[i][k] *= t; inverse[i*4+k] *= t; } for (j = 0; j < 4; j++) { if (j != i) { t = temp[j][i]; for (k = 0; k < 4; k++) { temp[j][k] -= temp[i][k]*t; inverse[j*4+k] -= inverse[i*4+k]*t; } } } } return GL_TRUE; } ///-------------------------------------------------------------------------- /// ///-------------------------------------------------------------------------- /// ///-------------------------------------------------------------------------- /// int main(int argc, char *argv[]) { srand(time(0)); GLdouble mat[16], inv[16], r[16]; printf("Test identity matrix:\n"); __gluMakeIdentityd(mat); printMatrix(mat); printf("Using old function:\n"); if (old__gluInvertMatrixd(mat, inv) == GL_FALSE) { printf("this is singular matrix\n"); } else { printMatrix(inv); } printf("Using new function:\n"); if (__gluInvertMatrixd(mat, inv) == GL_FALSE) { printf("this is singular matrix\n"); } else { printMatrix(inv); } unsigned long long int x; unsigned long long int oc = 0; unsigned long long int nc = 0; int i = 0; for (i = 0; i < 100; i++) { x = rdtsc(); old__gluInvertMatrixd(mat, inv); oc += (rdtsc() - x); x = rdtsc(); __gluInvertMatrixd(mat, inv); nc += (rdtsc() - x); } printf("\nold function:%d\n", oc); printf("\nnew function:%d\n", nc); printf("\nrate %f\n", (double)(nc)/(double)(oc)); return 0; } ``` Brian Paul 2006-05-02 02:02:24 UTC `OK, I've checked in the new function. I did some tests and it looks good. Thanks.` David Moore 2007-08-26 22:47:32 UTC `I'd like to reopen this bug if possible. The new version of __gluInvertMatrixd() that has been in Mesa since 2006 seems to have some numerical instabilities. Specifically, it will only pivot when a diagonal element equals zero. However, if the value is vanishingly small but not quite zero, no pivot will be performed, and the result of the inverse can blow up in some cases. I have observed this in real applications. I will try to find an example and post it.` David Moore 2007-08-26 23:57:57 UTC ```Created attachment 11287 [details] Test program demonstrating instability of Shan Hao Bo's algorithm This code demonstrates a matrix inverse on a matrix with some near-zero elements. The result is shown for three inverse functions: 1) The pre-2006 Mesa implementation 2) The post-2006 implementation from Shan Hao Bo 3) A proposed replacement for both from me. As you can see, methods 1 and 3 agree on the answer, but 2 differs. It also shows timing characteristics for all three. Here is the output when run on my machine: Starting with matrix: -0.000000, -1.420000, 0.000000, 0.000000 1.710000, -0.000000, 0.302000, -0.000000 0.174000, -0.000000, -0.985000, 49.800000 0.174000, -0.000000, -0.985000, 50.000000 Inverse using original pre-2006 Mesa function: 0.000000, 0.567103, 43.468298, -43.294425 -0.704225, -0.000000, -0.000000, 0.000000 0.000000, 0.100179, -246.128443, 245.143929 0.000000, 0.000000, -5.000000, 5.000000 Inverse using Shan Hao Bo function (2006 onward): -4.000000, 0.567103, 43.468298, -43.294425 -0.704225, -0.000000, -0.000000, 0.000000 0.000000, 0.100179, -246.128443, 245.143929 0.000000, 0.000000, -5.000000, 5.000000 Inverse using Moore's proposed function: -0.000000, 0.567103, 43.468298, -43.294425 -0.704225, -0.000000, -0.000000, 0.000000 0.000000, 0.100179, -246.128443, 245.143929 0.000000, -0.000000, -5.000000, 5.000000 old function: 59094 ticks Shan Hao Bo function: 51597 ticks Moore function: 32796 ticks``` Brian Paul 2007-08-27 09:38:10 UTC ```Thanks, David. I've checked in your new code. I plugged it into core Mesa for testing and things ran as expected, though it seems a bit slower than Mesa's present matrix inversion code (which is special-cased for various types of matrices). ```

Use of freedesktop.org services, including Bugzilla, is subject to our Code of Conduct. How we collect and use information is described in our Privacy Policy.