473,324 Members | 2,248 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,324 software developers and data experts.

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 4841
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

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

Similar topics

2
by: Jan-Hendrik Huehne | last post by:
hello, iv've got a little problem! i need a script to invert a 4x4 matrix! has anybody such a script or do you know a download source??
15
by: christopher diggins | last post by:
Here is some code I wrote for Matrix multiplication for arbitrary dimensionality known at compile-time. I am curious how practical it is. For instance, is it common to know the dimensionality of...
4
by: Bradley Burton | last post by:
I'm using Allen Brown's code for audit logging (http://allenbrowne.com/AppAudit.html), but I'm having a problem. My aud table doesn't populate with the tracking info at all. I think it might be a...
19
by: snowdy | last post by:
I am using Interactive C with my handboard (68HC11) development system but I've got a problem that I am asking for help with. I am not new to C but learning all the time and I just cant see how to...
2
by: Praveen K | last post by:
I have a problem in communicating between the C# and the Excel Interop objects. The problem is something as described below. I use Microsoft Office-XP PIA dll’s as these dll’s were been...
6
by: amit kumar singh | last post by:
Hi everyone I am new to this site and also new to programming world can anybody help me writing C code for matrix multiplication (without using pointers) Amit
3
by: lancered | last post by:
Hi dear all, I am using Python2.4.2+NumPy1.0.1 to deal with a parameter estimation problem with the least square methods. During the calculations, I use NumPy package to deal with matrix...
3
by: PeteJ | last post by:
Hello all.. I wrote code in C to invert n x n matrices, where 1<n<39 and the matrices are guaranteed to be invertible. Prior to this task, I haven't really seen linear algebra since my sophmore...
6
by: sheff | last post by:
HI everybody, I was trying to use the following code to invert a matrix, everytime I try to run this code, the only error msg that comes up is '' inv030710.c:83: error: expected initializer...
0
by: DolphinDB | last post by:
Tired of spending countless mintues downsampling your data? Look no further! In this article, you’ll learn how to efficiently downsample 6.48 billion high-frequency records to 61 million...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
1
by: PapaRatzi | last post by:
Hello, I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...
1
by: CloudSolutions | last post by:
Introduction: For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.