/* Vladimir Kotlyar Sepember 1995
 *
 * Matrix manipulation
 *
 * NOTE: To be consistent with various descriptions of matrix
 * algorithms, all indices start from ONE
 *
 */

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

#include "alpha_int.h"
#include "alpha_matrix.h"

FILE *fdbg = stdout;

xa_vec_t *
xa_alloc_vector (arena_info *arena, int n)
{
  xa_vec x;

  assert (n >= 0);
  x = a_alloc0(arena, xa_vec_t);
  xa_size_vec (x) = n;
  xa_elts_vec (x) = (n > 0) ? a_nalloc0(arena,int,n) : 0;

  return x;
}

AlphaVector
a_make_vector (arena_info *arena, int n)
{
  xa_vec_t *	v;
  int		i;

  assert (n >= 0);

  v = xa_alloc_vector (arena, n);
  assert (v);
  a_size_vec(v) = n;
  for (i = 1; i <= n; i++) {
    a_elt_vec(v,i) = 0;
  }

  return v;
}

void
a_set_vec (AlphaVector v, int x)
{
  int		i;

  for (i = 1; i <= a_size_vec(v); i++) {
    a_elt_vec(v,i) = x;
  }
}

void
a_neg_vec (AlphaVector x)
{
  int		i;

  for (i = 1; i <= a_size_vec (x); i++) {
    a_elt_vec (x, i) = - a_elt_vec (x, i);
  }
}

void
a_saxpy (int a, AlphaVector x, int b, AlphaVector y, AlphaVector z)
{
  int		i;

  assert (a_size_vec(x) == a_size_vec(y));
  assert (a_size_vec(x) == a_size_vec(z));

  for (i = 1; i <= a_size_vec (x); i++) {
    a_elt_vec (z,i) = a * a_elt_vec(x,i) + b * a_elt_vec(y,i);
  }
}

void		
a_resize_vec (arena_info *arena, AlphaVector *px, int newn)
{
  AlphaVector	x = *px;
  AlphaVector 	xnew;
  int		n = a_length_vec (x);
  int		i;
  int		minn = (newn < n) ? newn : n;

  assert (newn >= 0);

  xnew = xa_alloc_vector (arena, newn);
  for (i = 1; i <= minn; i++) {
    a_elt_vec (xnew, i) = a_elt_vec (x, i);

#if !NDEBUG
    a_elt_vec (x, i) = (int) 0xbadabada;
#endif

  }

#if ! NDEBUG
  for (i = minn+1; i <= n; i++) {
    a_elt_vec (x, i) = (int) 0xbadbbadb;
  }
#endif

  *px = xnew;
}

void
a_copy_vec (AlphaVector dest, AlphaVector src)
{
  int		i;

  assert (a_size_vec(dest) == a_size_vec(src));

  for (i = 1; i <= a_size_vec(src); i++) {
    a_elt_vec(dest,i) = a_elt_vec(src,i);
  }
}

int		
a_is_zero_vec (AlphaVector x)
{
  int i;
  int is_zero = 1;

  for (i = 1; i <= a_length_vec (x); i++) {
    is_zero = is_zero && (0 == a_elt_vec(x,i));
  }

  return is_zero;
}

void
a_swap_vec (AlphaVector *px, AlphaVector *py)
{
  AlphaVector tmp;
  tmp = *px;
  *px = *py;
  *py = tmp;
}

int
a_gcd_vec (AlphaVector x)
{
  int		g;
  int		i;

  g = a_abs(a_elt_vec(x,1));
  for (i = 2; i <= a_length_vec(x); i++) {
    g = a_gcd(g, a_abs(a_elt_vec(x,i)));
  }

  return g;
}

void
a_extract_vec (AlphaVector x, int i1, int i2, AlphaVector y)
{
  int		n = a_length_vec (x);
  int		n1 = i2 - i1 + 1;
  int		j, j2;

  assert (n >= 0);
  assert (n1 >= 0);
  assert (i1 >= 1);
  assert (i2 <= n);

  for (j = i1, j2 = 1; j <= i2; j++, j2++) {
    a_elt_vec (y, j2) = a_elt_vec (x, j);
  }
}

void
a_read_vector (FILE *fp, AlphaVector x)
{
  int		i;
  int		rc;

  for (i = 1; i <= a_size_vec (x); i++) {
    rc = fscanf (fp, "%d", & a_elt_vec(x,i));
    assert (1 == rc);
  }
}

void
a_dump_vector (FILE *fp, AlphaVector x)
{
  int		i;

  fprintf (fp, "[");
  for (i = 1; i <= a_size_vec(x); i++) {
    fprintf (fp, " %d", a_elt_vec(x,i));
  }
  fprintf (fp, " ]'");

#ifndef NDEBUG
  fflush (fp);
#endif
}


/**************************************************/

xa_mx_t *
xa_alloc_matrix (arena_info *arena, int m, int n)
{
  xa_mx A;

  assert (n >= 0);
  assert (m >= 0);
  
  A = a_alloc0 (arena, xa_mx_t);
  xa_col_data_mx(A) = (n > 0) ? a_nalloc0(arena, xa_vec, n) : 0;
  xa_rows_mx (A) = m;
  xa_cols_mx (A) = n;

  return A;
}

AlphaMatrix
a_make_matrix (arena_info *arena, int m, int n)
{
  int		j;
  xa_mx_t *	A;

  assert (m >= 0);
  assert (n >= 0);

  A = xa_alloc_matrix (arena, m, n);
  assert (A);
  a_rows_mx(A) = m;
  a_cols_mx(A) = n;

  for (j = 1; j <= a_cols_mx(A); j++) {
    a_col_mx(A,j) = a_make_vector (arena, a_rows_mx(A));
  }

  return A;
}

void
a_copy_mx (AlphaMatrix dst, AlphaMatrix src)
{
  int		n = a_cols_mx (dst);
  int		m = a_rows_mx (dst);
  int		j;

  assert (n == a_cols_mx(src));
  assert (m == a_rows_mx(src));

  for (j = 1; j <= n; j++) {
    a_copy_vec (a_col_mx (dst,j), a_col_mx (src, j));
  }
}

void
a_set2ident (AlphaMatrix A)
{
  int		i;
  int		n = a_cols_mx (A);
  assert (n == a_rows_mx (A));

  for (i = 1; i <= n; i++) {
    a_set_vec (a_col_mx(A,i), 0);
    a_elt_vec (a_col_mx(A,i), i) = 1;;
  }
}

void
a_swap_col (AlphaMatrix A, int i, int j)
{
  AlphaVector 	tmp;

  tmp = a_col_mx(A,i); /* just a pointer swap !!! */
  a_col_mx(A,i) = a_col_mx(A,j);
  a_col_mx(A,j) = tmp;
}

void		
a_swap_row (AlphaMatrix A, int i, int j)
{
  int 		m = a_rows_mx (A);
  int 		n = a_cols_mx (A);
  int		k;
  int		tmp;

  assert (1 <= i && i <= m);
  assert (1 <= j && j <= m);

  for (k = 1; k <= n; k++) {
    tmp = a_elt_mx(A,i,k);
    a_elt_mx(A,i,k) = a_elt_mx(A,j,k);
    a_elt_mx(A,j,k) = tmp;
  }
}

void
a_matmult (arena_info *arena, AlphaMatrix A, AlphaMatrix B, AlphaMatrix C)
{
  int 		j;

  assert (a_rows_mx(A) == a_rows_mx(C));
  assert (a_cols_mx(A) == a_rows_mx(B));
  assert (a_cols_mx(B) == a_cols_mx(C));

  for (j = 1; j <= a_cols_mx(B); j++) {
    /* C(:,j) = A * B(:,j) */
    a_mvm (arena, A, a_col_mx(B,j), a_col_mx(C,j));
  }
}

void
a_mvm (arena_info *arena, AlphaMatrix A, AlphaVector x, AlphaVector y)
{
  int		j;
  int		i;
  AlphaVector 	col;
  AlphaVector	ytmp;

  assert (a_cols_mx(A) == a_size_vec(x));
  assert (a_rows_mx(A) == a_size_vec(y));

  ytmp = a_make_vector (arena, a_size_vec(y));

  for (i = 1; i <= a_size_vec(y); i++) {
    a_elt_vec(y,i) = 0;
  }

  for (j = 1; j <= a_cols_mx(A); j++) {
    /* y += x[j]*col */
    col = a_col_mx (A, j);
    a_copy_vec (ytmp, y);
    a_saxpy (1, ytmp, a_elt_vec(x,j), col, y);
  }

}

void
a_transpose (AlphaMatrix A, AlphaMatrix At)
{
  int		i,j;

  assert (a_cols_mx(A) == a_rows_mx(At));
  assert (a_rows_mx(A) == a_cols_mx(At));

  for (j = 1; j <= a_cols_mx(At); j++) {
    for (i = 1; i <= a_rows_mx(At); i++) {
      a_elt_mx (At,i,j) = a_elt_mx(A,j,i);
    }
  }
}

/* we allow i2(j2) to be one smaller that i1(j1) so that we can create
 * a matrix with zero rows(columns)
 */
void
a_extract (AlphaMatrix A, int i1, int i2, int j1, int j2, AlphaMatrix B)
{
  int		i, j, ib, jb;
  int		ni;
  int		nj;
  AlphaVector	colA, colB;

  ni = i2 - i1 + 1;
  nj = j2 - j1 + 1;

  assert (ni >= 0);
  assert (nj >= 0);
  assert (1 <= i1 && i2 <= a_rows_mx (A));
  assert (1 <= j1 && j2 <= a_cols_mx (A));

  assert (ni == a_rows_mx (B));
  assert (nj == a_cols_mx (B));

  for (j = j1, jb=1; j <= j2; j++, jb++) {
    colA = a_col_mx (A, j);
    colB = a_col_mx (B, jb);
    for (i = i1, ib=1; i <= i2; i++, ib++) {
      a_elt_vec (colB, ib) = a_elt_vec (colA, i);
    }
  }
}

void
a_replace (AlphaMatrix A, int i1, int i2, int j1, int j2, AlphaMatrix B)
{
  int		i, j, ib, jb;
  int		ni;
  int		nj;
  AlphaVector	colA, colB;

  assert (i1 <= i2);
  assert (j1 <= j2);
  assert (1 <= i1 && i2 <= a_rows_mx (A));
  assert (1 <= j1 && j2 <= a_cols_mx (A));

  ni = i2 - i1 + 1;
  nj = j2 - j1 + 1;

  assert (ni == a_rows_mx (B));
  assert (nj == a_cols_mx (B));

  for (j = j1, jb=1; j <= j2; j++, jb++) {
    colA = a_col_mx (A, j);
    colB = a_col_mx (B, jb);
    for (i = i1, ib=1; i <= i2; i++, ib++) {
      a_elt_vec (colA, i) = a_elt_vec (colB, ib);
    }
  }
}

void
a_add_cols (arena_info *arena, AlphaMatrix *pA, int c)
{
  xa_mx_t	*Anew;
  xa_mx_t	*A = *pA;
  int 		m = a_rows_mx (A);
  int		n = a_cols_mx (A);
  int		i;


  assert (c > 0);
  Anew = xa_alloc_matrix (arena, m, n + c);

  for (i = 1; i <= n; i++) {
    a_col_mx (Anew,i) = a_col_mx (A, i);
  }

  for (i = n+1; i <= n+c; i++) {
    a_col_mx (Anew, i) = a_make_vector (arena, m);
  }

  *pA = Anew;
}

int
a_is_zero_matrix (AlphaMatrix A)
{
    int m = a_rows_mx (A);
    int n = a_cols_mx (A);
    int i, j;

    for (i = 1; i <= m; i++) {
	for (j = 1; j <= n; j++) {
	    if (a_elt_mx (A, i, j) != 0) {
		return 0;
	    }
	}
    }

    return 1;
}

int
a_is_ident_matrix (AlphaMatrix A)
{
    int m = a_rows_mx (A);
    int n = a_cols_mx (A);
    int i, j;

    for (i = 1; i <= m; i++) {
	if (a_elt_mx(A,i,i) != 1) {
	    return 0;
	}

	for (j = 1; j <= n; j++) {
	    if (i == j) continue;

	    if (a_elt_mx (A, i, j) != 0) {
		return 0;
	    }
	}
    }

    return 1;
}

void		
a_trunc_cols (arena_info *arena, AlphaMatrix *pA, int c)
{
  xa_mx_t	*Anew;
  xa_mx_t	*A = *pA;
  int		n = a_cols_mx (A);
  int		m = a_rows_mx (A);
  int		newn = n - c;
  int		i;

  assert (c >= 0);

  if (c == 0) {
    return;
  }
  
  assert (c <= n);

  Anew = xa_alloc_matrix (arena, m, newn);
  for (i = 1; i <= newn; i++) {
    a_col_mx (Anew, i) = a_col_mx (A, i);
  }
  *pA = Anew;
}

/*
 *
 *   [ x y ] = [x y] * T
 *  T is 2 by 2
 *
 *   I.e. if T = [ t11 t12 ]
 *               [ t21 t22 ]
 *
 *  Then the result in xnew = xold * t11 + yold * t21
 *                     ynew = xold * t12 + yold * t22
 * */
static void
two_vec_transform (arena_info *arena, AlphaVector x, AlphaVector y, AlphaMatrix T)
{
  int		n = a_size_vec (x);
  AlphaMatrix 	srcm = a_make_matrix (arena, n, 2);
  AlphaMatrix 	dstm = a_make_matrix (arena, n, 2);

  assert (n == a_size_vec (y));

  a_copy_vec (a_col_mx(srcm,1), x);
  a_copy_vec (a_col_mx(srcm,2), y);
  a_matmult (arena, srcm, T, dstm);
  a_copy_vec (x, a_col_mx(dstm,1));
  a_copy_vec (y, a_col_mx(dstm,2));
}

static void
two_row_transform (arena_info *arena, AlphaMatrix A, int i1, int i2, AlphaMatrix T)
{
  int		n = a_cols_mx (A);
  int		m = a_rows_mx (A);
  AlphaMatrix	src, dest;
  int		j;

  assert (1 <= i1 && i1 <= m);
  assert (1 <= i2 && i2 <= m);
  assert (2 == a_cols_mx (T));
  assert (2 == a_rows_mx (T));

  src = a_make_matrix (arena, 2, n);
  dest = a_make_matrix (arena, 2, n);
  
  for (j = 1; j <= n; j++) {
    a_elt_mx (src, 1, j) = a_elt_mx (A, i1, j);
    a_elt_mx (src, 2, j) = a_elt_mx (A, i2, j);
  }

  a_matmult (arena, T, src, dest);

  for (j = 1; j <= n; j++) {
    a_elt_mx (A, i1, j) = a_elt_mx (dest, 1, j);
    a_elt_mx (A, i2, j) = a_elt_mx (dest, 2, j);
  }
}


/*
 * Given a matrix A compute the
 * echelon form: A * U = [ L 0 ]
 * IMPORTANT: A is being overwritten with [ L 0 ]
 */
void
a_lower_echelon (arena_info *arena, AlphaMatrix A, AlphaMatrix U, int *prank)
{
  int		i, k, j;
  int		a, b, s, t, d;
  int		n = a_cols_mx (A);
  int		m = a_rows_mx (A);
  AlphaMatrix	T = a_make_matrix (arena, 2, 2);

  assert (n == a_rows_mx (U));
  assert (n == a_cols_mx (U));

  /* 
   * we compute U by setting it to identity and then
   * applying to it the same column ops as to A
   */
  a_set2ident (U);

  k = 0;
  for (i = 1; (i <= m) && (k < n); i++) {
    for (j = k+2; j <= n; j++) {
      a = a_elt_mx (A, i, k+1);
      b = a_elt_mx (A, i, j);

      if (0 == a && 0 == b) 
	continue;

      d = a_egcd (a, b, &s, &t);

      /*
       * s * a + t * b == gcd(a,b) == d
       *
       * T = [ s -b/d ]
       *     [ t  a/d ]
       *
       * so that:   det(T) == 1
       * and
       *
       * [a b] * T = [d 0]
       *
       * This is sort of like an "integer Givens rotation"
       */

      a_elt_mx (T, 1, 1) = s;
      a_elt_mx (T, 1, 2) = -b/d;
      a_elt_mx (T, 2, 1) = t;
      a_elt_mx (T, 2, 2) = a/d;

      two_vec_transform (arena, a_col_mx(A,k+1), a_col_mx(A,j), T);
      two_vec_transform (arena, a_col_mx(U,k+1), a_col_mx(U,j), T);
    }
    if (0 != a_elt_mx(A,i,k+1)) {
      k = k+1;
    }
  }

  *prank = k;
}



int
a_solve_linear (arena_info *arena, AlphaMatrix A, AlphaVector B, 
		AlphaVector X0, AlphaMatrix H, int *pnull_dim)
{
  int		m = a_rows_mx (A);
  int		n = a_cols_mx (A);
  int		r;
  AlphaMatrix   U;
  AlphaMatrix	L;
  AlphaVector	Y0;
  int		j, j1, i, k;
  int		sum;
  int		result;

  assert (m >= 0);
  assert (n >= 0);
  assert (m == a_length_vec (B));
  assert (n == a_length_vec (X0));
  assert (n == a_cols_mx (H));
  assert (n == a_rows_mx (H));

  if (n == 0 && m > 0) {
      /* empty_matrix * x = non_empty
	 can not have solutions */
      return 0;
  }

  U = a_make_matrix (arena, n, n);
  L = a_make_matrix (arena, m, n);
  Y0 = a_make_vector (arena, n);
  a_copy_mx (L, A);
  a_lower_echelon (arena, L, U, &r);

  for (j = r+1, j1 = 1; j <= n; j++, j1++) {
    a_copy_vec (a_col_mx(H,j1), a_col_mx(U,j));
  }

  /* now columns H(:,1) to H(:,n-r) contain the null space basis */

  *pnull_dim = n - r;

  if (r == 0) {
      if (a_is_zero_vec (B)) {
	  a_set_vec (X0, 0);
	  result = 1;
      }
      else {
	  result = 0;
      }
  }
  else {

    if (a_solve_echelon (L, B, Y0)) {

#if !NDEBUG
      { 
	AlphaVector should_be_B = a_make_vector (arena, a_length_vec (B));
	AlphaVector diff = a_make_vector (arena, a_length_vec (B));
	a_mvm (arena, L, Y0, should_be_B);
	a_saxpy (1, B, -1, should_be_B, diff);
	assert (a_is_zero_vec (diff));
      }
#endif

      result = 1;
      a_mvm (arena, U, Y0, X0);
    }
    else {
      result = 0;
    }

  }

  return result;  
}

/*
 * find a particular solution to a lower schelon system
 */
int
a_solve_echelon (AlphaMatrix L, AlphaVector B, AlphaVector X)
{
  int		m = a_rows_mx (L);
  int		n = a_cols_mx (L);
  int		j, i, k;
  int		sum;
  int		result;

  assert (m > 0);
  assert (n > 0);
  assert (m == a_length_vec (B));
  assert (n == a_length_vec (X));

  k = 0;
  result = 1;
  for (i = 1; (0 != result) && (i <= m); i++) {
    if (a_elt_mx(L,i,k+1) != 0) {
      /* we can compute a next element of X */
      k = k+1;
      sum = 0;
      for (j = 1; j < k; j++) {
	sum += a_elt_mx(L,i,j) * a_elt_vec (X,j);
      }
      sum = a_elt_vec(B,i) - sum;
      if (0 != sum % a_elt_mx(L,i,k)) { /* no integer solution here */
	result = 0;
      }
      else {
	a_elt_vec (X,k) = sum / a_elt_mx(L,i,k);
      }
    }
    else {
      /* dependent row -> just check for consistency */
      sum = 0;
      for (j = 1; j <= k; j++) {
	sum += a_elt_mx(L,i,j) * a_elt_vec (X,j);
      }
      if (sum != a_elt_vec (B,i)) {
	result = 0;
      }
    }
  }

  return result;  
}

/*
 * Compute the decomposition  U*A*V = Sigma
 *
 * P.S. Sigma is diagonal of the same rank as A, but I am not
 * sure if it exactly SNF
 */
void
a_snf (arena_info *arena, AlphaMatrix A, AlphaMatrix U, AlphaMatrix V, int *prank)
{
  int		m = a_rows_mx (A);
  int		n = a_cols_mx (A);
  AlphaMatrix   T;
  int		i, i0, j0, j, k;
  int		nz_a;
  int		a, aa, b, d, s, t;
  int		was_nz;

  assert (m == a_rows_mx (U));
  assert (m == a_cols_mx (U));
  assert (n == a_rows_mx (V));
  assert (n == a_cols_mx (V));

  a_set2ident (U);
  a_set2ident (V);

  if (m == 0 || n == 0) {
    *prank = 0;
    return;
  }

  T = a_make_matrix (arena, 2, 2);

  /*
   * we compute U and V by applying the same operations to an identity 
   * as to A. 
   */
  
  for (k = 1; k <= m && k <= n; k++) {
    /* we work on submatrix A(k:m,k:n) */

    /* find a non-zero entry */
    nz_a = 0;
    i0 = j0 = k;
    for (i = k; (0 == nz_a) && (i <= m); i++) {
	for (j = k; (j <= n); j++) {
	    if (0 != a_elt_mx(A,i,j)) {
		nz_a = a_elt_mx(A,i,j);
		i0 = i;
		j0 = j;
		break; /* out of the inner loop */
	    }
	}
    }

    if (nz_a == 0) { /* we are done ! */
      break;
    }

    /* bring the smallest elt in the "corner" */
    a_swap_row (A, k, i0);
    a_swap_row (U, k, i0);
    a_swap_col (A, k, j0);
    a_swap_col (V, k, j0);

    /* make sure it is positive */
    if (a_elt_mx (A, k, k) < 0) {
      a_neg_vec (a_col_mx (A, k));
      a_neg_vec (a_col_mx (V, k));
    }

    /* keep reducing k-th row and k-th column until the diagonal 
       element does not change anymore */ 
    do {
      was_nz = 0;

      a = a_elt_mx (A, k, k);

      /* reduce k-th row */
      for (j = k+1; j <= n; j++) {
	b = a_elt_mx (A, k, j);
	if (0 == b) {
	  continue;
	}
	was_nz = 1;
	d = a_egcd (a, b, &s, &t);

	/* see comments about T in a_lower_echelon() */
	a_elt_mx (T, 1, 1) = s;
	a_elt_mx (T, 1, 2) = -b/d;
	a_elt_mx (T, 2, 1) = t;
	a_elt_mx (T, 2, 2) = a/d;

	two_vec_transform (arena, a_col_mx(A,k), a_col_mx(A,j), T);
	two_vec_transform (arena, a_col_mx(V,k), a_col_mx(V,j), T);

	a = d;
      }

#if !NDEBUG
      aa = a_elt_mx (A, k, k);
      assert (a == aa);
#endif

      /* reduce k-th column */
      for (j = k+1; j <= m; j++) {
	b = a_elt_mx (A, j, k);
	if (0 == b) {
	  continue;
	}
	was_nz = 1;
	d = a_egcd (a, b, &s, &t);

	/* careful! use T transposed here */
	a_elt_mx (T, 1, 1) = s;
	a_elt_mx (T, 2, 1) = -b/d;
	a_elt_mx (T, 1, 2) = t;
	a_elt_mx (T, 2, 2) = a/d;

	two_row_transform (arena, A, k, j, T);
	two_row_transform (arena, U, k, j, T);

	a = d;
      }
      
#if !NDEBUG
      aa = a_elt_mx (A, k, k);
      assert (a == aa);
#endif

    }
    while (was_nz);
  }

  *prank = k-1;

#if !NDEBUG
  {/* verify that 1s come first in the SNF */
    int non_unit;
    non_unit = 0;
    for (j = 1; j < k; j++) {
      a = a_elt_mx (A, j, j);
      assert (  !non_unit || (a != 1) );
      non_unit = (a != 1);
    }
  }
#endif
}

/*
 * solve U*x = b, where U is unit upper triangular
 * (diagonal can be +- 1)
 */
void
a_unit_backsolve (AlphaMatrix U, AlphaVector b, AlphaVector x)
{
    int n = a_rows_mx (U);
    int i, j, sum, d;

    assert(a_cols_mx(U) == n);
    assert(a_length_vec(b) == n);
    assert(a_length_vec(x) == n);

    for (i = n; i >= 1; i--) {
	d = a_elt_mx(U,i,i);
	assert (1 == d || -1 == d);
	sum = 0;
	for (j = i+1; j <= n; j++) {
	    sum += a_elt_mx(U,i,j) * a_elt_vec(x,j);
	}
	a_elt_vec(x,i) = (a_elt_vec(b,i) - sum) * d;
    }
}

/*
 * Given a matrix A and a unit-upper-tringular U compute 
 * U^{-1} * A by doing a sequence of backsolves
 */
void
a_block_unit_backsolve (AlphaMatrix U, AlphaMatrix B, AlphaMatrix X)
{
    int n = a_rows_mx (U);
    int m = a_cols_mx (B);
    int k;
    
    assert(a_cols_mx(U) == n);
    assert(a_rows_mx(B) == n);
    assert(a_rows_mx(X) == n);
    assert(a_cols_mx(X) == m);

    for (k = 1; k <= m; k++) {
	a_unit_backsolve (U, a_col_mx(B,k), a_col_mx(X,k));
    }
}

void
a_read_matrix (FILE *fp, AlphaMatrix A)
{
  int		i;
  int		j;
  int		rc;

  for (i = 1; i <= a_rows_mx(A); i++) {
    for (j = 1; j <= a_cols_mx(A); j++) {
      rc = fscanf (fp, "%d", & a_elt_mx (A, i, j));
      assert (1 == rc);
    }
  }
}

void
a_dump_matrix (FILE *fp, AlphaMatrix A)
{
  int		i;
  int		j;

  fprintf (fp, "[\n");
  for (i = 1; i <= a_rows_mx(A); i++) {
    for (j = 1; j <= a_cols_mx(A); j++) {
      fprintf (fp, "%d ", a_elt_mx (A, i, j));
    }
    putc ('\n', fp);
  }
  fprintf (fp, "]\n");

#ifndef NDEBUG
  fflush (fp);
#endif
}

void
a_print_flat_matrix (FILE *fp, AlphaMatrix A)
{
  int i, j;
  int m = a_rows_mx (A);
  int n = a_cols_mx (A);
  fprintf (fp, "[");
  for (i = 1; i <= m; i++) {
    for (j = 1; j <= n; j++) {
      fprintf (fp, " %d", a_elt_mx (A, i, j));
    }
    if (i != m) {
      fprintf (fp, ";");
    }
  }
  fprintf (fp, "]");
}

