[Solved] Solving linear equations via adjunt matrix in C [closed]

Question

The problem with getting uncommented code “from a senior (buddy) of mine“, is that you have absolutely no guidance and zero reference for the algorithms used in the code. Many times this leaves you in a worse position than you would be in if you had taken the time to learn the math involved and coded the solution yourself.

Why? Instead of getting the benefit of learning the math involved and then logically applying that knowledge in forming a numeric solution, you end up guessing at where the error is in the code and needlessly guessing, recompiling, and failing again.

Compounding that fact, the code you got from your buddy doesn’t even contain valid C. There is no “For” loop in C, there is no “If” statement and there is no “Printf” function. In fact none of the standard library functions and none of the operators start with a Capital Letter at all. The inclusion of the curses.h header is superfluous.

Further, you invoke Undefined Behavior by allowing the value of j to exceed the bounds of your array A in the loop:

for(i=0;i<3;i++){ 
   swap (&A[i][j], &A[j][i]);}

Why? Because the value of j is 3 following the nested loops where you fill A from cof_x and it is never reset prior to attempting the transpose.

(note: that turns out to be only one of the problems with your attempt to transpose the cofactor matrix)

One of you primary complaints is the code does not compile. Of course it doesn’t. However, if you enable compiler warnings the compiler would tell you the exact line where the problem occurs helping you fix the code so it does compile. For gcc, add -Wall -Wextra -pedantic to your compile string. For clang, add -Weverything and for VS (cl.exe) add /W3 (or /Wall for literally all warnings). Do not accept code until it compiles cleanly without warning.

Once you handle the basics, you can turn to the specifics of the algorithms. While the code forming the adjugate (or adjunct of squares) looked suspect, it in fact properly forms for adjugate matrix. The determinate is properly calculated. Where your code fails is computing the transpose the adjugate matrix by improperly swapping elements while iterating over the entire matrix.

You have two options to transpose. Either iterate over the entire bounds assigning A[j][i] = A[i][j]; or for a square matrix, iterate over only those elements that change swapping the values.

After correcting the transpose (and commenting the code with reference to the operation at issue), the code works fine, e.g.

#include <stdio.h>

#define NDIM 3  /* if you need a constant - define one (or more) */

void swap (float *a, float *b)
{
    float temp;
    temp = *a;
    *a = *b;
    *b = temp;
}

int main (void) {

    int i, j, cof_1, cof_2, cof_3, cof_4;
    float   a[NDIM][NDIM]   = {{0}},    /* coefficient matrix */
            c[NDIM]         =  {0},     /* solution vector */
            A[NDIM][NDIM]   = {{0}},    /* adjunct matrix */
            INV[NDIM][NDIM] = {{0}},    /* inverse of a */
            det             =  0.0,     /* determinant of a */
            X[NDIM]         =  {0},     /* linear system roots */
            B[NDIM]         =  {0};     /* product A * c */

    i = j = cof_1 = cof_2 = cof_3 = cof_4 = 0;  /* initialize */

    for (i = 0; i < NDIM; i++) {    /* validate input */
        printf ("Eq[%d] coefficients and constant: ", i + 1);
        if (scanf ("%f%f%f%f", &a[i][0], &a[i][1], &a[i][2], &c[i]) != 4) {
            fprintf (stderr, "error: invalid coefficient.\n");
            return 1;
        }
    }
    putchar ('\n');

    /* cofactor matrix from matrix of minors incorporating
     * patterned +/- sign to maxtix of minors.
     */
    for (i = 0; i < NDIM; i++)
        for (j = 0; j < NDIM; j++) {
            cof_1 = (i + 1) % NDIM;
            cof_2 = (j + 1) % NDIM;
            cof_3 = (i + 2) % NDIM;
            cof_4 = (j + 2) % NDIM;

            A[i][j]   = (a[cof_1][cof_2] * a[cof_3][cof_4]) -
                        (a[cof_1][cof_4] * a[cof_3][cof_2]);
        }

    /* compute determinant */
    for (i = 0; i < NDIM; i++)
        det += a[i][0] * A[i][0];

    /* transpose cofactor (adjunct) */
    for (i = 1; i < NDIM; i++)
        for (j = 0; j < i; j++)
            swap (&A[i][j], &A[j][i]);

    /* consistency check B */
    for (i = 0; i < NDIM; i++)
        for (j = 0; j < NDIM; j++)
            B[i] += A[i][j] * c[j];

    if (det == 0) { /* eliminate infinite and no solution cases */
        if (B[0] == 0 && B[1] == 0 && B[2] == 0) {
            printf ("The system is consistent and there "
                    "are infinitely many solutions.\n");
            return 1;
        } else
            printf ("There is no solution to this system.\n");
        return 2;
    }

    /* compute inverse of A (1/det * A) */
    for (i = 0; i < NDIM; i++)
        for (j = 0; j < NDIM; j++)
            INV[i][j] = A[i][j] / det;

    /* multiply constant vector by inverse */
    for (i = 0; i < NDIM; i++)
        for (j = 0; j < NDIM; j++)
            X[i] += INV[i][j] * c[j];

    if (X[0] == 0 && X[1] == 0 && X[2] == 0)
        /* The case for trivial solution */
        printf ("This system has a trivial solution.\n");
    else
        /* The case for unique solution */
        printf ("This system has a unique solution.\n\n");

    /* output solutions for linear system of eq. */
    for (i = 0; i < NDIM; i++)
        printf ("\t%c = %.2f\n", 'x' + i, X[i]);

    return 0;
}

Example Use/Output

Linear equations from a simple proof of a 30-60-90 triangle used for the example:

$ ./bin/matrixinv
Eq[1] coefficients and constant: 1 1 1 180
Eq[2] coefficients and constant: 5 -1 -1 0
Eq[3] coefficients and constant: 1 1 -1 0

This system has a unique solution.

        x = 30.00
        y = 60.00
        z = 90.00

Look things over and let me know if you have further questions.

Accepted Answer