/*
  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 = 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] = 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;
}