mutable-matrix.cpp

Go to the documentation of this file.

// Copyright 2004 Michael E. Stillman
 
#include "interface/mutable-matrix.h"
 
#include <assert.h>
#include <algorithm>
#include <vector>
 
#include "LLL.hpp"
#include "aring-zzp-ffpack.hpp"
#include "buffer.hpp"
#include "error.h"
#include "exceptions.hpp"
#include "finalize.hpp"
#include "fractionfreeLU.hpp"
#include "interface/gmp-util.h"
#include "interface/random.h"
#include "mat.hpp"
#include "matrix.hpp"
#include "relem.hpp"
#include "ring.hpp"
#include "ringelem.hpp"
#include "util.hpp"
 
MutableMatrix *IM2_MutableMatrix_identity(const Ring *R,
                                          int n,
                                          M2_bool is_dense)
{
  if (n < 0)
    {
      ERROR("expected nonnegative integer");
      return nullptr;
    }
  size_t nrows = static_cast<size_t>(n);
  return internMutableMatrix(MutableMatrix::identity(R, nrows, is_dense));
}
 
MutableMatrix /* or null */ *IM2_MutableMatrix_make(const Ring *R,
                                                    int nrows,
                                                    int ncols,
                                                    M2_bool is_dense)
{
  if (nrows < 0 || ncols < 0) {
    ERROR("expected nonnegative integers");
    return nullptr;
  }
 
  size_t nr = static_cast<size_t>(nrows);
  size_t nc = static_cast<size_t>(ncols);
  //  return R->makeMutableMatrix(nr,nc,is_dense);
  return internMutableMatrix(MutableMatrix::zero_matrix(R, nr, nc, is_dense));
}
 
MutableMatrix *IM2_MutableMatrix_from_matrix(const Matrix *M, M2_bool is_dense)
{
  return internMutableMatrix(MutableMatrix::from_matrix(M, is_dense));
}
 
const Matrix *IM2_MutableMatrix_to_matrix(const MutableMatrix *M)
{
  return M->to_matrix();
}
 
M2_string IM2_MutableMatrix_to_string(const MutableMatrix *M)
{
  buffer o;
  M->text_out(o);
  return o.to_string();
}
 
unsigned int rawMutableMatrixHash(const MutableMatrix *M) { return M->hash(); }
MutableMatrix /* or null */ *rawMutableMatrixPromote(const Ring *S,
                                                     const MutableMatrix *f)
{
  (void) S;
  (void) f;
  // Given a natural map i : R --> S
  // f is a matrix over R.
  // returns a matrix over S.
  ERROR("MutableMatrix promote not implemented yet");
  return nullptr;
#if 0  
  auto result = MutableMatrix::zero_matrix(S, f->n_rows(), f->n_cols(), f->is_dense());
  return result;
#endif
}
 
MutableMatrix /* or null */ *rawMutableMatrixLift(int *success_return,
                                                  const Ring *R,
                                                  const MutableMatrix *f)
{
  // Given a natural map i : R --> S
  // f is a matrix over S.
  // returns a matrix over R.
 
  // ERROR("MutableMatrix lift not implemented yet");
  (void) R;
  (void) f;
  *success_return = 0;
  return nullptr;
#if 0
  auto result = MutableMatrix::zero_matrix(R, f->n_rows(), f->n_cols(), f->is_dense());
  *success_return = 1;
  return result;
#endif
}
 
int IM2_MutableMatrix_n_rows(const MutableMatrix *M)
{
  size_t nrows = M->n_rows();
  return static_cast<int>(nrows);
}
 
int IM2_MutableMatrix_n_cols(const MutableMatrix *M)
{
  size_t ncols = M->n_cols();
  return static_cast<int>(ncols);
}
 
void rawMutableMatrixFillRandom(MutableMatrix *M, long nelems)
{
  int nrows = static_cast<int>(M->n_rows());
  int ncols = static_cast<int>(M->n_cols());
  const Ring *R = M->get_ring();
  for (long i = 0; i < nelems; i++)
    {
      int r = rawRandomInt(nrows);
      int c = rawRandomInt(ncols);
      ring_elem a = R->random();
      if (!R->is_zero(a)) M->set_entry(r, c, R->random());
    }
}
 
void rawMutableMatrixFillRandomDensity(MutableMatrix *M,
                                       double density,
                                       int special_type)
/* special_type: 0 is general, 1 is (strictly) upper triangular. */
{
  bool doing_fraction = false;
  int threshold = 0;
 
  int nrows = static_cast<int>(M->n_rows());
  int ncols = static_cast<int>(M->n_cols());
  const Ring *R = M->get_ring();
 
  if (density != 1.0)
    {
      doing_fraction = true;
      threshold = static_cast<int>(density * 10000);
    }
 
  if (special_type == 0)
    {
      for (int i = 0; i < ncols; i++)
        for (int j = 0; j < nrows; j++)
          {
            if (doing_fraction)
              {
                int32_t u = rawRandomInt((int32_t)10000);
                if (u > threshold) continue;
              }
            ring_elem a = R->random();
            if (!R->is_zero(a)) M->set_entry(j, i, a);
          }
    }
  else if (special_type == 1)
    {
      for (int i = 0; i < ncols; i++)
        {
          int top = (i >= nrows ? nrows : i);
          for (int j = 0; j < top; j++)
            {
              if (doing_fraction)
                {
                  int32_t u = rawRandomInt((int32_t)10000);
                  if (u > threshold) continue;
                }
              ring_elem a = R->random();
              if (!R->is_zero(a)) M->set_entry(j, i, a);
            }
        }
    }
}
 
const RingElement /* or null */ *
IM2_MutableMatrix_get_entry(const MutableMatrix *M, int r, int c)
{
  if (r < 0 || r >= M->n_rows())
    {
      ERROR("matrix row index %d out of range 0 .. %d", r, M->n_rows() - 1);
      return nullptr;
    }
  if (c < 0 || c >= M->n_cols())
    {
      ERROR("matrix column index %d out of range 0 .. %d", c, M->n_cols() - 1);
      return nullptr;
    }
  ring_elem result;
  M->get_entry(r, c, result);
  return RingElement::make_raw(M->get_ring(), result);
}
 
engine_RawRingElementArrayArrayOrNull
IM2_MutableMatrix_get_entries(const MutableMatrix *M)
{
  try
    {
      int ncols = M->n_cols();
      int nrows = M->n_rows();
      if(nrows < 0 || ncols < 0)
        {
          ERROR("internal error: matrix has a negative size %d by %d",
                nrows,
                ncols);
          return nullptr;
        }
      engine_RawRingElementArrayArray entries =
          getmemarraytype(engine_RawRingElementArrayArray, nrows);
      entries->len = nrows;
      for(int r = 0; r < nrows; r++)
        {
          engine_RawRingElementArray currRow =
              getmemarraytype(engine_RawRingElementArray, ncols);
          currRow->len = ncols;
          entries->array[r] = currRow;
        }
      for(int r = 0; r < nrows; r++)
        {
          for(int c = 0; c < ncols; c++)
            {
              ring_elem result;
              M->get_entry(r, c, result);
              entries->array[r]->array[c] =
                  RingElement::make_raw(M->get_ring(), result);
            }
        }
      return entries;
    } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
    }
  return nullptr;
}
 
 
M2_bool IM2_MutableMatrix_set_entry(MutableMatrix *M,
                                    int r,
                                    int c,
                                    const RingElement *a)
{
  const Ring *R = M->get_ring();
  if (R != a->get_ring())
    {
      ERROR("expected same ring");
      return 0;
    }
  if (r < 0 || r >= M->n_rows())
    {
      ERROR("row index %d out of range 0..%d", r, M->n_rows() - 1);
      return 0;
    }
  if (c < 0 || c >= M->n_cols())
    {
      ERROR("column index %d out of range 0..%d", c, M->n_cols() - 1);
      return 0;
    }
 
  M->set_entry(r, c, a->get_value());
  return 1;
}
 
M2_bool IM2_MutableMatrix_row_swap(MutableMatrix *M, int i, int j)
/* swap rows: row(i) <--> row(j) */
{
  if (i < 0 || j < 0 || i >= M->n_rows() || j >= M->n_rows())
    {
      ERROR("row index out of range");
      return 0;
    }
 
  M->interchange_rows(i, j);
  return 1;
}
 
M2_bool IM2_MutableMatrix_column_swap(MutableMatrix *M, int i, int j)
/* swap columns: column(i) <--> column(j) */
{
  if (i < 0 || j < 0 || i >= M->n_cols() || j >= M->n_cols())
    {
      ERROR("column index out of range");
      return 0;
    }
 
  M->interchange_columns(i, j);
  return 1;
}
 
M2_bool IM2_MutableMatrix_row_operation(MutableMatrix *M,
                                        int i,
                                        const RingElement *r,
                                        int j,
                                        M2_bool opposite_mult)
/* row(i) <- row(i) + r * row(j) */
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (R != r->get_ring())
    {
      ERROR("expected same ring");
      return 0;
    }
  if (i < 0 || j < 0 || i >= M->n_rows() || j >= M->n_rows())
    {
      ERROR("row index out of range");
      return 0;
    }
 
  M->row_op(i, r->get_value(), j);
  return 1;
}
 
M2_bool IM2_MutableMatrix_column_operation(MutableMatrix *M,
                                           int i,
                                           const RingElement *r,
                                           int j,
                                           M2_bool opposite_mult)
/* column(i) <- column(i) + r * column(j) */
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (R != r->get_ring())
    {
      ERROR("expected same ring");
      return 0;
    }
  if (i < 0 || j < 0 || i >= M->n_cols() || j >= M->n_cols())
    {
      ERROR("column index out of range");
      return 0;
    }
 
  M->column_op(i, r->get_value(), j);
  return 1;
}
 
M2_bool IM2_MutableMatrix_row_scale(MutableMatrix *M,
                                    const RingElement *r,
                                    int i,
                                    M2_bool opposite_mult)
/* row(i) <- r * row(i) */
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (R != r->get_ring())
    {
      ERROR("expected same ring");
      return 0;
    }
  if (i < 0 || i >= M->n_rows())
    {
      ERROR("row index out of range");
      return 0;
    }
  M->scale_row(i, r->get_value());
  return 1;
}
 
M2_bool IM2_MutableMatrix_column_scale(MutableMatrix *M,
                                       const RingElement *r,
                                       int i,
                                       M2_bool opposite_mult)
/* column(i) <- r * column(i) */
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (R != r->get_ring())
    {
      ERROR("expected same ring");
      return 0;
    }
  if (i < 0 || i >= M->n_cols())
    {
      ERROR("column index out of range");
      return 0;
    }
  M->scale_column(i, r->get_value());
  return 1;
}
 
M2_bool IM2_MutableMatrix_insert_columns(MutableMatrix *M, int i, int n_to_add)
/* Insert n_to_add columns directly BEFORE column i. */
{
  return M->insert_columns(i, n_to_add);
}
 
M2_bool IM2_MutableMatrix_insert_rows(MutableMatrix *M, int i, int n_to_add)
/* Insert n_to_add rows directly BEFORE row i. */
{
  return M->insert_rows(i, n_to_add);
}
 
M2_bool IM2_MutableMatrix_delete_columns(MutableMatrix *M, int i, int j)
/* Delete columns i .. j from M */ { return M->delete_columns(i, j); }
M2_bool IM2_MutableMatrix_delete_rows(MutableMatrix *M, int i, int j)
/* Delete rows i .. j from M */ { return M->delete_rows(i, j); }
M2_bool IM2_MutableMatrix_column_2by2(MutableMatrix *M,
                                      int c1,
                                      int c2,
                                      const RingElement *a1,
                                      const RingElement *a2,
                                      const RingElement *b1,
                                      const RingElement *b2,
                                      M2_bool opposite_mult)
/* column(c1) <- a1 * column(c1) + a2 * column(c2)
   column(c2) <- b1 * column(c1) + b2 * column(c2)
*/
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (a1->get_ring() != R || a2->get_ring() != R || b1->get_ring() != R ||
      b2->get_ring() != R)
    {
      ERROR("expected elements in the same ring");
      return 0;
    }
  return M->column2by2(c1,
                       c2,
                       a1->get_value(),
                       a2->get_value(),
                       b1->get_value(),
                       b2->get_value());
}
 
M2_bool IM2_MutableMatrix_row_2by2(MutableMatrix *M,
                                   int r1,
                                   int r2,
                                   const RingElement *a1,
                                   const RingElement *a2,
                                   const RingElement *b1,
                                   const RingElement *b2,
                                   M2_bool opposite_mult)
/* row(r1) <- a1 * row(r1) + a2 * row(r2)
   row(r2) <- b1 * row(r1) + b2 * row(r2)
*/
{
  (void) opposite_mult;
  const Ring *R = M->get_ring();
  if (a1->get_ring() != R || a2->get_ring() != R || b1->get_ring() != R ||
      b2->get_ring() != R)
    {
      ERROR("expected elements in the same ring");
      return 0;
    }
  return M->row2by2(r1,
                    r2,
                    a1->get_value(),
                    a2->get_value(),
                    b1->get_value(),
                    b2->get_value());
}
 
M2_bool IM2_MutableMatrix_sort_columns(MutableMatrix *M, int lo, int hi)
/* Returns false if M is not mutable, or lo, or hi are out of range */
{
  (void) M;
  (void) lo;
  (void) hi;
  ERROR("not re-implemented yet");
  return false;
}
 
M2_bool IM2_MutableMatrix_row_permute(MutableMatrix *M,
                                      int start,
                                      M2_arrayint perm)
/* if perm = [p0 .. pr], then row(start + i) --> row(start + pi), and
   all other rows are unchanged.  p0 .. pr should be a permutation of 0..r */
{
  size_t nrows = M->n_rows();
  if (start < 0 || start + perm->len > nrows)
    {
      ERROR("row indices out of range");
      return false;
    }
  for (int i = 0; i < perm->len; i++)
    {
      int r = start + perm->array[i];
      if (r < 0 || r >= nrows)
        {
          ERROR("row indices out of range");
          return false;
        }
    }
  return M->row_permute(start, perm);
}
 
M2_bool IM2_MutableMatrix_column_permute(MutableMatrix *M,
                                         int start,
                                         M2_arrayint perm)
/* if perm = [p0 .. pr], then column(start + i) --> column(start + pi), and
   all other rows are unchanged.  p0 .. pr should be a permutation of 0..r */
{
  size_t ncols = M->n_cols();
  if (start < 0 || start + perm->len > ncols)
    {
      ERROR("column indices out of range");
      return false;
    }
  for (int i = 0; i < perm->len; i++)
    {
      int r = start + perm->array[i];
      if (r < 0 || r >= ncols)
        {
          ERROR("column indices out of range");
          return false;
        }
    }
  return M->column_permute(start, perm);
}
 
const RingElement *IM2_Matrix_dot_product(const MutableMatrix *M,
                                          int c1,
                                          int c2)
{
  ring_elem a;
  M->dot_product(c1, c2, a);
  return RingElement::make_raw(M->get_ring(), a);
}
 
M2_bool IM2_MutableMatrix_is_zero(const MutableMatrix *M)
{
  return M->is_zero();
}
 
M2_bool IM2_MutableMatrix_is_equal(const MutableMatrix *M,
                                   const MutableMatrix *N)
/* This checks that the entries of M,N are the same */
{
  return M->is_equal(N);
}
 
MutableMatrix /* or null */ *rawMutableMatrixTranspose(MutableMatrix *M)
{
  try
    {
      return internMutableMatrix(M->transpose());
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
M2_bool rawMutableMatrixIsDense(const MutableMatrix *M)
{
  return M->is_dense();
}
 
MutableMatrix *IM2_MutableMatrix_copy(MutableMatrix *M, M2_bool prefer_dense)
{
  return internMutableMatrix(M->copy(prefer_dense));
}
 
M2_bool IM2_MutableMatrix_set_values(MutableMatrix *M,
                                     M2_arrayint rows,
                                     M2_arrayint cols,
                                     engine_RawRingElementArray values)
{
  return M->set_values(rows, cols, values);
}
 
M2_bool IM2_MutableMatrix_reduce_by_pivots(MutableMatrix *M)
/* Using row and column operations, use unit pivots to reduce the matrix */
{
  M->reduce_by_pivots();
  return true;
}
 
MutableMatrix /* or null */ *IM2_MutableMatrix_submatrix(const MutableMatrix *M,
                                                         M2_arrayint rows,
                                                         M2_arrayint cols)
{
  return internMutableMatrix(M->submatrix(rows, cols));
}
 
MutableMatrix /* or null */ *IM2_MutableMatrix_submatrix1(
    const MutableMatrix *M,
    M2_arrayint cols)
{
  return internMutableMatrix(M->submatrix(cols));
}
 
/*******************************
 ** Computations ***************
 *******************************/
 
M2_arrayintOrNull IM2_FF_LU(MutableMatrix *M)
{
  return FF_LUComputation::DO(M);
}
 
#include <fplll-interface.hpp>
#include "ntl-interface.hpp"
 
M2_bool rawLLL(MutableMatrix *M,
               MutableMatrix /* or null */ *U,
               gmp_QQ threshold,
               int strategy)
{
  if (strategy == 0)
    {
      return LLLoperations::LLL(M, U, threshold);
    }
 
  if (strategy == 4)  // fplll, pfush...
    {
      return fp_LLL(M, U, strategy);
    }
 
  long a = mpz_get_si(mpq_numref(threshold));
  long b = mpz_get_si(mpq_denref(threshold));
  return ntl_LLL(M, U, a, b, strategy);
}
 
M2_bool IM2_SmithNormalForm(MutableMatrix *M)
{
  (void) M;
#ifdef DEVELOPMENT
#warning "implement smith"
#endif
  ERROR("not implemented yet");
  return 0;
}
 
M2_bool IM2_HermiteNormalForm(MutableMatrix *M)
{
  (void) M;
#ifdef DEVELOPMENT
#warning "implement hermite"
#endif
  ERROR("not implemented yet");
  return 0;
}
 
/***************************************************
 ***** Lapack routines for dense mutable matrices **
 ***************************************************/
 
/* Each of the following routines accepts honest MutableMatrix arguments,
   and returns false if there is an error.  The return values are placed into
   some of the (already existing) parameters of the routine */
 
M2_arrayintOrNull rawLU(const MutableMatrix *A,
                        MutableMatrix *L,
                        MutableMatrix *U)
{
  try
    {
      return A->LU(L, U);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
M2_arrayintOrNull rawLUincremental(M2_arrayintOrNull P,
                                   MutableMatrix *LU,
                                   const MutableMatrix* v,
                                   int i)
{
  try
    {
      // FIXME: can we not allocate new permutation array?
      std::vector<size_t> perm = M2_arrayint_to_stdvector<size_t>(P);
      return LU->LUincremental(perm, v, i);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
void rawTriangularSolve(MutableMatrix *Lv, /* modified in routine */
                        MutableMatrix *x, /* modified in routine */
                        int m,
                        int strategy)
{
  try
    {
      Lv->triangularSolve(x, m, strategy);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
  }
}

 
#if 0
M2_bool rawSolve(MutableMatrix *A,
                 MutableMatrix *b,
                 MutableMatrix *x)
{
  /* Check: A, b, x all have the same ring, either RR or CC */
  /* Check: if all of these are dense mutable matrices, then
     call the correct routine */
  /* Otherwise: give error:
     OR: make mutable matrices of the correct size, call the correct routine
     and afterwards, copy to x. */
  try {
    const Ring *R = A->get_ring();
    if (R != b->get_ring() || R != x->get_ring())
      {
        ERROR("expected matrices with same base ring");
        return false;
      }
    if (A->n_rows() != b->n_rows())
      {
        ERROR("expected matrices with the same number of rows");
        return false;
      }
    return A->solve(b,x);
  }
  catch (const exc::engine_error& e) {
    ERROR(e.what());
    return false;
  }
}
#endif
 
M2_bool rawEigenvalues(MutableMatrix *A,
                       MutableMatrix *eigenvalues,
                       M2_bool is_symm_or_hermitian)
{
  try
    {
      return A->eigenvalues(eigenvalues, is_symm_or_hermitian);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
M2_bool rawEigenvectors(MutableMatrix *A,
                        MutableMatrix *eigenvalues,
                        MutableMatrix *eigenvectors,
                        M2_bool is_symm_or_hermitian)
{
  try
    {
      return A->eigenvectors(eigenvalues, eigenvectors, is_symm_or_hermitian);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
M2_bool rawSVD(MutableMatrix *A,
               MutableMatrix *Sigma,
               MutableMatrix *U,
               MutableMatrix *VT,
               M2_bool use_divide_and_conquer)
{
  try
    {
      return A->SVD(Sigma, U, VT, use_divide_and_conquer);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
M2_bool rawLeastSquares(MutableMatrix *A,
                        MutableMatrix *b,
                        MutableMatrix *x, /* return value: argument modified */
                        M2_bool assume_full_rank)
/* Case 1: A is a dense matrix over RR.  Then so are b,x.
   Case 2: A is a dense matrix over CC.  Then so are b,x. */
{
  try
    {
      const Ring *R = A->get_ring();
      if (R != b->get_ring() || R != x->get_ring())
        {
          ERROR("expected matrices with same base ring");
          return false;
        }
      if (A->n_rows() != b->n_rows())
        {
          ERROR("expected matrices with the same number of rows");
          return false;
        }
      return A->least_squares(b, x, assume_full_rank);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
M2_bool rawQR(const MutableMatrix *A,
              MutableMatrix *Q,
              MutableMatrix *R,
              M2_bool return_QR)
{
  try
    {
      return A->QR(Q, R, return_QR);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}

// Fast linear algebra routines //
 
long rawLinAlgRank(MutableMatrix *M)
{
  try
    {
      return M->rank();
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return -1;
  }
}
 
const RingElement *rawLinAlgDeterminant(MutableMatrix *A)
{
  try
    {
      return A->determinant();
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
MutableMatrix *rawLinAlgInverse(MutableMatrix *A)
{
  try
    {
      MutableMatrix *B = A->invert();
      if (B==nullptr) ERROR("matrix not invertible");
      return internMutableMatrix(B);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
MutableMatrix *rawLinAlgRREF(MutableMatrix *A)
{
  try
    {
      return internMutableMatrix(A->rowReducedEchelonForm());
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
M2_arrayintOrNull rawLinAlgRankProfile(MutableMatrix *A, M2_bool row_profile)
{
  try
    {
      return A->rankProfile(row_profile);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
MutableMatrix *rawLinAlgNullSpace(MutableMatrix *A)
{
  try
    {
      return internMutableMatrix(A->nullSpace());
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
MutableMatrix *rawLinAlgSolve(const MutableMatrix *A,
                              const MutableMatrix *B,
                              int *success)
{
  try
    {
      *success = 1;
      MutableMatrix *result = A->solveLinear(B);
      if (result != nullptr)
        {
          return internMutableMatrix(result);
        }
      else
        {
          return nullptr;
        }
  } catch (const exc::engine_error& e)
    {
      *success = 0;
      ERROR(e.what());
      return nullptr;
  }
}
 
MutableMatrix *rawLinAlgSolveInvertible(const MutableMatrix *A,
                                        const MutableMatrix *B,
                                        int *success)
{
  try
    {
      *success = 1;
      MutableMatrix *result = A->solveInvertible(B);
      if (result != nullptr)
        {
          return internMutableMatrix(result);
        }
      else
        {
          return nullptr;
        }
  } catch (const exc::engine_error& e)
    {
      *success = 0;
      ERROR(e.what());
      return nullptr;
  }
}
 
M2_bool rawLinAlgAddMult(MutableMatrix *C,
                         const MutableMatrix *A,
                         const MutableMatrix *B)
{
  try
    {
      C->addMultipleTo(A, B);
      return true;
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
M2_bool rawLinAlgSubMult(MutableMatrix *C,
                         const MutableMatrix *A,
                         const MutableMatrix *B)
{
  try
    {
      C->subtractMultipleTo(A, B);
      return true;
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return false;
  }
}
 
/* Note: the following routine is *not* called by the front end, as
   the * operator for MutableMatrix is implemented directly in d/engine.dd
*/
MutableMatrix * /* or null */ rawLinAlgMult(const MutableMatrix *A,
                                            const MutableMatrix *B)
{
  try
    {
      return internMutableMatrix(A->mult(B));
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
engine_RawRingElementArray convertRingelemsToArray(
    const Ring *R,
    std::vector<M2::ARingZZpFFPACK::ElementType> &elems)
{
  size_t len = elems.size();
  engine_RawRingElementArray result =
      getmemarraytype(engine_RawRingElementArray, len);
  result->len = static_cast<int>(len);
  for (size_t i = 0; i < len; i++)
    result->array[i] = RingElement::make_raw(R, static_cast<int>(elems[i]));
 
  return result;
}
 
engine_RawRingElementArrayOrNull rawLinAlgCharPoly(MutableMatrix *A)
// returns an array whose coefficients give the characteristic polynomial of the
// square matrix A
{
  (void) A;
#if 1
#if 0
  const Ring *R = A->get_ring();
  typedef DMat<M2::ARingZZpFFPACK> DMatZZp;
  MutableMat<DMatZZp> *B = A->cast_to_MutableMat<DMatZZp>();
  if (B == 0)
    {
      ERROR("expected a dense mutable matrix over the ffpack finite field");
      return 0;
    }
  M2::ARingZZpFFPACK::ElementType *elemsA = B->get_Mat()->array();
  std::vector<M2::ARingZZpFFPACK::ElementType> charpoly;
 
  // CharPoly isn't there any more (?)
  FFPACK::CharPoly(B->get_Mat()->ring().field(), charpoly, A->n_rows(),
  elemsA, A->n_rows());
 
  for (size_t i = 0; i < charpoly.size(); i++) std::cout << charpoly[i] << " ";
  std::cout << std::endl;
  return convertRingelemsToArray(R, charpoly);
#else
  return nullptr;
#endif
#else
  ERROR("not implemented: configure M2 with --enable-ffpack-fflas");
  return 0;
#endif
}
 
engine_RawRingElementArrayOrNull rawLinAlgMinPoly(MutableMatrix *A)
// returns an array whose coefficients give the minimal polynomial of the square
// matrix A
{
  (void) A;
#if 1
#if 0
  const Ring *R = A->get_ring();
  typedef DMat<M2::ARingZZpFFPACK> DMatZZp;
  MutableMat<DMatZZp> *B = A->cast_to_MutableMat<DMatZZp>();
  if (B == 0)
    {
      ERROR("expected a dense mutable matrix over the ffpack finite field");
      return 0;
    }
  typedef M2::ARingZZpFFPACK::ElementType Element;
  typedef std::vector<Element> Polynomial;
 
  Element *elemsA = B->get_Mat()->array();
  size_t n = B->n_rows();
  Polynomial minpoly(n);
 
  // the following are uninitialized!  This is scratch space for use in the
  // algorithm, apparently...
  Element *X = new Element[n * (n + 1)];
  size_t *P = new size_t[n];
 
  // this is in ffpack 2.3 but not in 2.2, and it might be the wrong name now, was MinPoly before
  FFPACK::Protected::Hybrid_KGF_LUK_MinPoly(
       B->get_Mat()->ring().field(), minpoly, n, elemsA, n, X, n, P);
 
  delete[] P;
  delete[] X;
 
  for (size_t i = 0; i < minpoly.size(); i++) std::cout << minpoly[i] << " ";
  std::cout << std::endl;
  return convertRingelemsToArray(R, minpoly);
#else
  return nullptr;
#endif
#else
  ERROR("not implemented: configure M2 with --enable-ffpack-fflas");
  return 0;
#endif
}

// Support for RRR and CCC operations //
 
const Matrix /* or null */ *rawMatrixClean(gmp_RR epsilon, const Matrix *M)
{
  try
    {
      if (M->get_ring()->get_precision() == 0)
        {
          ERROR("expected ring over an RR or CC");
          return nullptr;
        }
      return M->clean(epsilon);
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
const RingElement /* or null */ *rawRingElementClean(gmp_RR epsilon,
                                                     const RingElement *f)
{
  const Ring *R = f->get_ring();
  if (R->get_precision() == 0)
    {
      ERROR("expected ring over an RR or CC");
      return nullptr;
    }
  return RingElement::make_raw(R, R->zeroize_tiny(epsilon, f->get_value()));
}
MutableMatrix /* or null */ *rawMutableMatrixClean(gmp_RR epsilon,
                                                   MutableMatrix *M)
{
  /* modifies M in place */
  try
    {
      if (M->get_ring()->get_precision() == 0)
        {
          ERROR("expected ring over an RR or CC");
          return nullptr;
        }
      M->clean(epsilon);
      return M;
  } catch (const exc::engine_error& e)
    {
      ERROR(e.what());
      return nullptr;
  }
}
 
static gmp_RRmutable get_norm_start(gmp_RR p, const Ring *R)
{
  if (R->get_precision() == 0)
    {
      ERROR("expected ring over an RR or CC");
      return nullptr;
    }
  gmp_RRmutable norm = getmemstructtype(gmp_RRmutable);
  mpfr_init2(norm, mpfr_get_prec(p));
  mpfr_ui_div(norm, 1, p, MPFR_RNDN);
  if (!mpfr_zero_p(norm))
    {
      ERROR("Lp norm only implemented for p = infinity");
      mpfr_clear(norm);
      return nullptr;
    }
  return norm;
}
 
gmp_RRorNull rawMatrixNorm(gmp_RR p, const Matrix *M) { return M->norm(p); }
gmp_RRorNull rawRingElementNorm(gmp_RR p, const RingElement *f)
{
  gmp_RRmutable norm = get_norm_start(p, f->get_ring());
  if (!norm) return nullptr;  // error already given.
  f->get_ring()->increase_maxnorm(norm, f->get_value());
  return moveTo_gmpRR(norm);
}
 
gmp_RRorNull rawMutableMatrixNorm(gmp_RR p, const MutableMatrix *M)
{
  (void) p;
  return M->norm();
}

 
// Local Variables:
// compile-command: "make -C $M2BUILDDIR/Macaulay2/e x-mutablemat.o "
// indent-tabs-mode: nil
// End: