First Blast from the Past (1997) Computing the Inverse Image of a Matrix by James Pate Williams, Jr.

/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.3.5 (Inverse Image Matrix). Let M be
  an m by n matrix and V be an m by r matrix, where
  n <= m. This algorithm either outputs a message
  saying that some column vector of V is not in the
  image of M, or outputs an n by r matrix X such
  that V = MX." -Henri Cohen- See "A Course in Com-
  putational Algebraic Number Theory" by Henri
  Cohen pages 60-61. We specialize to the field Q.
*/

#include <stdio.h>
#include <stdlib.h>

double** create_matrix(long m, long n)
{
    double** matrix = (double**)calloc(m, sizeof(double*));
    long i;

    if (!matrix) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from create_matrix\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        matrix[i] = (double*)calloc(n, sizeof(double));
        if (!matrix[i]) {
            fprintf(stderr, "fatal error\ninsufficient memory\n");
            fprintf(stderr, "from create_matrix\n");
            exit(1);
        }
    }
    return matrix;
}

void delete_matrix(long m, double** matrix)
{
    long i;

    for (i = 0; i < m; i++) free(matrix[i]);
    free(matrix);
}

void inverse_image_matrix(long m, long n, long r,
    double** M, double** V,
    double** X)
{
    double ck, d, sum, t;
    double** B = create_matrix(m, r);
    int found;
    long i, j, k, l;

    for (i = 0; i < m; i++)
        for (j = 0; j < r; j++)
            B[i][j] = V[i][j];
    for (j = 0; j < n; j++) {
        found = 0, i = j;
        while (!found && i < m) {
            found = M[i][j] != 0;
            if (!found) i++;
        }
        if (!found) {
            fprintf(stderr, "fatal error\nnot linearly independent\n");
            fprintf(stderr, "from inverse_image_matrix\n");
            exit(1);
        }
        if (i > j) {
            for (l = 0; l < n; l++)
                t = M[i][l], M[i][l] = M[j][l], M[j][l] = t;
            for (l = 0; l < r; l++)
                t = B[i][l], B[i][l] = B[j][l], B[j][l] = t;
        }
        d = 1.0 / M[j][j];
        for (k = j + 1; k < m; k++) {
            ck = d * M[k][j];
            for (l = j + 1; l < n; l++)
                M[k][l] -= ck * M[j][l];
            for (l = 0; l < r; l++)
                B[k][l] -= ck * B[j][l];
        }
    }
    for (i = n - 1; i >= 0; i--) {
        for (k = 0; k < r; k++) {
            sum = 0;
            for (j = i + 1; j < n; j++)
                sum += M[i][j] * X[j][k];
            X[i][k] = (B[i][k] - sum) / M[i][i];
        }
    }
    for (k = n + 1; k < m; k++) {
        for (j = 0; j < r; j++) {
            sum = 0;
            for (i = 0; i < n; i++)
                sum += M[k][i] * X[i][j];
            if (sum != B[k][j]) {
                fprintf(stderr, "fatal error\ncolumn not in image\n");
                fprintf(stderr, "from inverse_image_matrix\n");
                exit(1);
            }
        }
    }
    delete_matrix(m, B);
}

void matrix_multiply(long m, long n, long r,
    double** a, double** b,
    double** c)
    /* c = a * b */
{
    double sum;
    long i, j, k;

    for (i = 0; i < m; i++) {
        for (j = 0; j < r; j++) {
            sum = 0.0;
            for (k = 0; k < n; k++)
                sum += a[i][k] * b[k][j];
            c[i][j] = sum;
        }
    }
}

void print_matrix(long m, long n, double** a)
{
    long i, j;

    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%+10.6lf ", a[i][j]);
        printf("\n");
    }
}

int main(void)
{
    long i, j, m = 4, n = 4, r = 4;
    double** c = create_matrix(m, n);
    double** M = create_matrix(m, n);
    double** V = create_matrix(m, r);
    double** X = create_matrix(n, r);

    for (i = 0; i < m; i++) {
        c[i][i] = M[i][i] = 2.0;
        if (i > 0)
            c[i][i - 1] = M[i][i - 1] = -1.0;
        if (i < m - 1)
            c[i][i + 1] = M[i][i + 1] = -1.0;
    }
    for (i = 0; i < m; i++)
        for (j = 0; j < r; j++)
            V[i][j] = i + j + 1;
    printf("M is as follows:\n");
    print_matrix(m, n, M);
    printf("V is as follows:\n");
    print_matrix(m, r, V);
    inverse_image_matrix(m, n, r, M, V, X);
    printf("X is as follows:\n");
    print_matrix(n, r, X);
    matrix_multiply(m, n, r, c, X, M);
    printf("MX is as follows:\n");
    print_matrix(m, r, M);
    delete_matrix(m, c);
    delete_matrix(m, M);
    delete_matrix(m, V);
    delete_matrix(n, X);
    return 0;
}