A Simple Simultaneous Equation Solver

Listing 1: The matrix inversion function

/* =================================================================== */
/* augvrt - inverts an nr by nc matrix a, and computes its determinant */
/* =================================================================== */
# if defined(__STDC__) || defined(__PROTO__)
double augvrt(double *a, int nr, int nc, int nd)
# else
double augvrt(a, nr, nc, nd)
double *a;
int     nr, nc, nd;
# endif
/* ------------------------------------------------------------ */
/*-     Explanation of Calling Sequence Parameters              */
/* ------------------------------------------------------------ */
/*-     nr - no. rows in a                                      */
/*-     nc - no. cols in a (nc >= nr)                           */
/*-     nd - order of matrix that holds matrix a (nd >= nc)     */
/*-          (matrix to be inverted can be imbedded)            */
/*-                                                             */
/*-     if x = a-inv*b is being solved, there will be nc-nr     */
/*-     solution vectors in x                                   */
/* ------------------------------------------------------------ */
{
    double  det_val = 1.0;              /* determinant value */

    register
    int     i, j, k;
    double  r;
    int     ii, ij, ir, jr, ki, kj;

    ii  = 0;                            /* index of diagonal */
    ir  = 0;                            /* index of row */

    for (i = 0; i < nr; ++i, ir += nd, ii = ir + i)
    {
        det_val *= a[ii];               /* new value of determinant */

        /* ------------------------------------ */
        /* if value is zero, might be underflow */
        /* ------------------------------------ */

        if (det_val == 0.0)
        {
            if (a[ii] == 0.0)
            {
                break;                  /* error - exit now */
            }
            else                        /* must be underflow */
            {
                det_val = 1.0e-30;
            }
        }

        /* --------------- */
        /* Calculate Pivot */
        /* --------------- */
        r = 1.0 / a[ii];

        a[ii] = 1.0;

        ij = ir;                        /* index of pivot row */

        for (j = 0; j < nc; ++j)
        {
            a[ij] = r * a[ij];
            ++ij;                       /* index of next row element */
        }

        ki = i;                         /* index of ith column */
        jr = 0;                         /* index of jth row */

        for (k = 0; k < nr; ++k, ki += nd, jr += nd)
        {
            if (i != k && a[ki] != 0.0)
            {
                r = a[ki];              /* pivot target */

                a[ki] = 0.0;

                ij = ir;                /* index of pivot row */
                kj = jr;                /* index of jth row */

                /* -------------------------------------------- */
                /* subtract multiples of pivot row from jth row */
                /* -------------------------------------------- */
                for (j = 0; j < nc; ++j)
                {
                    a[kj] -= r * a[ij];
                    ++ij, ++kj;
                }
            }
        }
    }

    return(det_val);
}                                       /* augvrt */