Thread: Trying to make this code faster & Cramer

  1. #1
    Registered User
    Join Date
    Oct 2004
    Posts
    17

    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
    faster way to go about this.

    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;						
    				YDATA[0] = yadjusted[c];
    
    				//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;
    				YDATA[1] = yadjusted[k];
    
    				//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;
    				YDATA[2] = yadjusted[k];
    
    	/*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;
    Last edited by just2peachy; 12-03-2004 at 08:12 AM.

  2. #2
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,660
    Can you edit your post for width?
    In particular, just put actual code in code tags - long sentences in code tags are bad

  3. #3
    Registered User
    Join Date
    Oct 2004
    Posts
    17

    Sorry

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

  4. #4
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote 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
    faster way to go about this.

    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
    Last edited by Dave Evans; 12-03-2004 at 04:12 PM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Some help with make my programs faster
    By Sshakey6791 in forum C++ Programming
    Replies: 11
    Last Post: 12-11-2008, 01:41 PM
  2. No clue how to make a code to solve problems!
    By ctnzn in forum C Programming
    Replies: 8
    Last Post: 10-16-2008, 02:59 AM
  3. Explain this C code in english
    By soadlink in forum C Programming
    Replies: 16
    Last Post: 08-31-2006, 12:48 AM
  4. Obfuscated Code Contest
    By Stack Overflow in forum Contests Board
    Replies: 51
    Last Post: 01-21-2005, 04:17 PM
  5. << !! Posting Code? Read this First !! >>
    By kermi3 in forum C Programming
    Replies: 0
    Last Post: 10-03-2002, 03:04 PM