470,810 Members | 1,611 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 470,810 developers. It's quick & easy.

Matrix Inversion C code problem

Expand|Select|Wrap|Line Numbers
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <malloc.h>
  4. #include <stdlib.h>
  5. #include <time.h>
  6.  
  7.  
  8. struct SMatrix
  9. {
  10.     double** pValues;
  11.     int nR,nC;
  12. };
  13.  
  14. #define MAX_COLUMN 30
  15. #define MAX_ROW 30
  16.  
  17. void menu();
  18. int input_operation();
  19. void evaluate(int operation, float matrix[][MAX_ROW]);
  20. void addition(float matrix[][MAX_ROW]);
  21. void scalar_multiplication(float matrix[][MAX_ROW]);
  22. void matrix_multiplication(float matrix[][MAX_ROW]);
  23. void transpose(float matrix[][MAX_ROW]);
  24. struct SMatrix* inverse(float matrix[][MAX_ROW]);
  25. struct SMatrix* invert(float matrix[][MAX_ROW], int n);
  26. struct SMatrix* echelon(float m[(2*MAX_COLUMN)][MAX_ROW], int n);
  27. int find_max_row(float m[(2*MAX_COLUMN)][MAX_ROW], int cur_column, int start_row, int end);
  28. void reduce_echelon(float m[(2*MAX_COLUMN)][MAX_ROW], int index, int n);
  29. void reverse_echelon(float m[][MAX_ROW], int index, int n);
  30. void canonical_form(float m[][MAX_ROW], int index, int n);
  31. float get_product(float a[][MAX_ROW], float b[][MAX_ROW], int cur_col, int cur_row, int p);
  32. int get_size();
  33. float get_values();
  34. int input_scalar();
  35. void print_matrix(float m[][MAX_ROW], int column, int row);
  36.  
  37.  
  38. /*
  39.  *    This function is used to obtain the inverse of a matrix nxn. If the matrix input is not a
  40.  * square matrix, it prints out a message indicating that matrix is not invertible.
  41.  *    Arguments:
  42.  *        matrix:    an array where the values of the matrix are stored
  43.  *        size: an array where the size of the matrix are temporarily stored when obtained from the user
  44.  *    column, row: the size of the matrix
  45.  *        i, cur_column, cur_row:    indices, counters
  46.  *
  47.  */
  48.  
  49. struct SMatrix* inverse(struct SMatrix* A)
  50. {
  51.     int i,j;
  52.     float matrix[30][MAX_ROW];
  53.     for (i=0;i<3;i++)
  54.     {
  55.         for (j=0;j<3;j++)
  56.         {
  57.             matrix[i][j]=A->pValues[i][j];
  58.         }
  59.     }
  60.     int column=A->nR;
  61.     return invert(matrix, column);    
  62. }
  63.  
  64. /*
  65.  *    This function is a subfunction of inverse. It creates an augmented matrix M[A|I] used to
  66.  * work with Gaussian Elimination to get the inverse of a matrix where A is the matrix and 
  67.  * I is the identity matrix.
  68.  * Arguments:
  69.  *        matrix:     an array where the values of the matrix input are stored
  70.  *        m:            an array where the values of the augmented matrix are stored
  71.  *        i:            the identity matrix
  72.  *
  73.  */
  74.  
  75. struct SMatrix* invert(float matrix[][MAX_ROW], int n)
  76. {
  77.     float m[(2*MAX_COLUMN)][MAX_ROW];        /*augmented matrix*/
  78.     float i[MAX_COLUMN][MAX_ROW];                /*identity*/
  79.     int cur_row;
  80.     int cur_column;
  81.  
  82.     /*identity*/
  83.     for (cur_row = 0; cur_row < n; cur_row++)
  84.     {
  85.         for (cur_column = 0; cur_column < n; cur_column++)
  86.         {
  87.             if (cur_row == cur_column)
  88.             {
  89.                 i[cur_column][cur_row] = 1;
  90.             }
  91.             else
  92.             {
  93.                 i[cur_column][cur_row] = 0;
  94.             }
  95.         }
  96.     }
  97.  
  98.     /*augmented matrix, 1st half*/
  99.     for (cur_row = 0; cur_row < n; cur_row++)
  100.     {
  101.         for (cur_column = 0; cur_column < n; cur_column++)
  102.         {
  103.             m[cur_column][cur_row] = matrix[cur_column][cur_row];
  104.         }
  105.     }
  106.     /*augmented matrix, 2nd half*/
  107.     for (cur_row = 0; cur_row < n; cur_row++)
  108.     {
  109.         for (cur_column = n; cur_column < (2 * n); cur_column++)
  110.         {
  111.             m[cur_column][cur_row] = i[(cur_column - n)][cur_row];
  112.         }
  113.     }
  114.  
  115.     return echelon(m, n);
  116. }
  117.  
  118. /*
  119.  *    This is a subfunction of inverse where the augmented matrix is being reduced to echelon
  120.  * form. If through the process of row operations a row zeros is obtained at the left side of
  121.  * the augmented matrix or if any of the diagonal values is zero, it prints out an error 
  122.  * message that the matrix has no inverse.
  123.  *    Arguments:    
  124.  *        m: the augmented matrix
  125.  *        n: the size of the matrix being inverted
  126.  *        inverse: an array to store the values of the inverted matrix
  127.  *        temp: a temporary storage
  128.  */
  129.  
  130. struct SMatrix* echelon(float m[(2*MAX_COLUMN)][MAX_ROW], int n)
  131. {
  132.     int i;                                        /*counter*/
  133.     int cur_row;
  134.     int cur_column;
  135.     float inverse[MAX_COLUMN][MAX_ROW];
  136.    float temp;  
  137.  
  138.  
  139.     /*pivoting and swapping of rows*/
  140.     for (i = 0; i < (n-1); i++)
  141.     {
  142.         for (cur_row = i; cur_row < n; cur_row++)
  143.         {
  144.             int row_index_max = find_max_row(m, i, cur_row, n);
  145.  
  146.             for (cur_column = 0; cur_column < (2*n); cur_column++)
  147.             {
  148.                 temp = m[cur_column][cur_row];
  149.                 m[cur_column][cur_row] = m[cur_column][row_index_max];
  150.                 m[cur_column][row_index_max] = temp;
  151.             }
  152.  
  153.         }
  154.         reduce_echelon(m, i, n);
  155.  
  156.         if (m[i][i] == 0)
  157.         {
  158.             printf("\nMatrix has no inverse!\n\n");
  159.         }
  160.  
  161.     }
  162.  
  163.     /*reverse row operations*/
  164.     for (i = (n-1); i > 0;i--)
  165.     {
  166.  
  167.         reverse_echelon(m, i, n);
  168.  
  169.     }
  170.  
  171.     /*reduces m to canonical form*/
  172.     for (i = 0; i < n; i++)
  173.     {
  174.         canonical_form(m, i, n);
  175.     }
  176.  
  177.     /*gets the inverse*/
  178.     for (cur_row = 0; cur_row < n; cur_row++)
  179.     {
  180.         for (cur_column = 0; cur_column < n; cur_column++)
  181.         {
  182.             inverse[cur_column][cur_row] = m[(cur_column + n)][cur_row];
  183.         }
  184.     }
  185.  
  186.     int r,t;
  187.     struct SMatrix* Y;
  188.     Y->nC=3;
  189.     Y->nR=3;
  190.     for (r=0;r<3;r++)
  191.     {
  192.         for (t=0;t<3;t++)
  193.         {
  194.             Y->pValues[r][t]=inverse[r][t];
  195.         }
  196.     }
  197.         return (Y);
  198.  
  199. }
  200.  
  201. /*
  202.  * This function is used to obtain the row which has the leading non-zero entry, an important
  203.  * parameter in doing elementary row operations in obtaining the echelon form of an augmented
  204.  * matrix.
  205.  *    Arguments:
  206.  *     m: the augmented matrix
  207.  *        cur_column, start_row, end: indices and counters
  208.  *     index_max_row: the index of the row which has the leading non-zero entry
  209.  * Returns a value.
  210.  *
  211.  */
  212.  
  213. int find_max_row(float m[(2*MAX_COLUMN)][MAX_ROW], int cur_column, int start_row, int end)
  214. {
  215.     int index_max_row = start_row;
  216.     int cur_row;
  217.  
  218.     for (cur_row = start_row; cur_row < end; cur_row++)
  219.     {
  220.         if (m[cur_column][cur_row] > m[cur_column][index_max_row])
  221.         {
  222.             index_max_row = cur_row;
  223.         }
  224.     }
  225.     return index_max_row;
  226. }
  227.  
  228. /*
  229.  *    This function is used to get the echelon form of the augmented matrix.
  230.  * Arguments:
  231.  *        m: the augmented matrix
  232.  *        index: the current index of the values of the augmented matrix being reduced to echelon
  233.  *        n: the size of the matrix nxn being inverted
  234.  *        temp: temporary storage
  235.  *        factor: the factor used to reduce the values of the matrix into echelon
  236.  *
  237.  */
  238.  
  239. void reduce_echelon(float m[(2*MAX_COLUMN)][MAX_ROW], int index, int n)
  240. {
  241.     int cur_row;
  242.     int cur_column;
  243.     float temp;
  244.     float factor;
  245.  
  246.     for (cur_row = index; cur_row < (n-1); cur_row++)
  247.     {
  248.         factor = m[index][(cur_row + 1)]/m[index][index];
  249.         for (cur_column = index; cur_column < (n * 2); cur_column++)
  250.         {
  251.             temp = m[cur_column][(cur_row + 1)] - (factor * m[cur_column][index]);
  252.             m[cur_column][(cur_row + 1)] = temp;
  253.         }
  254.     }
  255. }
  256.  
  257. /*
  258.  * This function is used for the reverse elementary row operations.
  259.  *
  260.  */
  261.  
  262. void reverse_echelon(float m[(2*MAX_COLUMN)][MAX_ROW], int index, int n)
  263. {
  264.     int cur_row;
  265.     int cur_column;
  266.     float temp;
  267.     float factor;
  268.  
  269.     for (cur_row = index; cur_row > 0; cur_row--)
  270.     {
  271.         factor = m[index][(cur_row - 1)]/m[index][index];
  272.         for (cur_column = index; cur_column < (2 * n); cur_column++)
  273.         {
  274.             temp = m[cur_column][(cur_row - 1)] - (factor * m[cur_column][index]);
  275.             m[cur_column][(cur_row -1)] = temp;
  276.         }
  277.     }    
  278. }
  279.  
  280. /*
  281.  * This function is used to obtain the canonical form of the reduced echolon form of the 
  282.  * augmented matrix.
  283.  *
  284.  */
  285.  
  286. void canonical_form(float m[(2*MAX_COLUMN)][MAX_ROW], int index, int n)
  287. {
  288.     int cur_column;
  289.     float temp;
  290.     float factor;
  291.  
  292.     factor = m[index][index];
  293.  
  294.     for (cur_column = index; cur_column < (n*2); cur_column++)
  295.     {
  296.         temp = m[cur_column][index]/factor; 
  297.         m[cur_column][index] = temp;
  298.     }
  299. }
  300.  
  301. /*
  302.  * This function is used to obtain the value of ijth entry of the product matrix.
  303.  *    Arguments:
  304.  *        a, b: arrays where the values of the matrices A and B are stored
  305.  *        cur_col: the index of the current column value of the product
  306.  *        p: equivalent to the row of A and column of B
  307.  *     sum: the value of the ijth entry of the product of matrices A and B
  308.  *        product: temporary storage of the product of the ith row of A and jth column of B 
  309.  * Returns a value.
  310.  *
  311.  */
  312.  
  313. float get_product(float a[][MAX_ROW], float b[][MAX_ROW], int cur_col, int cur_row, int p)
  314. {
  315.     float sum = 0;
  316.     float product;
  317.     int i;
  318.  
  319.     for (i = 0; i < p; i++)
  320.     {
  321.         product = a[cur_col][i] * b[i][cur_row];
  322.         sum = sum + product;
  323.     }
  324.     return sum;
  325. }
  326.  
  327. /*
  328.  * This function is used to get the size of a matrix.
  329.  * Arguments:
  330.  *     n:        a variable where to store the size input by the user
  331.  * Returns a value.
  332.  *
  333.  */
  334.  
  335. int get_size()
  336. {
  337.     int n;
  338.  
  339.     scanf("%d", &n);
  340.     return n;
  341. }
  342.  
  343. /*
  344.  * This function is used to get the values of a matrix. It scans values in fractional form 
  345.  * (+|-)value_1/value_2 or in decimal form (+|-)value_1.value_2 where:
  346.  *            +,- is the sign
  347.  *            value_1 is the numerator (fractional form) and whole number (decimal form)
  348.  *       value_2 is the denominator (fractional form) and decimal value up to 2 digits (decimal form) 
  349.  * Arguments:
  350.  *        matrix:        an array to store the values of the matrix
  351.  *        cur_column, cur_row: the indices where to store the value in the array
  352.  *
  353.  */
  354.  
  355. float get_values()
  356. {
  357.  
  358.     int n;
  359.  
  360.     scanf("%d", &n);
  361.     return n;
  362. /*    
  363.     if (dummy == SLASH)
  364.     {
  365.         value = (float)(value_1/value_2);
  366.     }    
  367.     else if (dummy == POINT)
  368.     {
  369.         value = (float)((value_1 * 100)+ value_2)/ 100;
  370.     }
  371.  
  372.     if (sign == POS)
  373.     {
  374.         value = value;
  375.     }
  376.     else if (sign == NEG)
  377.     {
  378.         value = -(value);
  379.     }
  380. */    
  381. //    matrix[cur_column][cur_row] = value;
  382. }
  383.  
  384. /*
  385.  * This function displays the resulting matrix derived from the operations.
  386.  * Arguments:
  387.  *        msg:        the message from which the matrix was derived
  388.  *        matrix:    an array where the values of the matrix are stored
  389.  *        column, row: the size of the matrix
  390.  *
  391.  */
  392.  
  393. void print_matrix(float matrix[][MAX_ROW], int column, int row)
  394. {
  395.     int cur_row;
  396.     int cur_column;
  397.  
  398.     printf("Size: %d %d\n", column, row);
  399.     for (cur_row = 0; cur_row < row; cur_row++)
  400.     {
  401.         for (cur_column = 0; cur_column < column; cur_column++)
  402.         {
  403.             if (matrix[cur_column][cur_row] >= 0)
  404.             {
  405.                 printf("+%f ", matrix[cur_column][cur_row]);
  406.             }
  407.             else
  408.             {
  409.                 printf("%f ", matrix[cur_column][cur_row]);
  410.             }
  411.         }
  412.         printf("\n");
  413.     }
  414. }
  415.  
  416. void PrintMatrix(struct SMatrix* pM) 
  417. {
  418.     int r,c;
  419.     for (r=0;r<pM->nR;r++)
  420.     {
  421.         for(c=0;c<pM->nC;c++)
  422.         {
  423.             printf("%5.5f\t\t", pM->pValues[r][c]);
  424.         }
  425.         printf("\n");
  426.     }
  427. }
  428.  
  429. /*
  430.  * This function obtains the scalar constant where the value can be any integer (+/-).
  431.  * Arguments:
  432.  *     value:    where the integer is stored
  433.  *        sign:        where the sign of the integer is stored
  434.  *     POS:        constant character equivalent to (+) sign
  435.  *        NEG:        constant character equivalent to (-) sign
  436.  * Returns a value.
  437.  *
  438.  */
  439.  
  440.  
  441.  
  442. enum MFILLTYPE{MFT_NONE,  MFT_ZEORS, MFT_ONES,  MFT_RAND,  MFT_JEDENTITY};
  443. struct SMatrix* CreateMatrix(int nR,int nC,enum MFILLTYPE fillType)
  444. {
  445.     int r,c;
  446.     struct SMatrix* pRM;
  447.     pRM =(struct SMatrix*)malloc(sizeof(struct SMatrix));
  448.     pRM->nR = nR;
  449.     pRM->nC = nC;
  450.     //Allocate the matrix,  the matrix is an array of pointer to rows pointers
  451.     pRM->pValues = (double**)malloc(nR * sizeof(double*));
  452.     //Allocate the rows,  each row will point to all row elements for(r=0;r<nR;r++)
  453.     for (r=0;r<nR;r++)
  454.         pRM->pValues[r] = (double*)malloc(nC * sizeof(double));
  455.     if(fillType!=MFT_NONE)
  456.     {
  457.         for(r=0;r<nR;r++)
  458.             for(c=0;c<nC;c++)
  459.             {
  460.                 switch(fillType) 
  461.                 {
  462.                     case MFT_ZEORS:
  463.                         pRM->pValues[r][c] = 0; break; 
  464.                     case MFT_ONES:
  465.                         pRM->pValues[r][c] = 1; break;
  466.                     case MFT_RAND:
  467.                         pRM->pValues[r][c] = (rand()%100)/100.0; break;
  468.  
  469.  
  470.                 }
  471.             }
  472.     }
  473.     return pRM;
  474. }
  475.  
  476. void main ()
  477. {
  478.     struct SMatrix *A,*B;
  479.     A=CreateMatrix(3,3,MFT_RAND);
  480.     B=inverse(A);
  481.     PrintMatrix(B);
  482.  
  483. }
Jul 4 '10 #1
6 4535
weaknessforcats
9,208 Expert Mod 8TB
What is your question exactly?

It helps to state your problem and show only the code that applies to it rather than posting the whole thing and leaving it to us to figure out what might be your question.
Jul 4 '10 #2
line 197, it refuses to return struct SMatrix Y
Jul 4 '10 #3
here

Jul 4 '10 #4
problem solved
i forgot to use the CreateMatrix function :D
Jul 4 '10 #5
weaknessforcats
9,208 Expert Mod 8TB
This code in the echelon function:

Expand|Select|Wrap|Line Numbers
  1. struct SMatrix* Y; 
  2.     Y->nC=3; 
  3.     Y->nR=3; 
shows Y being defined. The address in Y is junk.

Then junk is used to access nC. Result: crash. My compiler won't compile the code.

You need a valid address in Y.

BTW: There are mixtures of float and double in the code. A conversion of double to float doesn't work without possible loss of data. I suggest you use double throughout. Personlly, I would not use float unless there was a technical reason for me to do so that I could write down on paper.
Jul 4 '10 #6
thnx i'll change them :)
Jul 4 '10 #7

Post your reply

Sign in to post your reply or Sign up for a free account.

Similar topics

2 posts views Thread by Jan-Hendrik Huehne | last post: by
15 posts views Thread by christopher diggins | last post: by
4 posts views Thread by Bradley Burton | last post: by
19 posts views Thread by snowdy | last post: by
3 posts views Thread by lancered | last post: by
6 posts views Thread by sheff | last post: by
reply views Thread by mihailmihai484 | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.