Bug 6748

*** Bug 6749 has been marked as a duplicate of this bug. ***

This looks like a fairly significant change to the function.  Have you
thoroughly exercised this change?

(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;

(In reply to comment #2)
> This looks like a fairly significant change to the function.  Have you
> thoroughly exercised this change?

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//============================================================================/
/
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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;
}

OK, I've checked in the new function.  I did some tests and it looks good.  Thanks.

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.

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

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).