Thread: Trying to make this code faster & Cramer

1. Trying to make this code faster & Cramer's Rule

I'm trying to fit a parabola to a segment of data using the two end points and a point near the center.

I've decided to use Cramer's rule because I think that it's much
easier to only worry with determinants instead of having to find
the inverse of the matrix. I was wondering if there is a better or

Also, since my determinants are only for 3x3 matrices would it be
a better idea to simply crank out the equation:
Code:
```double m00 = array[0][0];
double m01 = array[0][1];
double m02 = array[0][2];
double m10 = array[1][0];
double m11 = array[1][1];
double m12 = array[1][2];
double m20 = array[2][0];
double m21 = array[2][1];
double m22 = array[2][2];

double determinant = m00*m11*m22 + m01*m12*m20 + m02*m10*m21 -
m02*m11*m20 - m01*m10*m22 - m00*m12*m21;```
I realize that readability is in favor of the matrix_det function, but I'm concerned with the
time factor. This code needs to be fairly efficient. This particular loop of code is
executed anywhere from 100 to 1000 times per minute.
My concerns stem from stepping through the code and my inexperience with pointer
objects and memory control. I understand how the array of pointers is the matrix
structure that is passed to the function, which is a similar version of Prelude's determinant
function.
When I step through the code it jumps through all kinds of strange code in malloc.c,
dbgheap.c, new.cpp, memset.asm . . . would it be any less time consuming to just create
an array and carefully create a mapping of a 2d array to a 1d array?

Ex: 2D[0][0]->1D[0] . . . 2D[2][2]->1D[9].
(Or is this concern pointless) I think that the readability will be even more convoluted,
but from a speed standpoint is this really any faster?

Here's the current version of the code:

Code:
``` //parameters[3] is set at 10 for now
if ( consecutiveFlags[c] > parameters[3] )
{

//*******Set up for solving the system of equations
double **ARRAY, **ARRAY2;
ARRAY  = new double *[Q];	//Q defines the square matrix dimensions
ARRAY2 = new double *[Q];

double YDATA[3];

for ( i = 0; i < Q; i++ )
{
ARRAY[i]  = new double[Q];
ARRAY2[i] = new double[Q];
}//end set up storage

//set up for first row of matrix		//grab pt @ beginning of segment
ARRAY[0][0] = zlen[c]*zlen[c];
ARRAY[0][1] = zlen[c];
ARRAY[0][2] = 1.0;

//set up for second row of matrix	//grab pt near middle of
k = c + consecutiveFlags[c]/2-1;		//failed segment
ARRAY[1][0] = zlen[k]*zlen[k];
ARRAY[1][1] = zlen[k];
ARRAY[1][2] = 1.0;

//set up for third row of matrix		//grab pt @ end of failed segment
k = c + consecutiveFlags[c]-1;
ARRAY[2][0] = zlen[k]*zlen[k];
ARRAY[2][1] = zlen[k];
ARRAY[2][2] = 1.0;

/*To find the coefficients of the parabola for the segment correction we must solve:
[zlen1^2  zlen1  1][A]  [y_front ]
[zlen2^2  zlen2  1][B]= [y_middle]   or   ARRAY * Coefficients = YDATA
[zlen3^2  zlen3  1][C]  [y_end   ]

(Use Cramer's Rule)
Step 1)Check for det of ARRAY, if zero break out of loop b/c division by zero
Step 2)Set up ARRAY2 with the original ARRAY matrix data but substituting YDATA
transpose for each column corresponding to the position of the desired
Coefficient and then calculating the determinant.
Step 3)Coefficients will be equal to to ratio of the determinant of ARRAY2 to the
determinant of ARRAY.
*/
//Step 1:
double detARRAY = matrix_det ( ARRAY, 3 );
if ( detARRAY == 0 )
break;				//error
else
;//nothing			//continue

//Step 2:
double Coefficients[3];

for ( int m = 0; m < Q; m++ )				//m iterates the Coefficients
{	for ( j = 0; j < Q; j++ )				//j iterates the columns
{	for( i = 0; i < Q; i++ )			//i iterates the rows
{	if ( j==m )
ARRAY2[i][j] = YDATA[i];	//fill appropriate column with ydata
else
ARRAY2[i][j] = ARRAY[i][j];//copy original array in other columns
}
}//end for populate ARRAY2

double detARRAY2 = matrix_det ( ARRAY2, 3 );
Coefficients[m] = detARRAY2/detARRAY;
}//end for populate Coefficients

//created an array object to model a matrix
//must delete it to avoid a memory leak
for ( i = 0; i < Q; i++ )
{
delete [] ARRAY[i];
delete [] ARRAY2[i];
}//delete storage containers

delete [] ARRAY;
delete [] ARRAY2;

}//end if find the parabola for the segment

double matrix_det ( double **in_matrix, int n )
{
int i, j, k;
double **matrix;
double det = 1;

matrix = new double *[n];

for ( i = 0; i < n; i++ )				//setting up pointers for rows i, [0,n-1]
matrix[i] = new double[n];

for ( i = 0; i < n; i++ )
{
for ( j = 0; j < n; j++ )
matrix[i][j] = in_matrix[i][j];	//setting up pointers for columns j, [0,n-1]
}//end pointer setup

for ( k = 0; k < n; k++ )
{
if ( matrix[k][k] == 0 ) 			//check for zeros along the diagonal
{
bool ok = false;

for ( j = k; j < n; j++ ) 		//ensures that if a zero exists along the diagonal
{									//at least one value in that column has to be non-zero
if ( matrix[j][k] != 0 )
ok = true;
}//end for check if det will be zero

if ( !ok )						//if column of zeros exist, then the matrix has det=0
return 0;

for ( i = k; i < n; i++ )
{
double temp;
temp = matrix[i][k];
matrix[i][k] = matrix[i][j];
matrix[i][j] = temp;
//std::swap ( matrix[i][j], matrix[i][k] );
}//end for swap column k for column j, the last column

det = -det;

}//end if a zero occurs along diagonal

det *= matrix[k][k];

if ( k + 1 < n )
{
for ( i = k + 1; i < n; i++ )
{
for ( j = k + 1; j < n; j++ )
{
matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j] / matrix[k][k];
}//end for
}//end for
}//end if

}//end for k's

for ( i = 0; i < n; i++ )
delete [] matrix[i];

delete [] matrix;

return det;```

2. Can you edit your post for width?
In particular, just put actual code in code tags - long sentences in code tags are bad

3. Sorry

I was in the middle of posting and it jumped to submit by accident

4. Originally Posted by just2peachy
I'm trying to fit a parabola to a segment of data using the two end points and a point near the center.

I've decided to use Cramer's rule because I think that it's much
easier to only worry with determinants instead of having to find
the inverse of the matrix. I was wondering if there is a better or

Also, since my determinants are only for 3x3 matrices would it be
a better idea to simply crank out the equation:
A couple of things come to mind:

1. If your program is dealing only with 3x3 matrices/vectors, you can optimize your functions by using fixed size arrays (no new or delete required).

2. For debug purposes, just use the direct calculation that you have shown for the determinants. You can tell at a glance that the formula is correct. You can put the determinant calculation in a function. (The function will simply perform the operations on the matrix elements; no copying of matrices etc. required.) You could consider specific optimizations for the three of the four determinant calculations that have a column of all '1'. You could consider putting the calculations in-line rather than in functions, etc., etc.

3. If performance is still a problem you should realize that Cramer's rule (requires four determinant calculations plus a couple more multiplication/add operations) is not the most efficient way to solve systems of three equations. Look into methods with the name "Gauss" in them (maybe throw in a Cholesky or Crout also). Properly implemented, Gauss reduction with row (and, maybe, column) pivoting can also minimize roundoff error, and still require fewer operations than direct calculation of four determinants.

Regards,

Dave

Popular pages Recent additions