[pspp] / src / language / stats / matrix.c

/* PSPP - a program for statistical analysis.
   Copyright (C) 2021 Free Software Foundation, Inc.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>. */

#include <config.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <limits.h>
#include <math.h>
#include <uniwidth.h>

#include "data/any-reader.h"
#include "data/any-writer.h"
#include "data/casereader.h"
#include "data/casewriter.h"
#include "data/data-in.h"
#include "data/data-out.h"
#include "data/dataset.h"
#include "data/dictionary.h"
#include "data/file-handle-def.h"
#include "language/command.h"
#include "language/data-io/data-reader.h"
#include "language/data-io/data-writer.h"
#include "language/data-io/file-handle.h"
#include "language/lexer/format-parser.h"
#include "language/lexer/lexer.h"
#include "language/lexer/variable-parser.h"
#include "libpspp/array.h"
#include "libpspp/assertion.h"
#include "libpspp/compiler.h"
#include "libpspp/hmap.h"
#include "libpspp/i18n.h"
#include "libpspp/misc.h"
#include "libpspp/str.h"
#include "libpspp/string-array.h"
#include "libpspp/stringi-set.h"
#include "libpspp/u8-line.h"
#include "math/random.h"
#include "output/driver.h"
#include "output/output-item.h"
#include "output/pivot-table.h"

#include "gl/c-ctype.h"
#include "gl/minmax.h"
#include "gl/xsize.h"

#include "gettext.h"
#define _(msgid) gettext (msgid)
#define N_(msgid) (msgid)

struct matrix_var
  {
    struct hmap_node hmap_node;
    char *name;
    gsl_matrix *value;
  };

struct msave_common
  {
    /* Configuration. */
    struct file_handle *outfile;
    struct string_array variables;
    struct string_array fnames;
    struct string_array snames;
    bool has_factors;
    bool has_splits;
    size_t n_varnames;

    /* Execution state. */
    struct dictionary *dict;
    struct casewriter *writer;
  };

struct read_file
  {
    struct file_handle *file;
    struct dfm_reader *reader;
    char *encoding;
  };

struct write_file
  {
    struct file_handle *file;
    struct dfm_writer *writer;
    char *encoding;
    struct u8_line *held;
  };

struct matrix_state
  {
    struct dataset *dataset;
    struct session *session;
    struct lexer *lexer;
    struct hmap vars;
    bool in_loop;
    struct file_handle *prev_save_outfile;
    struct msave_common *common;

    struct file_handle *prev_read_file;
    struct read_file **read_files;
    size_t n_read_files;

    struct file_handle *prev_write_file;
    struct write_file **write_files;
    size_t n_write_files;
  };

static struct matrix_var *
matrix_var_lookup (struct matrix_state *s, struct substring name)
{
  struct matrix_var *var;

  HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
                           utf8_hash_case_substring (name, 0), &s->vars)
    if (!utf8_sscasecmp (ss_cstr (var->name), name))
      return var;

  return NULL;
}

static struct matrix_var *
matrix_var_create (struct matrix_state *s, struct substring name)
{
  struct matrix_var *var = xmalloc (sizeof *var);
  *var = (struct matrix_var) { .name = ss_xstrdup (name) };
  hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
  return var;
}

static void
matrix_var_set (struct matrix_var *var, gsl_matrix *value)
{
  gsl_matrix_free (var->value);
  var->value = value;
}
\f
#define MATRIX_FUNCTIONS \
    F(ABS, m_m) \
    F(ALL, d_m) \
    F(ANY, d_m) \
    F(ARSIN, m_m) \
    F(ARTAN, m_m) \
    F(BLOCK, m_any) \
    F(CHOL, m_m) \
    F(CMIN, m_m) \
    F(CMAX, m_m) \
    F(COS, m_m) \
    F(CSSQ, m_m) \
    F(CSUM, m_m) \
    F(DESIGN, m_m) \
    F(DET, d_m) \
    F(DIAG, m_m) \
    F(EVAL, m_m) \
    F(EXP, m_m) \
    F(GINV, m_m) \
    F(GRADE, m_m) \
    F(GSCH, m_m) \
    F(IDENT, IDENT) \
    F(INV, m_m) \
    F(KRONEKER, m_mm) \
    F(LG10, m_m) \
    F(LN, m_m) \
    F(MAGIC, m_d) \
    F(MAKE, m_ddd) \
    F(MDIAG, m_v) \
    F(MMAX, d_m) \
    F(MMIN, d_m) \
    F(MOD, m_md) \
    F(MSSQ, d_m) \
    F(MSUM, d_m) \
    F(NCOL, d_m) \
    F(NROW, d_m) \
    F(RANK, d_m) \
    F(RESHAPE, m_mdd) \
    F(RMAX, m_m) \
    F(RMIN, m_m) \
    F(RND, m_m) \
    F(RNKORDER, m_m) \
    F(RSSQ, m_m) \
    F(RSUM, m_m) \
    F(SIN, m_m) \
    F(SOLVE, m_mm) \
    F(SQRT, m_m) \
    F(SSCP, m_m) \
    F(SVAL, m_m) \
    F(SWEEP, m_md) \
    F(T, m_m) \
    F(TRACE, d_m) \
    F(TRANSPOS, m_m) \
    F(TRUNC, m_m) \
    F(UNIFORM, m_dd)

enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };

typedef gsl_matrix *matrix_proto_m_d (double);
typedef gsl_matrix *matrix_proto_m_dd (double, double);
typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
typedef double matrix_proto_d_m (gsl_matrix *);
typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
typedef gsl_matrix *matrix_proto_IDENT (double, double);

#define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
MATRIX_FUNCTIONS
#undef F
\f
/* Matrix expressions. */

struct matrix_expr
  {
    enum matrix_op
      {
        /* Functions. */
#define F(NAME, PROTOTYPE) MOP_F_##NAME,
        MATRIX_FUNCTIONS
#undef F

        /* Elementwise and scalar arithmetic. */
        MOP_NEGATE,             /* unary - */
        MOP_ADD_ELEMS,          /* + */
        MOP_SUB_ELEMS,          /* - */
        MOP_MUL_ELEMS,          /* &* */
        MOP_DIV_ELEMS,          /* / and &/ */
        MOP_EXP_ELEMS,          /* &** */
        MOP_SEQ,                /* a:b */
        MOP_SEQ_BY,             /* a:b:c */

        /* Matrix arithmetic. */
        MOP_MUL_MAT,            /* * */
        MOP_EXP_MAT,            /* ** */

        /* Relational. */
        MOP_GT,                 /* > */
        MOP_GE,                 /* >= */
        MOP_LT,                 /* < */
        MOP_LE,                 /* <= */
        MOP_EQ,                 /* = */
        MOP_NE,                 /* <> */

        /* Logical. */
        MOP_NOT,                /* NOT */
        MOP_AND,                /* AND */
        MOP_OR,                 /* OR */
        MOP_XOR,                /* XOR */

        /* {}. */
        MOP_PASTE_HORZ,         /* a, b, c, ... */
        MOP_PASTE_VERT,         /* a; b; c; ... */
        MOP_EMPTY,              /* {} */

        /* Sub-matrices. */
        MOP_VEC_INDEX,          /* x(y) */
        MOP_VEC_ALL,            /* x(:) */
        MOP_MAT_INDEX,          /* x(y,z) */
        MOP_ROW_INDEX,          /* x(y,:) */
        MOP_COL_INDEX,          /* x(:,z) */

        /* Literals. */
        MOP_NUMBER,
        MOP_VARIABLE,

        /* Oddball stuff. */
        MOP_EOF,                /* EOF('file') */
      }
    op;

    union
      {
        struct
          {
            struct matrix_expr **subs;
            size_t n_subs;
          };

        double number;
        struct matrix_var *variable;
        struct read_file *eof;
      };
  };

static void
matrix_expr_destroy (struct matrix_expr *e)
{
  if (!e)
    return;

  switch (e->op)
    {
#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
MATRIX_FUNCTIONS
#undef F
    case MOP_NEGATE:
    case MOP_ADD_ELEMS:
    case MOP_SUB_ELEMS:
    case MOP_MUL_ELEMS:
    case MOP_DIV_ELEMS:
    case MOP_EXP_ELEMS:
    case MOP_SEQ:
    case MOP_SEQ_BY:
    case MOP_MUL_MAT:
    case MOP_EXP_MAT:
    case MOP_GT:
    case MOP_GE:
    case MOP_LT:
    case MOP_LE:
    case MOP_EQ:
    case MOP_NE:
    case MOP_NOT:
    case MOP_AND:
    case MOP_OR:
    case MOP_XOR:
    case MOP_EMPTY:
    case MOP_PASTE_HORZ:
    case MOP_PASTE_VERT:
    case MOP_VEC_INDEX:
    case MOP_VEC_ALL:
    case MOP_MAT_INDEX:
    case MOP_ROW_INDEX:
    case MOP_COL_INDEX:
      for (size_t i = 0; i < e->n_subs; i++)
        matrix_expr_destroy (e->subs[i]);
      free (e->subs);
      break;

    case MOP_NUMBER:
    case MOP_VARIABLE:
    case MOP_EOF:
      break;
    }
  free (e);
}

static struct matrix_expr *
matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
                         size_t n_subs)
{
  struct matrix_expr *e = xmalloc (sizeof *e);
  *e = (struct matrix_expr) {
    .op = op,
    .subs = xmemdup (subs, n_subs * sizeof *subs),
    .n_subs = n_subs
  };
  return e;
}

static struct matrix_expr *
matrix_expr_create_0 (enum matrix_op op)
{
  struct matrix_expr *sub;
  return matrix_expr_create_subs (op, &sub, 0);
}

static struct matrix_expr *
matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
{
  return matrix_expr_create_subs (op, &sub, 1);
}

static struct matrix_expr *
matrix_expr_create_2 (enum matrix_op op,
                      struct matrix_expr *sub0, struct matrix_expr *sub1)
{
  struct matrix_expr *subs[] = { sub0, sub1 };
  return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
}

static struct matrix_expr *
matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
                      struct matrix_expr *sub1, struct matrix_expr *sub2)
{
  struct matrix_expr *subs[] = { sub0, sub1, sub2 };
  return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
}

static struct matrix_expr *
matrix_expr_create_number (double number)
{
  struct matrix_expr *e = xmalloc (sizeof *e);
  *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
  return e;
}

static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);

static struct matrix_expr *
matrix_parse_curly_comma (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_or_xor (s);
  if (!lhs)
    return NULL;

  while (lex_match (s->lexer, T_COMMA))
    {
      struct matrix_expr *rhs = matrix_parse_or_xor (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
    }
  return lhs;
}

static struct matrix_expr *
matrix_parse_curly_semi (struct matrix_state *s)
{
  if (lex_token (s->lexer) == T_RCURLY)
    return matrix_expr_create_0 (MOP_EMPTY);

  struct matrix_expr *lhs = matrix_parse_curly_comma (s);
  if (!lhs)
    return NULL;

  while (lex_match (s->lexer, T_SEMICOLON))
    {
      struct matrix_expr *rhs = matrix_parse_curly_comma (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
    }
  return lhs;
}

#define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M)                     \
  for (size_t Y = 0; Y < (M)->size1; Y++)                       \
    for (size_t X = 0; X < (M)->size2; X++)                     \
      for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)

static bool
is_vector (const gsl_matrix *m)
{
  return m->size1 <= 1 || m->size2 <= 1;
}

static gsl_vector
to_vector (gsl_matrix *m)
{
  return (m->size1 == 1
          ? gsl_matrix_row (m, 0).vector
          : gsl_matrix_column (m, 0).vector);
}


static gsl_matrix *
matrix_eval_ABS (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = fabs (*d);
  return m;
}

static double
matrix_eval_ALL (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    if (*d == 0.0)
      return 0.0;
  return 1.0;
}

static double
matrix_eval_ANY (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    if (*d != 0.0)
      return !.0;
  return 0.0;
}

static gsl_matrix *
matrix_eval_ARSIN (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = asin (*d);
  return m;
}

static gsl_matrix *
matrix_eval_ARTAN (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = atan (*d);
  return m;
}

static gsl_matrix *
matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
{
  size_t r = 0;
  size_t c = 0;
  for (size_t i = 0; i < n; i++)
    {
      r += m[i]->size1;
      c += m[i]->size2;
    }
  gsl_matrix *block = gsl_matrix_calloc (r, c);
  r = c = 0;
  for (size_t i = 0; i < n; i++)
    {
      for (size_t y = 0; y < m[i]->size1; y++)
        for (size_t x = 0; x < m[i]->size2; x++)
          gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
      r += m[i]->size1;
      c += m[i]->size2;
    }
  return block;
}

static gsl_matrix *
matrix_eval_CHOL (gsl_matrix *m)
{
  if (!gsl_linalg_cholesky_decomp1 (m))
    {
      for (size_t y = 0; y < m->size1; y++)
        for (size_t x = y + 1; x < m->size2; x++)
          gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));

      for (size_t y = 0; y < m->size1; y++)
        for (size_t x = 0; x < y; x++)
          gsl_matrix_set (m, y, x, 0);
      return m;
    }
  else
    {
      msg (SE, _("Input to CHOL function is not positive-definite."));
      return NULL;
    }
}

static gsl_matrix *
matrix_eval_col_extremum (gsl_matrix *m, bool min)
{
  if (m->size1 <= 1)
    return m;
  else if (!m->size2)
    return gsl_matrix_alloc (1, 0);

  gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
  for (size_t x = 0; x < m->size2; x++)
    {
      double ext = gsl_matrix_get (m, 0, x);
      for (size_t y = 1; y < m->size1; y++)
        {
          double value = gsl_matrix_get (m, y, x);
          if (min ? value < ext : value > ext)
            ext = value;
        }
      gsl_matrix_set (cext, 0, x, ext);
    }
  return cext;
}

static gsl_matrix *
matrix_eval_CMAX (gsl_matrix *m)
{
  return matrix_eval_col_extremum (m, false);
}

static gsl_matrix *
matrix_eval_CMIN (gsl_matrix *m)
{
  return matrix_eval_col_extremum (m, true);
}

static gsl_matrix *
matrix_eval_COS (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = cos (*d);
  return m;
}

static gsl_matrix *
matrix_eval_col_sum (gsl_matrix *m, bool square)
{
  if (m->size1 == 0)
    return m;
  else if (!m->size2)
    return gsl_matrix_alloc (1, 0);

  gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
  for (size_t x = 0; x < m->size2; x++)
    {
      double sum = 0;
      for (size_t y = 0; y < m->size1; y++)
        {
          double d = gsl_matrix_get (m, y, x);
          sum += square ? pow2 (d) : d;
        }
      gsl_matrix_set (result, 0, x, sum);
    }
  return result;
}

static gsl_matrix *
matrix_eval_CSSQ (gsl_matrix *m)
{
  return matrix_eval_col_sum (m, true);
}

static gsl_matrix *
matrix_eval_CSUM (gsl_matrix *m)
{
  return matrix_eval_col_sum (m, false);
}

static int
compare_double_3way (const void *a_, const void *b_)
{
  const double *a = a_;
  const double *b = b_;
  return *a < *b ? -1 : *a > *b;
}

static gsl_matrix *
matrix_eval_DESIGN (gsl_matrix *m)
{
  double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
  gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
  gsl_matrix_transpose_memcpy (&m2, m);

  for (size_t y = 0; y < m2.size1; y++)
    qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);

  size_t *n = xcalloc (m2.size1, sizeof *n);
  size_t n_total = 0;
  for (size_t i = 0; i < m2.size1; i++)
    {
      double *row = tmp + m2.size2 * i;
      for (size_t j = 0; j < m2.size2; )
        {
          size_t k;
          for (k = j + 1; k < m2.size2; k++)
            if (row[j] != row[k])
              break;
          row[n[i]++] = row[j];
          j = k;
        }

      if (n[i] <= 1)
        msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
      else
        n_total += n[i];
    }

  gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
  size_t x = 0;
  for (size_t i = 0; i < m->size2; i++)
    {
      if (n[i] <= 1)
        continue;

      const double *unique = tmp + m2.size2 * i;
      for (size_t j = 0; j < n[i]; j++, x++)
        {
          double value = unique[j];
          for (size_t y = 0; y < m->size1; y++)
            gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
        }
    }

  free (n);
  free (tmp);

  return result;
}

static double
matrix_eval_DET (gsl_matrix *m)
{
  gsl_permutation *p = gsl_permutation_alloc (m->size1);
  int signum;
  gsl_linalg_LU_decomp (m, p, &signum);
  gsl_permutation_free (p);
  return gsl_linalg_LU_det (m, signum);
}

static gsl_matrix *
matrix_eval_DIAG (gsl_matrix *m)
{
  gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
  for (size_t i = 0; i < diag->size1; i++)
    gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
  return diag;
}

static bool
is_symmetric (const gsl_matrix *m)
{
  if (m->size1 != m->size2)
    return false;

  for (size_t y = 0; y < m->size1; y++)
    for (size_t x = 0; x < y; x++)
      if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
        return false;

  return true;
}

static int
compare_double_desc (const void *a_, const void *b_)
{
  const double *a = a_;
  const double *b = b_;
  return *a > *b ? -1 : *a < *b;
}

static gsl_matrix *
matrix_eval_EVAL (gsl_matrix *m)
{
  if (!is_symmetric (m))
    {
      msg (SE, _("Argument of EVAL must be symmetric."));
      return NULL;
    }

  gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
  gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
  gsl_vector v_eval = to_vector (eval);
  gsl_eigen_symm (m, &v_eval, w);
  gsl_eigen_symm_free (w);

  assert (v_eval.stride == 1);
  qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);

  return eval;
}

static gsl_matrix *
matrix_eval_EXP (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = exp (*d);
  return m;
}

/* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
   marked as:

   Charl Linssen <charl@itfromb.it>
   Feb 2016
   PUBLIC DOMAIN */
static gsl_matrix *
matrix_eval_GINV (gsl_matrix *A)
{
  size_t n = A->size1;
  size_t m = A->size2;
  bool swap = m > n;
  gsl_matrix *tmp_mat = NULL;
  if (swap)
    {
      /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
      tmp_mat = gsl_matrix_alloc (m, n);
      gsl_matrix_transpose_memcpy (tmp_mat, A);
      A = tmp_mat;
      size_t i = m;
      m = n;
      n = i;
    }

  /* Do SVD. */
  gsl_matrix *V = gsl_matrix_alloc (m, m);
  gsl_vector *u = gsl_vector_alloc (m);

  gsl_vector *tmp_vec = gsl_vector_alloc (m);
  gsl_linalg_SV_decomp (A, V, u, tmp_vec);
  gsl_vector_free (tmp_vec);

  /* Compute Σ⁻¹. */
  gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
  gsl_matrix_set_zero (Sigma_pinv);
  double cutoff = 1e-15 * gsl_vector_max (u);

  for (size_t i = 0; i < m; ++i)
    {
      double x = gsl_vector_get (u, i);
      gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
    }

  /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
  gsl_matrix *U = gsl_matrix_calloc (n, n);
  for (size_t i = 0; i < n; i++)
    for (size_t j = 0; j < m; j++)
      gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));

  /* two dot products to obtain pseudoinverse */
  gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);

  gsl_matrix *A_pinv;
  if (swap)
    {
      A_pinv = gsl_matrix_alloc (n, m);
      gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
    }
  else
    {
      A_pinv = gsl_matrix_alloc (m, n);
      gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
    }

  gsl_matrix_free (tmp_mat);
  gsl_matrix_free (tmp_mat2);
  gsl_matrix_free (U);
  gsl_matrix_free (Sigma_pinv);
  gsl_vector_free (u);
  gsl_matrix_free (V);

  return A_pinv;
}

struct grade
  {
    size_t y, x;
    double value;
  };

static int
grade_compare_3way (const void *a_, const void *b_)
{
  const struct grade *a = a_;
  const struct grade *b = b_;

  return (a->value < b->value ? -1
          : a->value > b->value ? 1
          : a->y < b->y ? -1
          : a->y > b->y ? 1
          : a->x < b->x ? -1
          : a->x > b->x);
}

static gsl_matrix *
matrix_eval_GRADE (gsl_matrix *m)
{
  size_t n = m->size1 * m->size2;
  struct grade *grades = xmalloc (n * sizeof *grades);

  size_t i = 0;
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
  qsort (grades, n, sizeof *grades, grade_compare_3way);

  for (size_t i = 0; i < n; i++)
    gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);

  free (grades);

  return m;
}

static double
dot (gsl_vector *a, gsl_vector *b)
{
  double result = 0.0;
  for (size_t i = 0; i < a->size; i++)
    result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
  return result;
}

static double
norm2 (gsl_vector *v)
{
  double result = 0.0;
  for (size_t i = 0; i < v->size; i++)
    result += pow2 (gsl_vector_get (v, i));
  return result;
}

static double
norm (gsl_vector *v)
{
  return sqrt (norm2 (v));
}

static gsl_matrix *
matrix_eval_GSCH (gsl_matrix *v)
{
  if (v->size2 < v->size1)
    {
      msg (SE, _("GSCH requires its argument to have at least as many columns "
                 "as rows, but it has dimensions %zu×%zu."),
           v->size1, v->size2);
      return NULL;
    }
  if (!v->size1 || !v->size2)
    return v;

  gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
  size_t ux = 0;
  for (size_t vx = 0; vx < v->size2; vx++)
    {
      gsl_vector u_i = gsl_matrix_column (u, ux).vector;
      gsl_vector v_i = gsl_matrix_column (v, vx).vector;

      gsl_vector_memcpy (&u_i, &v_i);
      for (size_t j = 0; j < ux; j++)
        {
          gsl_vector u_j = gsl_matrix_column (u, j).vector;
          double scale = dot (&u_j, &u_i) / norm2 (&u_j);
          for (size_t k = 0; k < u_i.size; k++)
            gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
                                      - scale * gsl_vector_get (&u_j, k)));
        }

      double len = norm (&u_i);
      if (len > 1e-15)
        {
          gsl_vector_scale (&u_i, 1.0 / len);
          if (++ux >= v->size1)
            break;
        }
    }

  if (ux < v->size1)
    {
      msg (SE, _("%zu×%zu argument to GSCH contains only "
                 "%zu linearly independent columns."),
           v->size1, v->size2, ux);
      gsl_matrix_free (u);
      return NULL;
    }

  u->size2 = v->size1;
  return u;
}

static gsl_matrix *
matrix_eval_IDENT (double s1, double s2)
{
  if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
    {
      msg (SE, _("Arguments to IDENT must be non-negative integers."));
      return NULL;
    }

  gsl_matrix *m = gsl_matrix_alloc (s1, s2);
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = x == y;
  return m;
}

static void
invert_matrix (gsl_matrix *x)
{
  gsl_permutation *p = gsl_permutation_alloc (x->size1);
  int signum;
  gsl_linalg_LU_decomp (x, p, &signum);
  gsl_linalg_LU_invx (x, p);
  gsl_permutation_free (p);
}

static gsl_matrix *
matrix_eval_INV (gsl_matrix *m)
{
  invert_matrix (m);
  return m;
}

static gsl_matrix *
matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
{
  gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
                                    a->size2 * b->size2);
  size_t y = 0;
  for (size_t ar = 0; ar < a->size1; ar++)
    for (size_t br = 0; br < b->size1; br++, y++)
      {
        size_t x = 0;
        for (size_t ac = 0; ac < a->size2; ac++)
          for (size_t bc = 0; bc < b->size2; bc++, x++)
            {
              double av = gsl_matrix_get (a, ar, ac);
              double bv = gsl_matrix_get (b, br, bc);
              gsl_matrix_set (k, y, x, av * bv);
            }
      }
  return k;
}

static gsl_matrix *
matrix_eval_LG10 (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = log10 (*d);
  return m;
}

static gsl_matrix *
matrix_eval_LN (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = log (*d);
  return m;
}

static void
matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
{
  /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
  size_t y = 0;
  size_t x = n / 2;
  for (size_t i = 1; i <= n * n; i++)
    {
      gsl_matrix_set (m, y, x, i);

      size_t y1 = !y ? n - 1 : y - 1;
      size_t x1 = x + 1 >= n ? 0 : x + 1;
      if (gsl_matrix_get (m, y1, x1) == 0)
        {
          y = y1;
          x = x1;
        }
      else
        y = y + 1 >= n ? 0 : y + 1;
    }
}

static void
magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
{
  double a = gsl_matrix_get (m, y1, x1);
  double b = gsl_matrix_get (m, y2, x2);
  gsl_matrix_set (m, y1, x1, b);
  gsl_matrix_set (m, y2, x2, a);
}

static void
matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
{
  size_t x, y;

  /* A. Umar, "On the Construction of Even Order Magic Squares",
     https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
  x = y = 0;
  for (size_t i = 1; i <= n * n / 2; i++)
    {
      gsl_matrix_set (m, y, x, i);
      if (++y >= n)
        {
          y = 0;
          x++;
        }
    }

  x = n - 1;
  y = 0;
  for (size_t i = n * n; i > n * n / 2; i--)
    {
      gsl_matrix_set (m, y, x, i);
      if (++y >= n)
        {
          y = 0;
          x--;
        }
    }

  for (size_t y = 0; y < n; y++)
    for (size_t x = 0; x < n / 2; x++)
      {
        unsigned int d = gsl_matrix_get (m, y, x);
        if (d % 2 != (y < n / 2))
          magic_exchange (m, y, x, y, n - x - 1);
      }

  size_t y1 = n / 2;
  size_t y2 = n - 1;
  size_t x1 = n / 2 - 1;
  size_t x2 = n / 2;
  magic_exchange (m, y1, x1, y2, x1);
  magic_exchange (m, y1, x2, y2, x2);
}

static void
matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
{
  /* A. Umar, "On the Construction of Even Order Magic Squares",
     https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
  size_t x, y;

  x = y = 0;
  for (size_t i = 1; ; i++)
    {
      gsl_matrix_set (m, y, x, i);
      if (++y == n / 2 - 1)
        y += 2;
      else if (y >= n)
        {
          y = 0;
          if (++x >= n / 2)
            break;
        }
    }

  x = n - 1;
  y = 0;
  for (size_t i = n * n; ; i--)
    {
      gsl_matrix_set (m, y, x, i);
      if (++y == n / 2 - 1)
        y += 2;
      else if (y >= n)
        {
          y = 0;
          if (--x < n / 2)
            break;
        }
    }
  for (size_t y = 0; y < n; y++)
    if (y != n / 2 - 1 && y != n / 2)
      for (size_t x = 0; x < n / 2; x++)
        {
          unsigned int d = gsl_matrix_get (m, y, x);
          if (d % 2 != (y < n / 2))
            magic_exchange (m, y, x, y, n - x - 1);
        }

  size_t a0 = (n * n - 2 * n) / 2 + 1;
  for (size_t i = 0; i < n / 2; i++)
    {
      size_t a = a0 + i;
      gsl_matrix_set (m, n / 2, i, a);
      gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
    }
  for (size_t i = 0; i < n / 2; i++)
    {
      size_t a = a0 + i + n / 2;
      gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
      gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
    }
  for (size_t x = 1; x < n / 2; x += 2)
    magic_exchange (m, n / 2, x, n / 2 - 1, x);
  for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
    magic_exchange (m, n / 2, x, n / 2 - 1, x);
  size_t x1 = n / 2 - 2;
  size_t x2 = n / 2 + 1;
  size_t y1 = n / 2 - 2;
  size_t y2 = n / 2 + 1;
  magic_exchange (m, y1, x1, y2, x1);
  magic_exchange (m, y1, x2, y2, x2);
}

static gsl_matrix *
matrix_eval_MAGIC (double n_)
{
  if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
    {
      msg (SE, _("MAGIC argument must be an integer 3 or greater."));
      return NULL;
    }
  size_t n = n_;

  gsl_matrix *m = gsl_matrix_calloc (n, n);
  if (n % 2)
    matrix_eval_MAGIC_odd (m, n);
  else if (n % 4)
    matrix_eval_MAGIC_singly_even (m, n);
  else
    matrix_eval_MAGIC_doubly_even (m, n);
  return m;
}

static gsl_matrix *
matrix_eval_MAKE (double r, double c, double value)
{
  if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
    {
      msg (SE, _("Size arguments to MAKE must be integers."));
      return NULL;
    }

  gsl_matrix *m = gsl_matrix_alloc (r, c);
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = value;
  return m;
}

static gsl_matrix *
matrix_eval_MDIAG (gsl_vector *v)
{
  gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
  gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
  gsl_vector_memcpy (&diagonal, v);
  return m;
}

static double
matrix_eval_MMAX (gsl_matrix *m)
{
  return gsl_matrix_max (m);
}

static double
matrix_eval_MMIN (gsl_matrix *m)
{
  return gsl_matrix_min (m);
}

static gsl_matrix *
matrix_eval_MOD (gsl_matrix *m, double divisor)
{
  if (divisor == 0.0)
    {
      msg (SE, _("Divisor argument to MOD function must be nonzero."));
      return NULL;
    }

  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = fmod (*d, divisor);
  return m;
}

static double
matrix_eval_MSSQ (gsl_matrix *m)
{
  double mssq = 0.0;
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    mssq += *d * *d;
  return mssq;
}

static double
matrix_eval_MSUM (gsl_matrix *m)
{
  double msum = 0.0;
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    msum += *d;
  return msum;
}

static double
matrix_eval_NCOL (gsl_matrix *m)
{
  return m->size2;
}

static double
matrix_eval_NROW (gsl_matrix *m)
{
  return m->size1;
}

static double
matrix_eval_RANK (gsl_matrix *m)
{
  gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
  gsl_linalg_QR_decomp (m, tau);
  gsl_vector_free (tau);

  return gsl_linalg_QRPT_rank (m, -1);
}

static gsl_matrix *
matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
{
  if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
    {
      msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
      return NULL;
    }
  size_t r = r_;
  size_t c = c_;
  if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
    {
      msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
                 "product of matrix dimensions."));
      return NULL;
    }

  gsl_matrix *dst = gsl_matrix_alloc (r, c);
  size_t y1 = 0;
  size_t x1 = 0;
  MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
    {
      gsl_matrix_set (dst, y1, x1, *d);
      if (++x1 >= c)
        {
          x1 = 0;
          y1++;
        }
    }
  return dst;
}

static gsl_matrix *
matrix_eval_row_extremum (gsl_matrix *m, bool min)
{
  if (m->size2 <= 1)
    return m;
  else if (!m->size1)
    return gsl_matrix_alloc (0, 1);

  gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
  for (size_t y = 0; y < m->size1; y++)
    {
      double ext = gsl_matrix_get (m, y, 0);
      for (size_t x = 1; x < m->size2; x++)
        {
          double value = gsl_matrix_get (m, y, x);
          if (min ? value < ext : value > ext)
            ext = value;
        }
      gsl_matrix_set (rext, y, 0, ext);
    }
  return rext;
}

static gsl_matrix *
matrix_eval_RMAX (gsl_matrix *m)
{
  return matrix_eval_row_extremum (m, false);
}

static gsl_matrix *
matrix_eval_RMIN (gsl_matrix *m)
{
  return matrix_eval_row_extremum (m, true);
}

static gsl_matrix *
matrix_eval_RND (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = rint (*d);
  return m;
}

struct rank
  {
    size_t y, x;
    double value;
  };

static int
rank_compare_3way (const void *a_, const void *b_)
{
  const struct rank *a = a_;
  const struct rank *b = b_;

  return a->value < b->value ? -1 : a->value > b->value;
}

static gsl_matrix *
matrix_eval_RNKORDER (gsl_matrix *m)
{
  size_t n = m->size1 * m->size2;
  struct rank *ranks = xmalloc (n * sizeof *ranks);
  size_t i = 0;
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
  qsort (ranks, n, sizeof *ranks, rank_compare_3way);

  for (size_t i = 0; i < n; )
    {
      size_t j;
      for (j = i + 1; j < n; j++)
        if (ranks[i].value != ranks[j].value)
          break;

      double rank = (i + j + 1.0) / 2.0;
      for (size_t k = i; k < j; k++)
        gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);

      i = j;
    }

  free (ranks);

  return m;
}

static gsl_matrix *
matrix_eval_row_sum (gsl_matrix *m, bool square)
{
  if (m->size1 == 0)
    return m;
  else if (!m->size1)
    return gsl_matrix_alloc (0, 1);

  gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
  for (size_t y = 0; y < m->size1; y++)
    {
      double sum = 0;
      for (size_t x = 0; x < m->size2; x++)
        {
          double d = gsl_matrix_get (m, y, x);
          sum += square ? pow2 (d) : d;
        }
      gsl_matrix_set (result, y, 0, sum);
    }
  return result;
}

static gsl_matrix *
matrix_eval_RSSQ (gsl_matrix *m)
{
  return matrix_eval_row_sum (m, true);
}

static gsl_matrix *
matrix_eval_RSUM (gsl_matrix *m)
{
  return matrix_eval_row_sum (m, false);
}

static gsl_matrix *
matrix_eval_SIN (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = sin (*d);
  return m;
}

static gsl_matrix *
matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
{
  if (m1->size1 != m2->size1)
    {
      msg (SE, _("SOLVE requires its arguments to have the same number of "
                 "rows, but the first argument has dimensions %zu×%zu and "
                 "the second %zu×%zu."),
           m1->size1, m1->size2,
           m2->size1, m2->size2);
      return NULL;
    }

  gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
  gsl_permutation *p = gsl_permutation_alloc (m1->size1);
  int signum;
  gsl_linalg_LU_decomp (m1, p, &signum);
  for (size_t i = 0; i < m2->size2; i++)
    {
      gsl_vector bi = gsl_matrix_column (m2, i).vector;
      gsl_vector xi = gsl_matrix_column (x, i).vector;
      gsl_linalg_LU_solve (m1, p, &bi, &xi);
    }
  gsl_permutation_free (p);
  return x;
}

static gsl_matrix *
matrix_eval_SQRT (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    {
      if (*d < 0)
        {
          msg (SE, _("Argument to SQRT must be nonnegative."));
          return NULL;
        }
      *d = sqrt (*d);
    }
  return m;
}

static gsl_matrix *
matrix_eval_SSCP (gsl_matrix *m)
{
  gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
  gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
  return sscp;
}

static gsl_matrix *
matrix_eval_SVAL (gsl_matrix *m)
{
  gsl_matrix *tmp_mat = NULL;
  if (m->size2 > m->size1)
    {
      tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
      gsl_matrix_transpose_memcpy (tmp_mat, m);
      m = tmp_mat;
    }

  /* Do SVD. */
  gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
  gsl_vector *S = gsl_vector_alloc (m->size2);
  gsl_vector *work = gsl_vector_alloc (m->size2);
  gsl_linalg_SV_decomp (m, V, S, work);

  gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
  for (size_t i = 0; i < m->size2; i++)
    gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));

  gsl_matrix_free (V);
  gsl_vector_free (S);
  gsl_vector_free (work);
  gsl_matrix_free (tmp_mat);

  return vals;
}

static gsl_matrix *
matrix_eval_SWEEP (gsl_matrix *m, double d)
{
  if (d < 1 || d > SIZE_MAX)
    {
      msg (SE, _("Scalar argument to SWEEP must be integer."));
      return NULL;
    }
  size_t k = d - 1;
  if (k >= MIN (m->size1, m->size2))
    {
      msg (SE, _("Scalar argument to SWEEP must be integer less than or "
                 "equal to the smaller of the matrix argument's rows and "
                 "columns."));
      return NULL;
    }

  double m_kk = gsl_matrix_get (m, k, k);
  if (fabs (m_kk) > 1e-19)
    {
      gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
      MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
        {
          double m_ij = gsl_matrix_get (m, i, j);
          double m_ik = gsl_matrix_get (m, i, k);
          double m_kj = gsl_matrix_get (m, k, j);
          *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
                   : i != k ? -m_ik
                   : j != k ? m_kj
                   : 1.0) / m_kk;
        }
      return a;
    }
  else
    {
      for (size_t i = 0; i < m->size1; i++)
        {
          gsl_matrix_set (m, i, k, 0);
          gsl_matrix_set (m, k, i, 0);
        }
      return m;
    }
}

static double
matrix_eval_TRACE (gsl_matrix *m)
{
  double sum = 0;
  size_t n = MIN (m->size1, m->size2);
  for (size_t i = 0; i < n; i++)
    sum += gsl_matrix_get (m, i, i);
  return sum;
}

static gsl_matrix *
matrix_eval_T (gsl_matrix *m)
{
  return matrix_eval_TRANSPOS (m);
}

static gsl_matrix *
matrix_eval_TRANSPOS (gsl_matrix *m)
{
  if (m->size1 == m->size2)
    {
      gsl_matrix_transpose (m);
      return m;
    }
  else
    {
      gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
      gsl_matrix_transpose_memcpy (t, m);
      return t;
    }
}

static gsl_matrix *
matrix_eval_TRUNC (gsl_matrix *m)
{
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = trunc (*d);
  return m;
}

static gsl_matrix *
matrix_eval_UNIFORM (double r_, double c_)
{
  if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
    {
      msg (SE, _("Arguments to UNIFORM must be integers."));
      return NULL;
    }
  size_t r = r_;
  size_t c = c_;
  if (size_overflow_p (xtimes (r, xmax (c, 1))))
    {
      msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
      return NULL;
    }

  gsl_matrix *m = gsl_matrix_alloc (r, c);
  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
    *d = gsl_ran_flat (get_rng (), 0, 1);
  return m;
}

struct matrix_function
  {
    const char *name;
    enum matrix_op op;
    size_t min_args, max_args;
  };

static struct matrix_expr *matrix_parse_expr (struct matrix_state *);

static const struct matrix_function *
matrix_parse_function_name (const char *token)
{
  static const struct matrix_function functions[] =
    {
#define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
      MATRIX_FUNCTIONS
#undef F
    };
  enum { N_FUNCTIONS = sizeof functions / sizeof *functions };

  for (size_t i = 0; i < N_FUNCTIONS; i++)
    {
      if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
        return &functions[i];
    }
  return NULL;
}

static struct read_file *
read_file_create (struct matrix_state *s, struct file_handle *fh)
{
  for (size_t i = 0; i < s->n_read_files; i++)
    {
      struct read_file *rf = s->read_files[i];
      if (rf->file == fh)
        {
          fh_unref (fh);
          return rf;
        }
    }

  struct read_file *rf = xmalloc (sizeof *rf);
  *rf = (struct read_file) { .file = fh };

  s->read_files = xrealloc (s->read_files,
                            (s->n_read_files + 1) * sizeof *s->read_files);
  s->read_files[s->n_read_files++] = rf;
  return rf;
}

static struct dfm_reader *
read_file_open (struct read_file *rf)
{
  if (!rf->reader)
    rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
  return rf->reader;
}

static void
read_file_destroy (struct read_file *rf)
{
  if (rf)
    {
      fh_unref (rf->file);
      dfm_close_reader (rf->reader);
      free (rf->encoding);
      free (rf);
    }
}

static struct write_file *
write_file_create (struct matrix_state *s, struct file_handle *fh)
{
  for (size_t i = 0; i < s->n_write_files; i++)
    {
      struct write_file *wf = s->write_files[i];
      if (wf->file == fh)
        {
          fh_unref (fh);
          return wf;
        }
    }

  struct write_file *wf = xmalloc (sizeof *wf);
  *wf = (struct write_file) { .file = fh };

  s->write_files = xrealloc (s->write_files,
                             (s->n_write_files + 1) * sizeof *s->write_files);
  s->write_files[s->n_write_files++] = wf;
  return wf;
}

static struct dfm_writer *
write_file_open (struct write_file *wf)
{
  if (!wf->writer)
    wf->writer = dfm_open_writer (wf->file, wf->encoding);
  return wf->writer;
}

static void
write_file_destroy (struct write_file *wf)
{
  if (wf)
    {
      if (wf->held)
        {
          dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
                               wf->held->s.ss.length);
          u8_line_destroy (wf->held);
          free (wf->held);
        }

      fh_unref (wf->file);
      dfm_close_writer (wf->writer);
      free (wf->encoding);
      free (wf);
    }
}

static bool
matrix_parse_function (struct matrix_state *s, const char *token,
                       struct matrix_expr **exprp)
{
  *exprp = NULL;
  if (lex_next_token (s->lexer, 1) != T_LPAREN)
    return false;

  if (lex_match_id (s->lexer, "EOF"))
    {
      lex_get (s->lexer);
      struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
      if (!fh)
        return true;

      if (!lex_force_match (s->lexer, T_RPAREN))
        {
          fh_unref (fh);
          return true;
        }

      struct read_file *rf = read_file_create (s, fh);

      struct matrix_expr *e = xmalloc (sizeof *e);
      *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
      *exprp = e;
      return true;
    }

  const struct matrix_function *f = matrix_parse_function_name (token);
  if (!f)
    return false;

  lex_get_n (s->lexer, 2);

  struct matrix_expr *e = xmalloc (sizeof *e);
  *e = (struct matrix_expr) { .op = f->op, .subs = NULL };

  size_t allocated_subs = 0;
  do
    {
      struct matrix_expr *sub = matrix_parse_expr (s);
      if (!sub)
        goto error;

      if (e->n_subs >= allocated_subs)
        e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
      e->subs[e->n_subs++] = sub;
    }
  while (lex_match (s->lexer, T_COMMA));
  if (!lex_force_match (s->lexer, T_RPAREN))
    goto error;

  *exprp = e;
  return true;

error:
  matrix_expr_destroy (e);
  return true;
}

static struct matrix_expr *
matrix_parse_primary (struct matrix_state *s)
{
  if (lex_is_number (s->lexer))
    {
      double number = lex_number (s->lexer);
      lex_get (s->lexer);
      return matrix_expr_create_number (number);
    }
  else if (lex_is_string (s->lexer))
    {
      char string[sizeof (double)];
      buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
      lex_get (s->lexer);

      double number;
      memcpy (&number, string, sizeof number);
      return matrix_expr_create_number (number);
    }
  else if (lex_match (s->lexer, T_LPAREN))
    {
      struct matrix_expr *e = matrix_parse_or_xor (s);
      if (!e || !lex_force_match (s->lexer, T_RPAREN))
        {
          matrix_expr_destroy (e);
          return NULL;
        }
      return e;
    }
  else if (lex_match (s->lexer, T_LCURLY))
    {
      struct matrix_expr *e = matrix_parse_curly_semi (s);
      if (!e || !lex_force_match (s->lexer, T_RCURLY))
        {
          matrix_expr_destroy (e);
          return NULL;
        }
      return e;
    }
  else if (lex_token (s->lexer) == T_ID)
    {
      struct matrix_expr *retval;
      if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
        return retval;

      struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
      if (!var)
        {
          lex_error (s->lexer, _("Unknown variable %s."),
                     lex_tokcstr (s->lexer));
          return NULL;
        }
      lex_get (s->lexer);

      struct matrix_expr *e = xmalloc (sizeof *e);
      *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
      return e;
    }
  else if (lex_token (s->lexer) == T_ALL)
    {
      struct matrix_expr *retval;
      if (matrix_parse_function (s, "ALL", &retval))
        return retval;
    }

  lex_error (s->lexer, NULL);
  return NULL;
}

static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);

static bool
matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
{
  if (lex_match (s->lexer, T_COLON))
    {
      *indexp = NULL;
      return true;
    }
  else
    {
      *indexp = matrix_parse_expr (s);
      return *indexp != NULL;
    }
}

static struct matrix_expr *
matrix_parse_postfix (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_primary (s);
  if (!lhs || !lex_match (s->lexer, T_LPAREN))
    return lhs;

  struct matrix_expr *i0;
  if (!matrix_parse_index_expr (s, &i0))
    {
      matrix_expr_destroy (lhs);
      return NULL;
    }
  if (lex_match (s->lexer, T_RPAREN))
    return (i0
            ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
            : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
  else if (lex_match (s->lexer, T_COMMA))
    {
      struct matrix_expr *i1;
      if (!matrix_parse_index_expr (s, &i1)
          || !lex_force_match (s->lexer, T_RPAREN))
        {
          matrix_expr_destroy (lhs);
          matrix_expr_destroy (i0);
          matrix_expr_destroy (i1);
          return NULL;
        }
      return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
              : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
              : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
              : lhs);
    }
  else
    {
      lex_error_expecting (s->lexer, "`)'", "`,'");
      return NULL;
    }
}

static struct matrix_expr *
matrix_parse_unary (struct matrix_state *s)
{
  if (lex_match (s->lexer, T_DASH))
    {
      struct matrix_expr *lhs = matrix_parse_unary (s);
      return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
    }
  else if (lex_match (s->lexer, T_PLUS))
    return matrix_parse_unary (s);
  else
    return matrix_parse_postfix (s);
}

static struct matrix_expr *
matrix_parse_seq (struct matrix_state *s)
{
  struct matrix_expr *start = matrix_parse_unary (s);
  if (!start || !lex_match (s->lexer, T_COLON))
    return start;

  struct matrix_expr *end = matrix_parse_unary (s);
  if (!end)
    {
      matrix_expr_destroy (start);
      return NULL;
    }

  if (lex_match (s->lexer, T_COLON))
    {
      struct matrix_expr *increment = matrix_parse_unary (s);
      if (!increment)
        {
          matrix_expr_destroy (start);
          matrix_expr_destroy (end);
          return NULL;
        }
      return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
    }
  else
    return matrix_expr_create_2 (MOP_SEQ, start, end);
}

static struct matrix_expr *
matrix_parse_exp (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_seq (s);
  if (!lhs)
    return NULL;

  for (;;)
    {
      enum matrix_op op;
      if (lex_match (s->lexer, T_EXP))
        op = MOP_EXP_MAT;
      else if (lex_match_phrase (s->lexer, "&**"))
        op = MOP_EXP_ELEMS;
      else
        return lhs;

      struct matrix_expr *rhs = matrix_parse_seq (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (op, lhs, rhs);
    }
}

static struct matrix_expr *
matrix_parse_mul_div (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_exp (s);
  if (!lhs)
    return NULL;

  for (;;)
    {
      enum matrix_op op;
      if (lex_match (s->lexer, T_ASTERISK))
        op = MOP_MUL_MAT;
      else if (lex_match (s->lexer, T_SLASH))
        op = MOP_DIV_ELEMS;
      else if (lex_match_phrase (s->lexer, "&*"))
        op = MOP_MUL_ELEMS;
      else if (lex_match_phrase (s->lexer, "&/"))
        op = MOP_DIV_ELEMS;
      else
        return lhs;

      struct matrix_expr *rhs = matrix_parse_exp (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (op, lhs, rhs);
    }
}

static struct matrix_expr *
matrix_parse_add_sub (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_mul_div (s);
  if (!lhs)
    return NULL;

  for (;;)
    {
      enum matrix_op op;
      if (lex_match (s->lexer, T_PLUS))
        op = MOP_ADD_ELEMS;
      else if (lex_match (s->lexer, T_DASH))
        op = MOP_SUB_ELEMS;
      else if (lex_token (s->lexer) == T_NEG_NUM)
        op = MOP_ADD_ELEMS;
      else
        return lhs;

      struct matrix_expr *rhs = matrix_parse_mul_div (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (op, lhs, rhs);
    }
}

static struct matrix_expr *
matrix_parse_relational (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_add_sub (s);
  if (!lhs)
    return NULL;

  for (;;)
    {
      enum matrix_op op;
      if (lex_match (s->lexer, T_GT))
        op = MOP_GT;
      else if (lex_match (s->lexer, T_GE))
        op = MOP_GE;
      else if (lex_match (s->lexer, T_LT))
        op = MOP_LT;
      else if (lex_match (s->lexer, T_LE))
        op = MOP_LE;
      else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
        op = MOP_EQ;
      else if (lex_match (s->lexer, T_NE))
        op = MOP_NE;
      else
        return lhs;

      struct matrix_expr *rhs = matrix_parse_add_sub (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (op, lhs, rhs);
    }
}

static struct matrix_expr *
matrix_parse_not (struct matrix_state *s)
{
  if (lex_match (s->lexer, T_NOT))
    {
      struct matrix_expr *lhs = matrix_parse_not (s);
      return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
    }
  else
    return matrix_parse_relational (s);
}

static struct matrix_expr *
matrix_parse_and (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_not (s);
  if (!lhs)
    return NULL;

  while (lex_match (s->lexer, T_AND))
    {
      struct matrix_expr *rhs = matrix_parse_not (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
    }
  return lhs;
}

static struct matrix_expr *
matrix_parse_or_xor (struct matrix_state *s)
{
  struct matrix_expr *lhs = matrix_parse_and (s);
  if (!lhs)
    return NULL;

  for (;;)
    {
      enum matrix_op op;
      if (lex_match (s->lexer, T_OR))
        op = MOP_OR;
      else if (lex_match_id (s->lexer, "XOR"))
        op = MOP_XOR;
      else
        return lhs;

      struct matrix_expr *rhs = matrix_parse_and (s);
      if (!rhs)
        {
          matrix_expr_destroy (lhs);
          return NULL;
        }
      lhs = matrix_expr_create_2 (op, lhs, rhs);
    }
}

static struct matrix_expr *
matrix_parse_expr (struct matrix_state *s)
{
  return matrix_parse_or_xor (s);
}
\f
/* Expression evaluation. */

static double
matrix_op_eval (enum matrix_op op, double a, double b)
{
  switch (op)
    {
    case MOP_ADD_ELEMS: return a + b;
    case MOP_SUB_ELEMS: return a - b;
    case MOP_MUL_ELEMS: return a * b;
    case MOP_DIV_ELEMS: return a / b;
    case MOP_EXP_ELEMS: return pow (a, b);
    case MOP_GT: return a > b;
    case MOP_GE: return a >= b;
    case MOP_LT: return a < b;
    case MOP_LE: return a <= b;
    case MOP_EQ: return a == b;
    case MOP_NE: return a != b;
    case MOP_AND: return (a > 0) && (b > 0);
    case MOP_OR: return (a > 0) || (b > 0);
    case MOP_XOR: return (a > 0) != (b > 0);

#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
      MATRIX_FUNCTIONS
#undef F
    case MOP_NEGATE:
    case MOP_SEQ:
    case MOP_SEQ_BY:
    case MOP_MUL_MAT:
    case MOP_EXP_MAT:
    case MOP_NOT:
    case MOP_PASTE_HORZ:
    case MOP_PASTE_VERT:
    case MOP_EMPTY:
    case MOP_VEC_INDEX:
    case MOP_VEC_ALL:
    case MOP_MAT_INDEX:
    case MOP_ROW_INDEX:
    case MOP_COL_INDEX:
    case MOP_NUMBER:
    case MOP_VARIABLE:
    case MOP_EOF:
      NOT_REACHED ();
    }
  NOT_REACHED ();
}

static const char *
matrix_op_name (enum matrix_op op)
{
  switch (op)
    {
    case MOP_ADD_ELEMS: return "+";
    case MOP_SUB_ELEMS: return "-";
    case MOP_MUL_ELEMS: return "&*";
    case MOP_DIV_ELEMS: return "&/";
    case MOP_EXP_ELEMS: return "&**";
    case MOP_GT: return ">";
    case MOP_GE: return ">=";
    case MOP_LT: return "<";
    case MOP_LE: return "<=";
    case MOP_EQ: return "=";
    case MOP_NE: return "<>";
    case MOP_AND: return "AND";
    case MOP_OR: return "OR";
    case MOP_XOR: return "XOR";

#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
      MATRIX_FUNCTIONS
#undef F
    case MOP_NEGATE:
    case MOP_SEQ:
    case MOP_SEQ_BY:
    case MOP_MUL_MAT:
    case MOP_EXP_MAT:
    case MOP_NOT:
    case MOP_PASTE_HORZ:
    case MOP_PASTE_VERT:
    case MOP_EMPTY:
    case MOP_VEC_INDEX:
    case MOP_VEC_ALL:
    case MOP_MAT_INDEX:
    case MOP_ROW_INDEX:
    case MOP_COL_INDEX:
    case MOP_NUMBER:
    case MOP_VARIABLE:
    case MOP_EOF:
      NOT_REACHED ();
    }
  NOT_REACHED ();
}

static bool
is_scalar (const gsl_matrix *m)
{
  return m->size1 == 1 && m->size2 == 1;
}

static double
to_scalar (const gsl_matrix *m)
{
  assert (is_scalar (m));
  return gsl_matrix_get (m, 0, 0);
}

static gsl_matrix *
matrix_expr_evaluate_elementwise (enum matrix_op op,
                                  gsl_matrix *a, gsl_matrix *b)
{
  if (is_scalar (b))
    {
      double be = to_scalar (b);
      for (size_t r = 0; r < a->size1; r++)
        for (size_t c = 0; c < a->size2; c++)
          {
            double *ae = gsl_matrix_ptr (a, r, c);
            *ae = matrix_op_eval (op, *ae, be);
          }
      return a;
    }
  else if (is_scalar (a))
    {
      double ae = to_scalar (a);
      for (size_t r = 0; r < b->size1; r++)
        for (size_t c = 0; c < b->size2; c++)
          {
            double *be = gsl_matrix_ptr (b, r, c);
            *be = matrix_op_eval (op, ae, *be);
          }
      return b;
    }
  else if (a->size1 == b->size1 && a->size2 == b->size2)
    {
      for (size_t r = 0; r < a->size1; r++)
        for (size_t c = 0; c < a->size2; c++)
          {
            double *ae = gsl_matrix_ptr (a, r, c);
            double be = gsl_matrix_get (b, r, c);
            *ae = matrix_op_eval (op, *ae, be);
          }
      return a;
    }
  else
    {
      msg (SE, _("Operands to %s must have the same dimensions or one "
                 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
           matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
      return NULL;
    }
}

static gsl_matrix *
matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
{
  if (is_scalar (a) || is_scalar (b))
    return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);

  if (a->size2 != b->size1)
    {
      msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
                 "not conformable for multiplication."),
           a->size1, a->size2, b->size1, b->size2);
      return NULL;
    }

  gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
  return c;
}

static void
swap_matrix (gsl_matrix **a, gsl_matrix **b)
{
  gsl_matrix *tmp = *a;
  *a = *b;
  *b = tmp;
}

static void
mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
            gsl_matrix **tmp)
{
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
  swap_matrix (z, tmp);
}

static void
square_matrix (gsl_matrix **x, gsl_matrix **tmp)
{
  mul_matrix (x, *x, *x, tmp);
}

static gsl_matrix *
matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
{
  gsl_matrix *x = x_;
  if (x->size1 != x->size2)
    {
      msg (SE, _("Matrix exponentation with ** requires a square matrix on "
                 "the left-hand size, not one with dimensions %zu×%zu."),
           x->size1, x->size2);
      return NULL;
    }
  if (!is_scalar (b))
    {
      msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
                 "right-hand side, not a matrix with dimensions %zu×%zu."),
           b->size1, b->size2);
      return NULL;
    }
  double bf = to_scalar (b);
  if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
    {
      msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
                 "or outside the valid range."), bf);
      return NULL;
    }
  long int bl = bf;

  gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
  gsl_matrix *y = y_;
  gsl_matrix_set_identity (y);
  if (bl == 0)
    return y;

  gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
  gsl_matrix *t = t_;
  for (unsigned long int n = labs (bl); n > 1; n /= 2)
    if (n & 1)
      {
        mul_matrix (&y, x, y, &t);
        square_matrix (&x, &t);
      }
    else
      square_matrix (&x, &t);

  mul_matrix (&y, x, y, &t);
  if (bf < 0)
    invert_matrix (y);

  /* Garbage collection.

     There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
     a permutation of them.  We are returning one of them; that one must not be
     destroyed.  We must not destroy 'x_' because the caller owns it. */
  if (y != y_)
    gsl_matrix_free (y_);
  if (y != t_)
    gsl_matrix_free (t_);

  return y;
}

static gsl_matrix *
matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
                          gsl_matrix *by_)
{
  if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
    {
      msg (SE, _("All operands of : operator must be scalars."));
      return NULL;
    }

  long int start = to_scalar (start_);
  long int end = to_scalar (end_);
  long int by = by_ ? to_scalar (by_) : 1;

  if (!by)
    {
      msg (SE, _("The increment operand to : must be nonzero."));
      return NULL;
    }

  long int n = (end >= start && by > 0 ? (end - start + by) / by
                : end <= start && by < 0 ? (start - end - by) / -by
                : 0);
  gsl_matrix *m = gsl_matrix_alloc (1, n);
  for (long int i = 0; i < n; i++)
    gsl_matrix_set (m, 0, i, start + i * by);
  return m;
}

static gsl_matrix *
matrix_expr_evaluate_not (gsl_matrix *a)
{
  for (size_t r = 0; r < a->size1; r++)
    for (size_t c = 0; c < a->size2; c++)
      {
        double *ae = gsl_matrix_ptr (a, r, c);
        *ae = !(*ae > 0);
      }
  return a;
}

static gsl_matrix *
matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
{
  if (a->size1 != b->size1)
    {
      if (!a->size1 || !a->size2)
        return b;
      else if (!b->size1 || !b->size2)
        return a;

      msg (SE, _("All columns in a matrix must have the same number of rows, "
                 "but this tries to paste matrices with %zu and %zu rows."),
           a->size1, b->size1);
      return NULL;
    }

  gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
  for (size_t y = 0; y < a->size1; y++)
    {
      for (size_t x = 0; x < a->size2; x++)
        gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
      for (size_t x = 0; x < b->size2; x++)
        gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
    }
  return c;
}

static gsl_matrix *
matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
{
  if (a->size2 != b->size2)
    {
      if (!a->size1 || !a->size2)
        return b;
      else if (!b->size1 || !b->size2)
        return a;

      msg (SE, _("All rows in a matrix must have the same number of columns, "
                 "but this tries to stack matrices with %zu and %zu columns."),
           a->size2, b->size2);
      return NULL;
    }

  gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
  for (size_t x = 0; x < a->size2; x++)
    {
      for (size_t y = 0; y < a->size1; y++)
        gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
      for (size_t y = 0; y < b->size1; y++)
        gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
    }
  return c;
}

static gsl_vector *
matrix_to_vector (gsl_matrix *m)
{
  assert (m->owner);
  gsl_vector v = to_vector (m);
  assert (v.block == m->block || !v.block);
  assert (!v.owner);
  v.owner = 1;
  m->owner = 0;
  return xmemdup (&v, sizeof v);
}

enum index_type {
  IV_ROW,
  IV_COLUMN,
  IV_VECTOR
};

struct index_vector
  {
    size_t *indexes;
    size_t n;
  };
#define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }

static bool
matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
                               enum index_type index_type, size_t other_size,
                               struct index_vector *iv)
{
  if (m)
    {
      if (!is_vector (m))
        {
          switch (index_type)
            {
            case IV_VECTOR:
              msg (SE, _("Vector index must be scalar or vector, not a "
                         "%zu×%zu matrix."),
                   m->size1, m->size2);
              break;

            case IV_ROW:
              msg (SE, _("Matrix row index must be scalar or vector, not a "
                         "%zu×%zu matrix."),
                   m->size1, m->size2);
              break;

            case IV_COLUMN:
              msg (SE, _("Matrix column index must be scalar or vector, not a "
                         "%zu×%zu matrix."),
                   m->size1, m->size2);
              break;
            }
          return false;
        }

      gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
      *iv = (struct index_vector) {
        .indexes = xnmalloc (v.size, sizeof *iv->indexes),
        .n = v.size,
      };
      for (size_t i = 0; i < v.size; i++)
        {
          double index = gsl_vector_get (&v, i);
          if (index < 1 || index >= size + 1)
            {
              switch (index_type)
                {
                case IV_VECTOR:
                  msg (SE, _("Index %g is out of range for vector "
                             "with %zu elements."), index, size);
                  break;

                case IV_ROW:
                  msg (SE, _("%g is not a valid row index for "
                             "a %zu×%zu matrix."),
                       index, size, other_size);
                  break;

                case IV_COLUMN:
                  msg (SE, _("%g is not a valid column index for "
                             "a %zu×%zu matrix."),
                       index, other_size, size);
                  break;
                }

              free (iv->indexes);
              return false;
            }
          iv->indexes[i] = index - 1;
        }
      return true;
    }
  else
    {
      *iv = (struct index_vector) {
        .indexes = xnmalloc (size, sizeof *iv->indexes),
        .n = size,
      };
      for (size_t i = 0; i < size; i++)
        iv->indexes[i] = i;
      return true;
    }
}

static gsl_matrix *
matrix_expr_evaluate_vec_all (gsl_matrix *sm)
{
  if (!is_vector (sm))
    {
      msg (SE, _("Vector index operator may not be applied to "
                 "a %zu×%zu matrix."),
           sm->size1, sm->size2);
      return NULL;
    }

  return sm;
}

static gsl_matrix *
matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
{
  if (!matrix_expr_evaluate_vec_all (sm))
    return NULL;

  gsl_vector sv = to_vector (sm);
  struct index_vector iv;
  if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
    return NULL;

  gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
                                     sm->size1 == 1 ? iv.n : 1);
  gsl_vector dv = to_vector (dm);
  for (size_t dx = 0; dx < iv.n; dx++)
    {
      size_t sx = iv.indexes[dx];
      gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
    }
  free (iv.indexes);

  return dm;
}

static gsl_matrix *
matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
                                gsl_matrix *im1)
{
  struct index_vector iv0;
  if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
    return NULL;

  struct index_vector iv1;
  if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
                                      &iv1))
    {
      free (iv0.indexes);
      return NULL;
    }

  gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
  for (size_t dy = 0; dy < iv0.n; dy++)
    {
      size_t sy = iv0.indexes[dy];

      for (size_t dx = 0; dx < iv1.n; dx++)
        {
          size_t sx = iv1.indexes[dx];
          gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
        }
    }
  free (iv0.indexes);
  free (iv1.indexes);
  return dm;
}

#define F(NAME, PROTOTYPE)                              \
  static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
    const char *name, gsl_matrix *[], size_t,           \
    matrix_proto_##PROTOTYPE *);
MATRIX_FUNCTIONS
#undef F

static bool
check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
{
  if (!is_scalar (subs[index]))
    {
      msg (SE, _("Function %s argument %zu must be a scalar, "
                 "not a %zu×%zu matrix."),
           name, index + 1, subs[index]->size1, subs[index]->size2);
      return false;
    }
  return true;
}

static bool
check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
{
  if (!is_vector (subs[index]))
    {
      msg (SE, _("Function %s argument %zu must be a vector, "
                 "not a %zu×%zu matrix."),
           name, index + 1, subs[index]->size1, subs[index]->size2);
      return false;
    }
  return true;
}

static bool
to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
{
  for (size_t i = 0; i < n_subs; i++)
    {
      if (!check_scalar_arg (name, subs, i))
        return false;
      d[i] = to_scalar (subs[i]);
    }
  return true;
}

static gsl_matrix *
matrix_expr_evaluate_m_d (const char *name,
                          gsl_matrix *subs[], size_t n_subs,
                          matrix_proto_m_d *f)
{
  assert (n_subs == 1);

  double d;
  return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
}

static gsl_matrix *
matrix_expr_evaluate_m_dd (const char *name,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_dd *f)
{
  assert (n_subs == 2);

  double d[2];
  return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
}

static gsl_matrix *
matrix_expr_evaluate_m_ddd (const char *name,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_ddd *f)
{
  assert (n_subs == 3);

  double d[3];
  return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
}

static gsl_matrix *
matrix_expr_evaluate_m_m (const char *name UNUSED,
                          gsl_matrix *subs[], size_t n_subs,
                          matrix_proto_m_m *f)
{
  assert (n_subs == 1);
  return f (subs[0]);
}

static gsl_matrix *
matrix_expr_evaluate_m_md (const char *name UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_md *f)
{
  assert (n_subs == 2);
  if (!check_scalar_arg (name, subs, 1))
    return NULL;
  return f (subs[0], to_scalar (subs[1]));
}

static gsl_matrix *
matrix_expr_evaluate_m_mdd (const char *name UNUSED,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_mdd *f)
{
  assert (n_subs == 3);
  if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
    return NULL;
  return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
}

static gsl_matrix *
matrix_expr_evaluate_m_mm (const char *name UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_mm *f)
{
  assert (n_subs == 2);
  return f (subs[0], subs[1]);
}

static gsl_matrix *
matrix_expr_evaluate_m_v (const char *name,
                          gsl_matrix *subs[], size_t n_subs,
                          matrix_proto_m_v *f)
{
  assert (n_subs == 1);
  if (!check_vector_arg (name, subs, 0))
    return NULL;
  gsl_vector v = to_vector (subs[0]);
  return f (&v);
}

static gsl_matrix *
matrix_expr_evaluate_d_m (const char *name UNUSED,
                          gsl_matrix *subs[], size_t n_subs,
                          matrix_proto_d_m *f)
{
  assert (n_subs == 1);

  gsl_matrix *m = gsl_matrix_alloc (1, 1);
  gsl_matrix_set (m, 0, 0, f (subs[0]));
  return m;
}

static gsl_matrix *
matrix_expr_evaluate_m_any (const char *name UNUSED,
                          gsl_matrix *subs[], size_t n_subs,
                          matrix_proto_m_any *f)
{
  return f (subs, n_subs);
}

static gsl_matrix *
matrix_expr_evaluate_IDENT (const char *name,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_IDENT *f)
{
  assert (n_subs <= 2);

  double d[2];
  if (!to_scalar_args (name, subs, n_subs, d))
    return NULL;

  return f (d[0], d[n_subs - 1]);
}

static gsl_matrix *
matrix_expr_evaluate (const struct matrix_expr *e)
{
  if (e->op == MOP_NUMBER)
    {
      gsl_matrix *m = gsl_matrix_alloc (1, 1);
      gsl_matrix_set (m, 0, 0, e->number);
      return m;
    }
  else if (e->op == MOP_VARIABLE)
    {
      const gsl_matrix *src = e->variable->value;
      if (!src)
        {
          msg (SE, _("Uninitialized variable %s used in expression."),
               e->variable->name);
          return NULL;
        }

      gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
      gsl_matrix_memcpy (dst, src);
      return dst;
    }
  else if (e->op == MOP_EOF)
    {
      struct dfm_reader *reader = read_file_open (e->eof);
      gsl_matrix *m = gsl_matrix_alloc (1, 1);
      gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
      return m;
    }

  enum { N_LOCAL = 3 };
  gsl_matrix *local_subs[N_LOCAL];
  gsl_matrix **subs = (e->n_subs < N_LOCAL
                       ? local_subs
                       : xmalloc (e->n_subs * sizeof *subs));

  for (size_t i = 0; i < e->n_subs; i++)
    {
      subs[i] = matrix_expr_evaluate (e->subs[i]);
      if (!subs[i])
        {
          for (size_t j = 0; j < i; j++)
            gsl_matrix_free (subs[j]);
          if (subs != local_subs)
            free (subs);
          return NULL;
        }
    }

  gsl_matrix *result = NULL;
  switch (e->op)
    {
#define F(NAME, PROTOTYPE)                                              \
      case MOP_F_##NAME:                                                \
        result = matrix_expr_evaluate_##PROTOTYPE (#NAME,               \
                                                   subs, e->n_subs,     \
                                                   matrix_eval_##NAME); \
      break;
      MATRIX_FUNCTIONS
#undef F

    case MOP_NEGATE:
      gsl_matrix_scale (subs[0], -1.0);
      result = subs[0];
      break;

    case MOP_ADD_ELEMS:
    case MOP_SUB_ELEMS:
    case MOP_MUL_ELEMS:
    case MOP_DIV_ELEMS:
    case MOP_EXP_ELEMS:
    case MOP_GT:
    case MOP_GE:
    case MOP_LT:
    case MOP_LE:
    case MOP_EQ:
    case MOP_NE:
    case MOP_AND:
    case MOP_OR:
    case MOP_XOR:
      result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
      break;

    case MOP_NOT:
      result = matrix_expr_evaluate_not (subs[0]);
      break;

    case MOP_SEQ:
      result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
      break;

    case MOP_SEQ_BY:
      result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
      break;

    case MOP_MUL_MAT:
      result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
      break;

    case MOP_EXP_MAT:
      result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
      break;

    case MOP_PASTE_HORZ:
      result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
      break;

    case MOP_PASTE_VERT:
      result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
      break;

    case MOP_EMPTY:
      result = gsl_matrix_alloc (0, 0);
      break;

    case MOP_VEC_INDEX:
      result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
      break;

    case MOP_VEC_ALL:
      result = matrix_expr_evaluate_vec_all (subs[0]);
      break;

    case MOP_MAT_INDEX:
      result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
      break;

    case MOP_ROW_INDEX:
      result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
      break;

    case MOP_COL_INDEX:
      result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
      break;

    case MOP_NUMBER:
    case MOP_VARIABLE:
    case MOP_EOF:
      NOT_REACHED ();
    }

  for (size_t i = 0; i < e->n_subs; i++)
    if (subs[i] != result)
      gsl_matrix_free (subs[i]);
  if (subs != local_subs)
    free (subs);
  return result;
}

static bool
matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
                             double *d)
{
  gsl_matrix *m = matrix_expr_evaluate (e);
  if (!m)
    return false;

  if (!is_scalar (m))
    {
      msg (SE, _("Expression for %s must evaluate to scalar, "
                 "not a %zu×%zu matrix."),
           context, m->size1, m->size2);
      gsl_matrix_free (m);
      return false;
    }

  *d = to_scalar (m);
  gsl_matrix_free (m);
  return true;
}

static bool
matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
                              long int *integer)
{
  double d;
  if (!matrix_expr_evaluate_scalar (e, context, &d))
    return false;

  d = trunc (d);
  if (d < LONG_MIN || d > LONG_MAX)
    {
      msg (SE, _("Expression for %s is outside the integer range."), context);
      return false;
    }
  *integer = d;
  return true;
}
\f
struct matrix_lvalue
  {
    struct matrix_var *var;
    struct matrix_expr *indexes[2];
    size_t n_indexes;
  };

static void
matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
{
  if (lvalue)
    {
      for (size_t i = 0; i < lvalue->n_indexes; i++)
        matrix_expr_destroy (lvalue->indexes[i]);
      free (lvalue);
    }
}

static struct matrix_lvalue *
matrix_lvalue_parse (struct matrix_state *s)
{
  struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
  if (!lex_force_id (s->lexer))
    goto error;
  lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
  if (lex_next_token (s->lexer, 1) == T_LPAREN)
    {
      if (!lvalue->var)
        {
          msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
          free (lvalue);
          return NULL;
        }

      lex_get_n (s->lexer, 2);

      if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
        goto error;
      lvalue->n_indexes++;

      if (lex_match (s->lexer, T_COMMA))
        {
          if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
            goto error;
          lvalue->n_indexes++;
        }
      if (!lex_force_match (s->lexer, T_RPAREN))
        goto error;
    }
  else
    {
      if (!lvalue->var)
        lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
      lex_get (s->lexer);
    }
  return lvalue;

error:
  matrix_lvalue_destroy (lvalue);
  return NULL;
}

static bool
matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
                               enum index_type index_type, size_t other_size,
                               struct index_vector *iv)
{
  gsl_matrix *m;
  if (e)
    {
      m = matrix_expr_evaluate (e);
      if (!m)
        return false;
    }
  else
    m = NULL;

  bool ok = matrix_normalize_index_vector (m, size, index_type,
                                           other_size, iv);
  gsl_matrix_free (m);
  return ok;
}

static bool
matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
                             struct index_vector *iv, gsl_matrix *sm)
{
  gsl_vector dv = to_vector (lvalue->var->value);

  /* Convert source matrix 'sm' to source vector 'sv'. */
  if (!is_vector (sm))
    {
      msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
           sm->size1, sm->size2);
      return false;
    }
  gsl_vector sv = to_vector (sm);
  if (iv->n != sv.size)
    {
      msg (SE, _("Can't assign %zu-element vector "
                 "to %zu-element subvector."), sv.size, iv->n);
      return false;
    }

  /* Assign elements. */
  for (size_t x = 0; x < iv->n; x++)
    gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
  return true;
}

static bool
matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
                             struct index_vector *iv0,
                             struct index_vector *iv1, gsl_matrix *sm)
{
  gsl_matrix *dm = lvalue->var->value;

  /* Convert source matrix 'sm' to source vector 'sv'. */
  if (iv0->n != sm->size1)
    {
      msg (SE, _("Row index vector for assignment to %s has %zu elements "
                 "but source matrix has %zu rows."),
           lvalue->var->name, iv0->n, sm->size1);
      return false;
    }
  if (iv1->n != sm->size2)
    {
      msg (SE, _("Column index vector for assignment to %s has %zu elements "
                 "but source matrix has %zu columns."),
           lvalue->var->name, iv1->n, sm->size2);
      return false;
    }

  /* Assign elements. */
  for (size_t y = 0; y < iv0->n; y++)
    {
      size_t f0 = iv0->indexes[y];
      for (size_t x = 0; x < iv1->n; x++)
        {
          size_t f1 = iv1->indexes[x];
          gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
        }
    }
  return true;
}

/* Takes ownership of M. */
static bool
matrix_lvalue_assign (struct matrix_lvalue *lvalue,
                      struct index_vector *iv0, struct index_vector *iv1,
                      gsl_matrix *sm)
{
  if (!lvalue->n_indexes)
    {
      gsl_matrix_free (lvalue->var->value);
      lvalue->var->value = sm;
      return true;
    }
  else
    {
      bool ok = (lvalue->n_indexes == 1
                 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
                 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
      free (iv0->indexes);
      free (iv1->indexes);
      gsl_matrix_free (sm);
      return ok;
    }
}

static bool
matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
                        struct index_vector *iv0,
                        struct index_vector *iv1)
{
  *iv0 = INDEX_VECTOR_INIT;
  *iv1 = INDEX_VECTOR_INIT;
  if (!lvalue->n_indexes)
    return true;

  /* Validate destination matrix exists and has the right shape. */
  gsl_matrix *dm = lvalue->var->value;
  if (!dm)
    {
      msg (SE, _("Undefined variable %s."), lvalue->var->name);
      return false;
    }
  else if (dm->size1 == 0 || dm->size2 == 0)
    {
      msg (SE, _("Cannot index %zu×%zu matrix."),
           dm->size1, dm->size2);
      return false;
    }
  else if (lvalue->n_indexes == 1)
    {
      if (!is_vector (dm))
        {
          msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
               dm->size1, dm->size2, lvalue->var->name);
          return false;
        }
      return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
                                            MAX (dm->size1, dm->size2),
                                            IV_VECTOR, 0, iv0);
    }
  else
    {
      assert (lvalue->n_indexes == 2);
      if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
                                          IV_ROW, dm->size2, iv0))
        return false;

      if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
                                          IV_COLUMN, dm->size1, iv1))
        {
          free (iv0->indexes);
          return false;
        }
      return true;
    }
}

/* Takes ownership of M. */
static bool
matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
{
  struct index_vector iv0, iv1;
  if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
    {
      gsl_matrix_free (sm);
      return false;
    }

  return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
}

\f
/* Matrix command. */

struct matrix_cmds
  {
    struct matrix_cmd **commands;
    size_t n;
  };

static void matrix_cmds_uninit (struct matrix_cmds *);

struct matrix_cmd
  {
    enum matrix_cmd_type
      {
        MCMD_COMPUTE,
        MCMD_PRINT,
        MCMD_DO_IF,
        MCMD_LOOP,
        MCMD_BREAK,
        MCMD_DISPLAY,
        MCMD_RELEASE,
        MCMD_SAVE,
        MCMD_READ,
        MCMD_WRITE,
        MCMD_GET,
        MCMD_MSAVE,
        MCMD_MGET,
        MCMD_EIGEN,
        MCMD_SETDIAG,
        MCMD_SVD,
      }
    type;

    union
      {
        struct compute_command
          {
            struct matrix_lvalue *lvalue;
            struct matrix_expr *rvalue;
          }
        compute;

        struct print_command
          {
            struct matrix_expr *expression;
            bool use_default_format;
            struct fmt_spec format;
            char *title;
            int space;          /* -1 means NEWPAGE. */

            struct print_labels
              {
                struct string_array literals; /* CLABELS/RLABELS. */
                struct matrix_expr *expr;     /* CNAMES/RNAMES. */
              }
            *rlabels, *clabels;
          }
        print;

        struct do_if_command
          {
            struct do_if_clause
              {
                struct matrix_expr *condition;
                struct matrix_cmds commands;
              }
            *clauses;
            size_t n_clauses;
          }
        do_if;

        struct loop_command
          {
            /* Index. */
            struct matrix_var *var;
            struct matrix_expr *start;
            struct matrix_expr *end;
            struct matrix_expr *increment;

            /* Loop conditions. */
            struct matrix_expr *top_condition;
            struct matrix_expr *bottom_condition;

            /* Commands. */
            struct matrix_cmds commands;
          }
        loop;

        struct display_command
          {
            struct matrix_state *state;
            bool dictionary;
            bool status;
          }
        display;

        struct release_command
          {
            struct matrix_var **vars;
            size_t n_vars;
          }
        release;

        struct save_command
          {
            struct matrix_expr *expression;
            struct file_handle *outfile;
            struct string_array *variables;
            struct matrix_expr *names;
            struct stringi_set strings;
          }
        save;

        struct read_command
          {
            struct read_file *rf;
            struct matrix_lvalue *dst;
            struct matrix_expr *size;
            int c1, c2;
            enum fmt_type format;
            int w;
            int d;
            bool symmetric;
            bool reread;
          }
        read;

        struct write_command
          {
            struct write_file *wf;
            struct matrix_expr *expression;
            int c1, c2;
            enum fmt_type format;
            int w;
            int d;
            bool triangular;
            bool hold;          /* XXX */
          }
        write;

        struct get_command
          {
            struct matrix_lvalue *dst;
            struct file_handle *file;
            char *encoding;
            struct string_array variables;
            struct matrix_var *names;

            /* Treatment of missing values. */
            struct
              {
                enum
                  {
                    MGET_ERROR,  /* Flag error on command. */
                    MGET_ACCEPT, /* Accept without change, user-missing only. */
                    MGET_OMIT,   /* Drop this case. */
                    MGET_RECODE  /* Recode to 'substitute'. */
                  }
                treatment;
                double substitute; /* MGET_RECODE only. */
              }
            user, system;
          }
        get;

        struct msave_command
          {
            struct msave_common *common;
            char *varname_;
            struct matrix_expr *expr;
            const char *rowtype;
            struct matrix_expr *factors;
            struct matrix_expr *splits;
          }
         msave;

        struct mget_command
          {
            struct matrix_state *state;
            struct file_handle *file;
            struct stringi_set rowtypes;
          }
        mget;

        struct eigen_command
          {
            struct matrix_expr *expr;
            struct matrix_var *evec;
            struct matrix_var *eval;
          }
        eigen;

        struct setdiag_command
          {
            struct matrix_var *dst;
            struct matrix_expr *expr;
          }
        setdiag;

        struct svd_command
          {
            struct matrix_expr *expr;
            struct matrix_var *u;
            struct matrix_var *s;
            struct matrix_var *v;
          }
        svd;
      };
  };

static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
static bool matrix_cmd_execute (struct matrix_cmd *);
static void matrix_cmd_destroy (struct matrix_cmd *);

\f
static struct matrix_cmd *
matrix_parse_compute (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_COMPUTE,
    .compute = { .lvalue = NULL },
  };

  struct compute_command *compute = &cmd->compute;
  compute->lvalue = matrix_lvalue_parse (s);
  if (!compute->lvalue)
    goto error;

  if (!lex_force_match (s->lexer, T_EQUALS))
    goto error;

  compute->rvalue = matrix_parse_expr (s);
  if (!compute->rvalue)
    goto error;

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_compute (struct compute_command *compute)
{
  gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
  if (!value)
    return;

  matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
}
\f
static void
print_labels_destroy (struct print_labels *labels)
{
  if (labels)
    {
      string_array_destroy (&labels->literals);
      matrix_expr_destroy (labels->expr);
      free (labels);
    }
}

static void
parse_literal_print_labels (struct matrix_state *s,
                            struct print_labels **labelsp)
{
  lex_match (s->lexer, T_EQUALS);
  print_labels_destroy (*labelsp);
  *labelsp = xzalloc (sizeof **labelsp);
  while (lex_token (s->lexer) != T_SLASH
         && lex_token (s->lexer) != T_ENDCMD
         && lex_token (s->lexer) != T_STOP)
    {
      struct string label = DS_EMPTY_INITIALIZER;
      while (lex_token (s->lexer) != T_COMMA
             && lex_token (s->lexer) != T_SLASH
             && lex_token (s->lexer) != T_ENDCMD
             && lex_token (s->lexer) != T_STOP)
        {
          if (!ds_is_empty (&label))
            ds_put_byte (&label, ' ');

          if (lex_token (s->lexer) == T_STRING)
            ds_put_cstr (&label, lex_tokcstr (s->lexer));
          else
            {
              char *rep = lex_next_representation (s->lexer, 0, 0);
              ds_put_cstr (&label, rep);
              free (rep);
            }
          lex_get (s->lexer);
        }
      string_array_append_nocopy (&(*labelsp)->literals,
                                  ds_steal_cstr (&label));

      if (!lex_match (s->lexer, T_COMMA))
        break;
    }
}

static bool
parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
{
  lex_match (s->lexer, T_EQUALS);
  struct matrix_expr *e = matrix_parse_exp (s);
  if (!e)
    return false;

  print_labels_destroy (*labelsp);
  *labelsp = xzalloc (sizeof **labelsp);
  (*labelsp)->expr = e;
  return true;
}

static struct matrix_cmd *
matrix_parse_print (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_PRINT,
    .print = {
      .use_default_format = true,
    }
  };

  if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
    {
      size_t depth = 0;
      for (size_t i = 0; ; i++)
        {
          enum token_type t = lex_next_token (s->lexer, i);
          if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
            depth++;
          else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
            depth--;
          else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
            {
              if (i > 0)
                cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
              break;
            }
        }

      cmd->print.expression = matrix_parse_exp (s);
      if (!cmd->print.expression)
        goto error;
    }

  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "FORMAT"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (!parse_format_specifier (s->lexer, &cmd->print.format))
            goto error;
          cmd->print.use_default_format = false;
        }
      else if (lex_match_id (s->lexer, "TITLE"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (!lex_force_string (s->lexer))
            goto error;
          free (cmd->print.title);
          cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
          lex_get (s->lexer);
        }
      else if (lex_match_id (s->lexer, "SPACE"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (lex_match_id (s->lexer, "NEWPAGE"))
            cmd->print.space = -1;
          else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
            {
              cmd->print.space = lex_integer (s->lexer);
              lex_get (s->lexer);
            }
          else
            goto error;
        }
      else if (lex_match_id (s->lexer, "RLABELS"))
        parse_literal_print_labels (s, &cmd->print.rlabels);
      else if (lex_match_id (s->lexer, "CLABELS"))
        parse_literal_print_labels (s, &cmd->print.clabels);
      else if (lex_match_id (s->lexer, "RNAMES"))
        {
          if (!parse_expr_print_labels (s, &cmd->print.rlabels))
            goto error;
        }
      else if (lex_match_id (s->lexer, "CNAMES"))
        {
          if (!parse_expr_print_labels (s, &cmd->print.clabels))
            goto error;
        }
      else
        {
          lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
                               "RLABELS", "CLABELS", "RNAMES", "CNAMES");
          goto error;
        }

    }
  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static bool
matrix_is_integer (const gsl_matrix *m)
{
  for (size_t y = 0; y < m->size1; y++)
    for (size_t x = 0; x < m->size2; x++)
      {
        double d = gsl_matrix_get (m, y, x);
        if (d != floor (d))
          return false;
      }
  return true;
}

static double
matrix_max_magnitude (const gsl_matrix *m)
{
  double max = 0.0;
  for (size_t y = 0; y < m->size1; y++)
    for (size_t x = 0; x < m->size2; x++)
      {
        double d = fabs (gsl_matrix_get (m, y, x));
        if (d > max)
          max = d;
      }
  return max;
}

static bool
format_fits (struct fmt_spec format, double x)
{
  char *s = data_out (&(union value) { .f = x }, NULL,
                      &format, settings_get_fmt_settings ());
  bool fits = *s != '*' && !strchr (s, 'E');
  free (s);
  return fits;
}

static struct fmt_spec
default_format (const gsl_matrix *m, int *log_scale)
{
  *log_scale = 0;

  double max = matrix_max_magnitude (m);

  if (matrix_is_integer (m))
    for (int w = 1; w <= 12; w++)
      {
        struct fmt_spec format = { .type = FMT_F, .w = w };
        if (format_fits (format, -max))
          return format;
      }

  if (max >= 1e9 || max <= 1e-4)
    {
      char s[64];
      snprintf (s, sizeof s, "%.1e", max);

      const char *e = strchr (s, 'e');
      if (e)
        *log_scale = atoi (e + 1);
    }

  return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
}

static char *
trimmed_string (double d)
{
  struct substring s = ss_buffer ((char *) &d, sizeof d);
  ss_rtrim (&s, ss_cstr (" "));
  return ss_xstrdup (s);
}

static struct string_array *
print_labels_get (const struct print_labels *labels, size_t n,
                  const char *prefix, bool truncate)
{
  if (!labels)
    return NULL;

  struct string_array *out = xzalloc (sizeof *out);
  if (labels->literals.n)
    string_array_clone (out, &labels->literals);
  else if (labels->expr)
    {
      gsl_matrix *m = matrix_expr_evaluate (labels->expr);
      if (m && is_vector (m))
        {
          gsl_vector v = to_vector (m);
          for (size_t i = 0; i < v.size; i++)
            string_array_append_nocopy (out, trimmed_string (
                                          gsl_vector_get (&v, i)));
        }
      gsl_matrix_free (m);
    }

  while (out->n < n)
    {
      if (labels->expr)
        string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
                                                    out->n + 1));
      else
        string_array_append (out, "");
    }
  while (out->n > n)
    string_array_delete (out, out->n - 1);

  if (truncate)
    for (size_t i = 0; i < out->n; i++)
      {
        char *s = out->strings[i];
        s[strnlen (s, 8)] = '\0';
      }

  return out;
}

static void
matrix_cmd_print_space (int space)
{
  if (space < 0)
    output_item_submit (page_break_item_create ());
  for (int i = 0; i < space; i++)
    output_log ("%s", "");
}

static void
matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
                       struct fmt_spec format, int log_scale)
{
  matrix_cmd_print_space (print->space);
  if (print->title)
    output_log ("%s", print->title);
  if (log_scale != 0)
    output_log ("  10 ** %d   X", log_scale);

  struct string_array *clabels = print_labels_get (print->clabels,
                                                   m->size2, "col", true);
  if (clabels && format.w < 8)
    format.w = 8;
  struct string_array *rlabels = print_labels_get (print->rlabels,
                                                   m->size1, "row", true);

  if (clabels)
    {
      struct string line = DS_EMPTY_INITIALIZER;
      if (rlabels)
        ds_put_byte_multiple (&line, ' ', 8);
      for (size_t x = 0; x < m->size2; x++)
        ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
      output_log_nocopy (ds_steal_cstr (&line));
    }

  double scale = pow (10.0, log_scale);
  bool numeric = fmt_is_numeric (format.type);
  for (size_t y = 0; y < m->size1; y++)
    {
      struct string line = DS_EMPTY_INITIALIZER;
      if (rlabels)
        ds_put_format (&line, "%-8s", rlabels->strings[y]);

      for (size_t x = 0; x < m->size2; x++)
        {
          double f = gsl_matrix_get (m, y, x);
          char *s = (numeric
                     ? data_out (&(union value) { .f = f / scale}, NULL,
                                 &format, settings_get_fmt_settings ())
                     : trimmed_string (f));
          ds_put_format (&line, " %s", s);
          free (s);
        }
      output_log_nocopy (ds_steal_cstr (&line));
    }

  string_array_destroy (rlabels);
  free (rlabels);
  string_array_destroy (clabels);
  free (clabels);
}

static void
create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
                        const struct print_labels *print_labels, size_t n,
                        const char *name, const char *prefix)
{
  struct string_array *labels = print_labels_get (print_labels, n, prefix,
                                                  false);
  struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
  for (size_t i = 0; i < n; i++)
    pivot_category_create_leaf (
      d->root, (labels
                ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
                : pivot_value_new_integer (i + 1)));
  if (!labels)
    d->hide_all_labels = true;
  string_array_destroy (labels);
  free (labels);
}

static void
matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
                         struct fmt_spec format, int log_scale)
{
  struct pivot_table *table = pivot_table_create__ (
    pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
    "Matrix Print");

  create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
                          N_("Rows"), "row");
  create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
                          N_("Columns"), "col");

  struct pivot_footnote *footnote = NULL;
  if (log_scale != 0)
    {
      char *s = xasprintf ("× 10**%d", log_scale);
      footnote = pivot_table_create_footnote (
        table, pivot_value_new_user_text_nocopy (s));
    }

  double scale = pow (10.0, log_scale);
  bool numeric = fmt_is_numeric (format.type);
  for (size_t y = 0; y < m->size1; y++)
    for (size_t x = 0; x < m->size2; x++)
      {
        double f = gsl_matrix_get (m, y, x);
        struct pivot_value *v;
        if (numeric)
          {
            v = pivot_value_new_number (f / scale);
            v->numeric.format = format;
          }
        else
          v = pivot_value_new_user_text_nocopy (trimmed_string (f));
        if (footnote)
          pivot_value_add_footnote (v, footnote);
        pivot_table_put2 (table, y, x, v);
      }

  pivot_table_submit (table);
}

static void
matrix_cmd_execute_print (const struct print_command *print)
{
  if (print->expression)
    {
      gsl_matrix *m = matrix_expr_evaluate (print->expression);
      if (!m)
        return;

      int log_scale = 0;
      struct fmt_spec format = (print->use_default_format
                                ? default_format (m, &log_scale)
                                : print->format);

      if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
        matrix_cmd_print_text (print, m, format, log_scale);
      else
        matrix_cmd_print_tables (print, m, format, log_scale);

      gsl_matrix_free (m);
    }
  else
    {
      matrix_cmd_print_space (print->space);
      if (print->title)
        output_log ("%s", print->title);
    }
}
\f
/* DO IF. */

static bool
matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
                       const char *command_name,
                       const char *stop1, const char *stop2)
{
  lex_end_of_command (s->lexer);
  lex_discard_rest_of_command (s->lexer);

  size_t allocated = 0;
  for (;;)
    {
      while (lex_token (s->lexer) == T_ENDCMD)
        lex_get (s->lexer);

      if (lex_at_phrase (s->lexer, stop1)
          || (stop2 && lex_at_phrase (s->lexer, stop2)))
        return true;

      if (lex_at_phrase (s->lexer, "END MATRIX"))
        {
          msg (SE, _("Premature END MATRIX within %s."), command_name);
          return false;
        }

      struct matrix_cmd *cmd = matrix_parse_command (s);
      if (!cmd)
        return false;

      if (c->n >= allocated)
        c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
      c->commands[c->n++] = cmd;
    }
}

static bool
matrix_parse_do_if_clause (struct matrix_state *s,
                           struct do_if_command *ifc,
                           bool parse_condition,
                           size_t *allocated_clauses)
{
  if (ifc->n_clauses >= *allocated_clauses)
    ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
                               sizeof *ifc->clauses);
  struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
  *c = (struct do_if_clause) { .condition = NULL };

  if (parse_condition)
    {
      c->condition = matrix_parse_expr (s);
      if (!c->condition)
        return false;
    }

  return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
}

static struct matrix_cmd *
matrix_parse_do_if (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_DO_IF,
    .do_if = { .n_clauses = 0 }
  };

  size_t allocated_clauses = 0;
  do
    {
      if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
        goto error;
    }
  while (lex_match_phrase (s->lexer, "ELSE IF"));

  if (lex_match_id (s->lexer, "ELSE")
      && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
    goto error;

  if (!lex_match_phrase (s->lexer, "END IF"))
    NOT_REACHED ();
  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static bool
matrix_cmd_execute_do_if (struct do_if_command *cmd)
{
  for (size_t i = 0; i < cmd->n_clauses; i++)
    {
      struct do_if_clause *c = &cmd->clauses[i];
      if (c->condition)
        {
          double d;
          if (!matrix_expr_evaluate_scalar (c->condition,
                                            i ? "ELSE IF" : "DO IF",
                                            &d) || d <= 0)
            continue;
        }

      for (size_t j = 0; j < c->commands.n; j++)
        if (!matrix_cmd_execute (c->commands.commands[j]))
          return false;
      return true;
    }
  return true;
}
\f
static struct matrix_cmd *
matrix_parse_loop (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };

  struct loop_command *loop = &cmd->loop;
  if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
    {
      struct substring name = lex_tokss (s->lexer);
      loop->var = matrix_var_lookup (s, name);
      if (!loop->var)
        loop->var = matrix_var_create (s, name);

      lex_get (s->lexer);
      lex_get (s->lexer);

      loop->start = matrix_parse_expr (s);
      if (!loop->start || !lex_force_match (s->lexer, T_TO))
        goto error;

      loop->end = matrix_parse_expr (s);
      if (!loop->end)
        goto error;

      if (lex_match (s->lexer, T_BY))
        {
          loop->increment = matrix_parse_expr (s);
          if (!loop->increment)
            goto error;
        }
    }

  if (lex_match_id (s->lexer, "IF"))
    {
      loop->top_condition = matrix_parse_expr (s);
      if (!loop->top_condition)
        goto error;
    }

  bool was_in_loop = s->in_loop;
  s->in_loop = true;
  bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
                                   "END LOOP", NULL);
  s->in_loop = was_in_loop;
  if (!ok)
    goto error;

  if (!lex_match_phrase (s->lexer, "END LOOP"))
    NOT_REACHED ();

  if (lex_match_id (s->lexer, "IF"))
    {
      loop->bottom_condition = matrix_parse_expr (s);
      if (!loop->bottom_condition)
        goto error;
    }

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_loop (struct loop_command *cmd)
{
  long int value = 0;
  long int end = 0;
  long int increment = 1;
  if (cmd->var)
    {
      if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
          || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
          || (cmd->increment
              && !matrix_expr_evaluate_integer (cmd->increment, "BY",
                                                &increment)))
        return;

      if (increment > 0 ? value > end
          : increment < 0 ? value < end
          : true)
        return;
    }

  int mxloops = settings_get_mxloops ();
  for (int i = 0; i < mxloops; i++)
    {
      if (cmd->var)
        {
          if (cmd->var->value
              && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
            {
              gsl_matrix_free (cmd->var->value);
              cmd->var->value = NULL;
            }
          if (!cmd->var->value)
            cmd->var->value = gsl_matrix_alloc (1, 1);

          gsl_matrix_set (cmd->var->value, 0, 0, value);
        }

      if (cmd->top_condition)
        {
          double d;
          if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
                                            &d) || d <= 0)
            return;
        }

      for (size_t j = 0; j < cmd->commands.n; j++)
        if (!matrix_cmd_execute (cmd->commands.commands[j]))
          return;

      if (cmd->bottom_condition)
        {
          double d;
          if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
                                            "END LOOP IF", &d) || d > 0)
            return;
        }

      if (cmd->var)
        {
          value += increment;
          if (increment > 0 ? value > end : value < end)
            return;
        }
    }
}
\f
static struct matrix_cmd *
matrix_parse_break (struct matrix_state *s)
{
  if (!s->in_loop)
    {
      msg (SE, _("BREAK not inside LOOP."));
      return NULL;
    }

  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
  return cmd;
}
\f
static struct matrix_cmd *
matrix_parse_display (struct matrix_state *s)
{
  bool dictionary = false;
  bool status = false;
  while (lex_token (s->lexer) == T_ID)
    {
      if (lex_match_id (s->lexer, "DICTIONARY"))
        dictionary = true;
      else if (lex_match_id (s->lexer, "STATUS"))
        status = true;
      else
        {
          lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
          return NULL;
        }
    }
  if (!dictionary && !status)
    dictionary = status = true;

  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_DISPLAY,
    .display = {
      .state = s,
      .dictionary = dictionary,
      .status = status,
    }
  };
  return cmd;
}

static int
compare_matrix_var_pointers (const void *a_, const void *b_)
{
  const struct matrix_var *const *ap = a_;
  const struct matrix_var *const *bp = b_;
  const struct matrix_var *a = *ap;
  const struct matrix_var *b = *bp;
  return strcmp (a->name, b->name);
}

static void
matrix_cmd_execute_display (struct display_command *cmd)
{
  const struct matrix_state *s = cmd->state;

  struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
  pivot_dimension_create (
    table, PIVOT_AXIS_COLUMN, N_("Property"),
    N_("Rows"), N_("Columns"), N_("Size (kB)"));

  const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
  size_t n_vars = 0;

  const struct matrix_var *var;
  HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
    if (var->value)
      vars[n_vars++] = var;
  qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);

  struct pivot_dimension *rows = pivot_dimension_create (
    table, PIVOT_AXIS_ROW, N_("Variable"));
  for (size_t i = 0; i < n_vars; i++)
    {
      const struct matrix_var *var = vars[i];
      pivot_category_create_leaf (
        rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));

      size_t r = var->value->size1;
      size_t c = var->value->size2;
      double values[] = { r, c, r * c * sizeof (double) / 1024 };
      for (size_t j = 0; j < sizeof values / sizeof *values; j++)
        pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
    }
  pivot_table_submit (table);
}
\f
static struct matrix_cmd *
matrix_parse_release (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };

  size_t allocated_vars = 0;
  while (lex_token (s->lexer) == T_ID)
    {
      struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
      if (var)
        {
          if (cmd->release.n_vars >= allocated_vars)
            cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
                                            sizeof *cmd->release.vars);
          cmd->release.vars[cmd->release.n_vars++] = var;
        }
      else
        lex_error (s->lexer, _("Variable name expected."));
      lex_get (s->lexer);

      if (!lex_match (s->lexer, T_COMMA))
        break;
    }

  return cmd;
}

static void
matrix_cmd_execute_release (struct release_command *cmd)
{
  for (size_t i = 0; i < cmd->n_vars; i++)
    {
      struct matrix_var *var = cmd->vars[i];
      gsl_matrix_free (var->value);
      var->value = NULL;
    }
}
\f
static struct matrix_cmd *
matrix_parse_save (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_SAVE,
    .save = {
      .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
    }
  };

  struct save_command *save = &cmd->save;
  save->expression = matrix_parse_exp (s);
  if (!save->expression)
    goto error;

  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "OUTFILE"))
        {
          lex_match (s->lexer, T_EQUALS);

          fh_unref (save->outfile);
          save->outfile = (lex_match (s->lexer, T_ASTERISK)
                           ? fh_inline_file ()
                           : fh_parse (s->lexer, FH_REF_FILE, s->session));
          if (!save->outfile)
            goto error;
        }
      else if (lex_match_id (s->lexer, "VARIABLES"))
        {
          lex_match (s->lexer, T_EQUALS);

          char **names;
          size_t n;
          struct dictionary *d = dict_create (get_default_encoding ());
          bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
                                          PV_NO_SCRATCH | PV_NO_DUPLICATE);
          dict_unref (d);
          if (!ok)
            goto error;

          string_array_destroy (save->variables);
          if (!save->variables)
            save->variables = xmalloc (sizeof *save->variables);
          *save->variables = (struct string_array) {
            .strings = names,
            .n = n,
            .allocated = n,
          };
        }
      else if (lex_match_id (s->lexer, "NAMES"))
        {
          lex_match (s->lexer, T_EQUALS);
          matrix_expr_destroy (save->names);
          save->names = matrix_parse_exp (s);
          if (!save->names)
            goto error;
        }
      else if (lex_match_id (s->lexer, "STRINGS"))
        {
          lex_match (s->lexer, T_EQUALS);
          while (lex_token (s->lexer) == T_ID)
            {
              stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
              lex_get (s->lexer);
              if (!lex_match (s->lexer, T_COMMA))
                break;
            }
        }
      else
        {
          lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
                               "STRINGS");
          goto error;
        }
    }

  if (!save->outfile)
    {
      if (s->prev_save_outfile)
        save->outfile = fh_ref (s->prev_save_outfile);
      else
        {
          lex_sbc_missing ("OUTFILE");
          goto error;
        }
    }
  fh_unref (s->prev_save_outfile);
  s->prev_save_outfile = fh_ref (save->outfile);

  if (save->variables && save->names)
    {
      msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
      matrix_expr_destroy (save->names);
      save->names = NULL;
    }

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_save (const struct save_command *save)
{
  assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */

  gsl_matrix *m = matrix_expr_evaluate (save->expression);
  if (!m)
    return;

  bool ok = true;
  struct dictionary *dict = dict_create (get_default_encoding ());
  struct string_array names = { .n = 0 };
  if (save->variables)
    string_array_clone (&names, save->variables);
  else if (save->names)
    {
      gsl_matrix *nm = matrix_expr_evaluate (save->names);
      if (nm && is_vector (nm))
        {
          gsl_vector nv = to_vector (nm);
          for (size_t i = 0; i < nv.size; i++)
            {
              char *name = trimmed_string (gsl_vector_get (&nv, i));
              if (dict_id_is_valid (dict, name, true))
                string_array_append_nocopy (&names, name);
              else
                ok = false;
            }
        }
    }

  struct stringi_set strings;
  stringi_set_clone (&strings, &save->strings);

  for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
    {
      char tmp_name[64];
      const char *name;
      if (i < names.n)
        name = names.strings[i];
      else
        {
          snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
          name = tmp_name;
        }

      int width = stringi_set_delete (&strings, name) ? 8 : 0;
      struct variable *var = dict_create_var (dict, name, width);
      if (!var)
        {
          msg (SE, _("Duplicate variable name %s in SAVE statement."),
               name);
          ok = false;
        }
    }

  if (!stringi_set_is_empty (&strings))
    {
      const char *example = stringi_set_node_get_string (
        stringi_set_first (&strings));
      msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
                         "with that name was found on SAVE.",
                         "STRINGS specified %2$zu variables, including %1$s, "
                         "whose names were not found on SAVE.",
                         stringi_set_count (&strings)),
           example, stringi_set_count (&strings));
      ok = false;
    }
  stringi_set_destroy (&strings);
  string_array_destroy (&names);

  if (!ok)
    {
      gsl_matrix_free (m);
      dict_unref (dict);
      return;
    }

  struct casewriter *writer = any_writer_open (save->outfile, dict);
  if (!writer)
    {
      gsl_matrix_free (m);
      dict_unref (dict);
      return;
    }

  const struct caseproto *proto = dict_get_proto (dict);
  for (size_t y = 0; y < m->size1; y++)
    {
      struct ccase *c = case_create (proto);
      for (size_t x = 0; x < m->size2; x++)
        {
          double d = gsl_matrix_get (m, y, x);
          int width = caseproto_get_width (proto, x);
          union value *value = case_data_rw_idx (c, x);
          if (width == 0)
            value->f = d;
          else
            memcpy (value->s, &d, width);
        }
      casewriter_write (writer, c);
    }
  casewriter_destroy (writer);

  gsl_matrix_free (m);
  dict_unref (dict);
}
\f
static struct matrix_cmd *
matrix_parse_read (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_READ,
    .read = { .format = FMT_F },
  };

  struct file_handle *fh = NULL;
  char *encoding = NULL;
  struct read_command *read = &cmd->read;
  read->dst = matrix_lvalue_parse (s);
  if (!read->dst)
    goto error;

  int by = 0;
  int repetitions = 0;
  int record_width = 0;
  bool seen_format = false;
  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "FILE"))
        {
          lex_match (s->lexer, T_EQUALS);

          fh_unref (fh);
          fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
          if (!fh)
            goto error;
        }
      else if (lex_match_id (s->lexer, "ENCODING"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (!lex_force_string (s->lexer))
            goto error;

          free (encoding);
          encoding = ss_xstrdup (lex_tokss (s->lexer));

          lex_get (s->lexer);
        }
      else if (lex_match_id (s->lexer, "FIELD"))
        {
          lex_match (s->lexer, T_EQUALS);

          if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
            goto error;
          read->c1 = lex_integer (s->lexer);
          lex_get (s->lexer);
          if (!lex_force_match (s->lexer, T_TO)
              || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
            goto error;
          read->c2 = lex_integer (s->lexer) + 1;
          lex_get (s->lexer);

          record_width = read->c2 - read->c1;
          if (lex_match (s->lexer, T_BY))
            {
              if (!lex_force_int_range (s->lexer, "BY", 1,
                                        read->c2 - read->c1))
                goto error;
              by = lex_integer (s->lexer);
              lex_get (s->lexer);

              if (record_width % by)
                {
                  msg (SE, _("BY %d does not evenly divide record width %d."),
                       by, record_width);
                  goto error;
                }
            }
          else
            by = 0;
        }
      else if (lex_match_id (s->lexer, "SIZE"))
        {
          lex_match (s->lexer, T_EQUALS);
          read->size = matrix_parse_exp (s);
          if (!read->size)
            goto error;
        }
      else if (lex_match_id (s->lexer, "MODE"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (lex_match_id (s->lexer, "RECTANGULAR"))
            read->symmetric = false;
          else if (lex_match_id (s->lexer, "SYMMETRIC"))
            read->symmetric = true;
          else
            {
              lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
              goto error;
            }
        }
      else if (lex_match_id (s->lexer, "REREAD"))
        read->reread = true;
      else if (lex_match_id (s->lexer, "FORMAT"))
        {
          if (seen_format)
            {
              lex_sbc_only_once ("FORMAT");
              goto error;
            }
          seen_format = true;

          lex_match (s->lexer, T_EQUALS);

          if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
            goto error;

          const char *p = lex_tokcstr (s->lexer);
          if (c_isdigit (p[0]))
            {
              repetitions = atoi (p);
              p += strspn (p, "0123456789");
              if (!fmt_from_name (p, &read->format))
                {
                  lex_error (s->lexer, _("Unknown format %s."), p);
                  goto error;
                }
              lex_get (s->lexer);
            }
          else if (!fmt_from_name (p, &read->format))
            {
              struct fmt_spec format;
              if (!parse_format_specifier (s->lexer, &format))
                goto error;
              read->format = format.type;
              read->w = format.w;
              read->d = format.d;
            }
        }
      else
        {
          lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
                               "REREAD", "FORMAT");
          goto error;
        }
    }

  if (!read->c1)
    {
      lex_sbc_missing ("FIELD");
      goto error;
    }

  if (!read->dst->n_indexes && !read->size)
    {
      msg (SE, _("SIZE is required for reading data into a full matrix "
                 "(as opposed to a submatrix)."));
      goto error;
    }

  if (!fh)
    {
      if (s->prev_read_file)
        fh = fh_ref (s->prev_read_file);
      else
        {
          lex_sbc_missing ("FILE");
          goto error;
        }
    }
  fh_unref (s->prev_read_file);
  s->prev_read_file = fh_ref (fh);

  read->rf = read_file_create (s, fh);
  fh = NULL;
  if (encoding)
    {
      free (read->rf->encoding);
      read->rf->encoding = encoding;
      encoding = NULL;
    }

  /* Field width may be specified in multiple ways:

     1. BY on FIELD.
     2. The format on FORMAT.
     3. The repetition factor on FORMAT.

     (2) and (3) are mutually exclusive.

     If more than one of these is present, they must agree.  If none of them is
     present, then free-field format is used.
   */
  if (repetitions > record_width)
    {
      msg (SE, _("%d repetitions cannot fit in record width %d."),
           repetitions, record_width);
      goto error;
    }
  int w = (repetitions ? record_width / repetitions
           : read->w ? read->w
           : by);
  if (by && w != by)
    {
      if (repetitions)
        msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
                   "which implies field width %d, "
                   "but BY specifies field width %d."),
             repetitions, record_width, w, by);
      else
        msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
             w, by);
      goto error;
    }
  read->w = w;
  return cmd;

error:
  fh_unref (fh);
  matrix_cmd_destroy (cmd);
  free (encoding);
  return NULL;
}

static void
matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
                       gsl_matrix *m, struct substring p, size_t y, size_t x)
{
  const char *input_encoding = dfm_reader_get_encoding (reader);
  union value v;
  char *error = data_in (p, input_encoding, read->format,
                         settings_get_fmt_settings (), &v, 0, NULL);
  /* XXX report error if value is missing */
  if (error)
    msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
  else
    {
      gsl_matrix_set (m, y, x, v.f);
      if (read->symmetric && x != y)
        gsl_matrix_set (m, x, y, v.f);
    }
}

static bool
matrix_read_line (struct read_command *read, struct dfm_reader *reader,
                  struct substring *line)
{
  if (dfm_eof (reader))
    {
      msg (SE, _("Unexpected end of file reading matrix data."));
      return false;
    }
  dfm_expand_tabs (reader);
  *line = ss_substr (dfm_get_record (reader),
                     read->c1 - 1, read->c2 - read->c1);
  return true;
}

static void
matrix_read (struct read_command *read, struct dfm_reader *reader,
             gsl_matrix *m)
{
  for (size_t y = 0; y < m->size1; y++)
    {
      size_t nx = read->symmetric ? y + 1 : m->size2;

      struct substring line = ss_empty ();
      for (size_t x = 0; x < nx; x++)
        {
          struct substring p;
          if (!read->w)
            {
              for (;;)
                {
                  ss_ltrim (&line, ss_cstr (" ,"));
                  if (!ss_is_empty (line))
                    break;
                  if (!matrix_read_line (read, reader, &line))
                    return;
                  dfm_forward_record (reader);
                }

              ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
            }
          else
            {
              if (!matrix_read_line (read, reader, &line))
                return;
              size_t fields_per_line = (read->c2 - read->c1) / read->w;
              int f = x % fields_per_line;
              if (f == fields_per_line - 1)
                dfm_forward_record (reader);

              p = ss_substr (line, read->w * f, read->w);
            }

          matrix_read_set_field (read, reader, m, p, y, x);
        }

      if (read->w)
        dfm_forward_record (reader);
      else
        {
          ss_ltrim (&line, ss_cstr (" ,"));
          if (!ss_is_empty (line))
            msg (SW, _("Trailing garbage on line \"%.*s\""),
                 (int) line.length, line.string);
        }
    }
}

static void
matrix_cmd_execute_read (struct read_command *read)
{
  struct index_vector iv0, iv1;
  if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
    return;

  size_t size[2] = { SIZE_MAX, SIZE_MAX };
  if (read->size)
    {
      gsl_matrix *m = matrix_expr_evaluate (read->size);
      if (!m)
        return;

      if (!is_vector (m))
        {
          msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
                     "not a %zu×%zu matrix."), m->size1, m->size2);
          gsl_matrix_free (m);
          free (iv0.indexes);
          free (iv1.indexes);
          return;
        }

      gsl_vector v = to_vector (m);
      double d[2];
      if (v.size == 1)
        {
          d[0] = gsl_vector_get (&v, 0);
          d[1] = 1;
        }
      else if (v.size == 2)
        {
          d[0] = gsl_vector_get (&v, 0);
          d[1] = gsl_vector_get (&v, 1);
        }
      else
        {
          msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
                     "not a %zu×%zu matrix."),
               m->size1, m->size2),
          gsl_matrix_free (m);
          free (iv0.indexes);
          free (iv1.indexes);
          return;
        }
      gsl_matrix_free (m);

      if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
        {
          msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
          free (iv0.indexes);
          free (iv1.indexes);
          return;
        }

      size[0] = d[0];
      size[1] = d[1];
    }

  if (read->dst->n_indexes)
    {
      size_t submatrix_size[2];
      if (read->dst->n_indexes == 2)
        {
          submatrix_size[0] = iv0.n;
          submatrix_size[1] = iv1.n;
        }
      else if (read->dst->var->value->size1 == 1)
        {
          submatrix_size[0] = 1;
          submatrix_size[1] = iv0.n;
        }
      else
        {
          submatrix_size[0] = iv0.n;
          submatrix_size[1] = 1;
        }

      if (read->size)
        {
          if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
            {
              msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
                         "%zu×%zu."),
                   size[0], size[1],
                   submatrix_size[0], submatrix_size[1]);
              free (iv0.indexes);
              free (iv1.indexes);
              return;
            }
        }
      else
        {
          size[0] = submatrix_size[0];
          size[1] = submatrix_size[1];
        }
    }

  struct dfm_reader *reader = read_file_open (read->rf);
  if (read->reread)
    dfm_reread_record (reader, 1);

  if (read->symmetric && size[0] != size[1])
    {
      msg (SE, _("Cannot read non-square %zu×%zu matrix "
                 "using READ with MODE=SYMMETRIC."),
           size[0], size[1]);
      free (iv0.indexes);
      free (iv1.indexes);
      return;
    }
  gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
  matrix_read (read, reader, tmp);
  matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
}
\f
static struct matrix_cmd *
matrix_parse_write (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_WRITE,
    .write = { .format = FMT_F },
  };

  struct file_handle *fh = NULL;
  char *encoding = NULL;
  struct write_command *write = &cmd->write;
  write->expression = matrix_parse_exp (s);
  if (!write->expression)
    goto error;

  int by = 0;
  int repetitions = 0;
  int record_width = 0;
  bool seen_format = false;
  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "OUTFILE"))
        {
          lex_match (s->lexer, T_EQUALS);

          fh_unref (fh);
          fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
          if (!fh)
            goto error;
        }
      else if (lex_match_id (s->lexer, "ENCODING"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (!lex_force_string (s->lexer))
            goto error;

          free (encoding);
          encoding = ss_xstrdup (lex_tokss (s->lexer));

          lex_get (s->lexer);
        }
      else if (lex_match_id (s->lexer, "FIELD"))
        {
          lex_match (s->lexer, T_EQUALS);

          if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
            goto error;
          write->c1 = lex_integer (s->lexer);
          lex_get (s->lexer);
          if (!lex_force_match (s->lexer, T_TO)
              || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
            goto error;
          write->c2 = lex_integer (s->lexer) + 1;
          lex_get (s->lexer);

          record_width = write->c2 - write->c1;
          if (lex_match (s->lexer, T_BY))
            {
              if (!lex_force_int_range (s->lexer, "BY", 1,
                                        write->c2 - write->c1))
                goto error;
              by = lex_integer (s->lexer);
              lex_get (s->lexer);

              if (record_width % by)
                {
                  msg (SE, _("BY %d does not evenly divide record width %d."),
                       by, record_width);
                  goto error;
                }
            }
          else
            by = 0;
        }
      else if (lex_match_id (s->lexer, "MODE"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (lex_match_id (s->lexer, "RECTANGULAR"))
            write->triangular = false;
          else if (lex_match_id (s->lexer, "TRIANGULAR"))
            write->triangular = true;
          else
            {
              lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
              goto error;
            }
        }
      else if (lex_match_id (s->lexer, "HOLD"))
        write->hold = true;
      else if (lex_match_id (s->lexer, "FORMAT"))
        {
          if (seen_format)
            {
              lex_sbc_only_once ("FORMAT");
              goto error;
            }
          seen_format = true;

          lex_match (s->lexer, T_EQUALS);

          if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
            goto error;

          const char *p = lex_tokcstr (s->lexer);
          if (c_isdigit (p[0]))
            {
              repetitions = atoi (p);
              p += strspn (p, "0123456789");
              if (!fmt_from_name (p, &write->format))
                {
                  lex_error (s->lexer, _("Unknown format %s."), p);
                  goto error;
                }
              lex_get (s->lexer);
            }
          else if (!fmt_from_name (p, &write->format))
            {
              struct fmt_spec format;
              if (!parse_format_specifier (s->lexer, &format))
                goto error;
              write->format = format.type;
              write->w = format.w;
              write->d = format.d;
            }
        }
      else
        {
          lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
                               "HOLD", "FORMAT");
          goto error;
        }
    }

  if (!write->c1)
    {
      lex_sbc_missing ("FIELD");
      goto error;
    }

  if (!fh)
    {
      if (s->prev_write_file)
        fh = fh_ref (s->prev_write_file);
      else
        {
          lex_sbc_missing ("OUTFILE");
          goto error;
        }
    }
  fh_unref (s->prev_write_file);
  s->prev_write_file = fh_ref (fh);

  write->wf = write_file_create (s, fh);
  fh = NULL;
  if (encoding)
    {
      free (write->wf->encoding);
      write->wf->encoding = encoding;
      encoding = NULL;
    }

  /* Field width may be specified in multiple ways:

     1. BY on FIELD.
     2. The format on FORMAT.
     3. The repetition factor on FORMAT.

     (2) and (3) are mutually exclusive.

     If more than one of these is present, they must agree.  If none of them is
     present, then free-field format is used.
   */
  if (repetitions > record_width)
    {
      msg (SE, _("%d repetitions cannot fit in record width %d."),
           repetitions, record_width);
      goto error;
    }
  int w = (repetitions ? record_width / repetitions
           : write->w ? write->w
           : by);
  if (by && w != by)
    {
      if (repetitions)
        msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
                   "which implies field width %d, "
                   "but BY specifies field width %d."),
             repetitions, record_width, w, by);
      else
        msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
             w, by);
      goto error;
    }
  write->w = w;
  return cmd;

error:
  fh_unref (fh);
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_write (struct write_command *write)
{
  gsl_matrix *m = matrix_expr_evaluate (write->expression);
  if (!m)
    return;

  if (write->triangular && m->size1 != m->size2)
    {
      msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
                 "the matrix to be written has dimensions %zu×%zu."),
           m->size1, m->size2);
      gsl_matrix_free (m);
      return;
    }

  struct dfm_writer *writer = write_file_open (write->wf);
  if (!writer || !m->size1)
    {
      gsl_matrix_free (m);
      return;
    }

  const struct fmt_settings *settings = settings_get_fmt_settings ();
  struct fmt_spec format = {
    .type = write->format,
    .w = write->w ? write->w : 40,
    .d = write->d
  };
  struct u8_line *line = write->wf->held;
  for (size_t y = 0; y < m->size1; y++)
    {
      if (!line)
        {
          line = xmalloc (sizeof *line);
          u8_line_init (line);
        }
      size_t nx = write->triangular ? y + 1 : m->size2;
      int x0 = write->c1;
      for (size_t x = 0; x < nx; x++)
        {
          /* XXX string values */
          union value v = { .f = gsl_matrix_get (m, y, x) };
          char *s = (write->w
                     ? data_out (&v, NULL, &format, settings)
                     : data_out_stretchy (&v, NULL, &format, settings, NULL));
          size_t len = strlen (s);
          int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
          if (width + x0 > write->c2)
            {
              dfm_put_record_utf8 (writer, line->s.ss.string,
                                   line->s.ss.length);
              u8_line_clear (line);
              x0 = write->c1;
            }
          u8_line_put (line, x0, x0 + width, s, len);
          free (s);

          x0 += write->w ? write->w : width + 1;
        }

      if (y + 1 >= m->size1 && write->hold)
        break;
      dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
      u8_line_clear (line);
    }
  if (!write->hold)
    {
      u8_line_destroy (line);
      line = NULL;
    }
  write->wf->held = line;
  dfm_close_writer (writer);

  gsl_matrix_free (m);
}
\f
static struct matrix_cmd *
matrix_parse_get (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_GET,
    .get = {
      .user = { .treatment = MGET_ERROR },
      .system = { .treatment = MGET_ERROR },
    }
  };

  struct get_command *get = &cmd->get;
  get->dst = matrix_lvalue_parse (s);
  if (!get->dst)
    goto error;

  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "FILE"))
        {
          if (get->variables.n)
            {
              lex_error (s->lexer, _("FILE must precede VARIABLES"));
              goto error;
            }
          lex_match (s->lexer, T_EQUALS);

          fh_unref (get->file);
          get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
          if (!get->file)
            goto error;
        }
      else if (lex_match_id (s->lexer, "ENCODING"))
        {
          if (get->variables.n)
            {
              lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
              goto error;
            }
          lex_match (s->lexer, T_EQUALS);
          if (!lex_force_string (s->lexer))
            goto error;

          free (get->encoding);
          get->encoding = ss_xstrdup (lex_tokss (s->lexer));

          lex_get (s->lexer);
        }
      else if (lex_match_id (s->lexer, "VARIABLES"))
        {
          lex_match (s->lexer, T_EQUALS);

          struct dictionary *dict = NULL;
          if (!get->file)
            {
              dict = dataset_dict (s->dataset);
              if (dict_get_var_cnt (dict) == 0)
                {
                  lex_error (s->lexer, _("GET cannot read empty active file."));
                  goto error;
                }
            }
          else
            {
              struct casereader *reader = any_reader_open_and_decode (
                get->file, get->encoding, &dict, NULL);
              if (!reader)
                goto error;
              casereader_destroy (reader);
            }

          struct variable **vars;
          size_t n_vars;
          bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
                                     PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
          if (!ok)
            {
              dict_unref (dict);
              goto error;
            }

          string_array_clear (&get->variables);
          for (size_t i = 0; i < n_vars; i++)
            string_array_append (&get->variables, var_get_name (vars[i]));
          free (vars);
          dict_unref (dict);
        }
      else if (lex_match_id (s->lexer, "NAMES"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (!lex_force_id (s->lexer))
            goto error;

          struct substring name = lex_tokss (s->lexer);
          get->names = matrix_var_lookup (s, name);
          if (!get->names)
            get->names = matrix_var_create (s, name);
          lex_get (s->lexer);
        }
      else if (lex_match_id (s->lexer, "MISSING"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (lex_match_id (s->lexer, "ACCEPT"))
            get->user.treatment = MGET_ACCEPT;
          else if (lex_match_id (s->lexer, "OMIT"))
            get->user.treatment = MGET_OMIT;
          else if (lex_is_number (s->lexer))
            {
              get->user.treatment = MGET_RECODE;
              get->user.substitute = lex_number (s->lexer);
              lex_get (s->lexer);
            }
          else
            {
              lex_error (s->lexer, NULL);
              goto error;
            }
        }
      else if (lex_match_id (s->lexer, "SYSMIS"))
        {
          lex_match (s->lexer, T_EQUALS);
          if (lex_match_id (s->lexer, "OMIT"))
            get->user.treatment = MGET_OMIT;
          else if (lex_is_number (s->lexer))
            {
              get->user.treatment = MGET_RECODE;
              get->user.substitute = lex_number (s->lexer);
              lex_get (s->lexer);
            }
          else
            {
              lex_error (s->lexer, NULL);
              goto error;
            }
        }
      else
        {
          lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
                               "MISSING", "SYSMIS");
          goto error;
        }
    }
  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_get (struct get_command *get)
{
  assert (get->file);           /* XXX */

  struct dictionary *dict;
  struct casereader *reader = any_reader_open_and_decode (
    get->file, get->encoding, &dict, NULL);
  if (!reader)
    return;

  const struct variable **vars = xnmalloc (
    get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
    sizeof *vars);
  size_t n_vars = 0;

  if (get->variables.n)
    {
      for (size_t i = 0; i < get->variables.n; i++)
        {
          const char *name = get->variables.strings[i];
          const struct variable *var = dict_lookup_var (dict, name);
          if (!var)
            {
              msg (SE, _("GET: Data file does not contain variable %s."),
                   name);
              dict_unref (dict);
              free (vars);
              return;
            }
          if (!var_is_numeric (var))
            {
              msg (SE, _("GET: Variable %s is not numeric."), name);
              dict_unref (dict);
              free (vars);
              return;
            }
          vars[n_vars++] = var;
        }
    }
  else
    {
      for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
        {
          const struct variable *var = dict_get_var (dict, i);
          if (!var_is_numeric (var))
            {
              msg (SE, _("GET: Variable %s is not numeric."),
                   var_get_name (var));
              dict_unref (dict);
              free (vars);
              return;
            }
          vars[n_vars++] = var;
        }
    }

  size_t n_rows = 0;
  gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
  long long int casenum = 1;
  bool error = false;
  for (struct ccase *c = casereader_read (reader); c;
       c = casereader_read (reader), casenum++)
    {
      if (n_rows >= m->size1)
        {
          gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
          for (size_t y = 0; y < n_rows; y++)
            for (size_t x = 0; x < n_vars; x++)
              gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
          gsl_matrix_free (m);
          m = p;
        }

      bool keep = true;
      for (size_t x = 0; x < n_vars; x++)
        {
          const struct variable *var = vars[x];
          double d = case_num (c, var);
          if (d == SYSMIS)
            {
              if (get->system.treatment == MGET_RECODE)
                d = get->system.substitute;
              else if (get->system.treatment == MGET_OMIT)
                keep = false;
              else
                {
                  msg (SE, _("GET: Variable %s in case %lld "
                             "is system-missing."),
                       var_get_name (var), casenum);
                  error = true;
                }
            }
          else if (var_is_num_missing (var, d, MV_USER))
            {
              if (get->user.treatment == MGET_RECODE)
                d = get->user.substitute;
              else if (get->user.treatment == MGET_OMIT)
                keep = false;
              else if (get->user.treatment != MGET_ACCEPT)
                {
                  msg (SE, _("GET: Variable %s in case %lld has user-missing "
                             "value %g."),
                       var_get_name (var), casenum, d);
                  error = true;
                }
            }
          gsl_matrix_set (m, n_rows, x, d);
        }
      case_unref (c);
      if (error)
        break;
      if (keep)
        n_rows++;
    }
  casereader_destroy (reader);
  if (!error)
    {
      m->size1 = n_rows;
      matrix_lvalue_evaluate_and_assign (get->dst, m);
    }
  else
    gsl_matrix_free (m);
  dict_unref (dict);
  free (vars);
}
\f
static const char *
match_rowtype (struct lexer *lexer)
{
  static const char *rowtypes[] = {
    "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
  };
  size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;

  for (size_t i = 0; i < n_rowtypes; i++)
    if (lex_match_id (lexer, rowtypes[i]))
      return rowtypes[i];

  lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
  return NULL;
}

static bool
parse_var_names (struct lexer *lexer, struct string_array *sa)
{
  lex_match (lexer, T_EQUALS);

  string_array_clear (sa);

  struct dictionary *dict = dict_create (get_default_encoding ());
  dict_create_var_assert (dict, "ROWTYPE_", 8);
  dict_create_var_assert (dict, "VARNAME_", 8);
  char **names;
  size_t n_names;
  bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
                                  PV_NO_DUPLICATE | PV_NO_SCRATCH);
  dict_unref (dict);

  if (ok)
    {
      string_array_clear (sa);
      sa->strings = names;
      sa->n = sa->allocated = n_names;
    }
  return ok;
}

static void
msave_common_uninit (struct msave_common *common)
{
  if (common)
    {
      fh_unref (common->outfile);
      string_array_destroy (&common->variables);
      string_array_destroy (&common->fnames);
      string_array_destroy (&common->snames);
    }
}

static bool
compare_variables (const char *keyword,
                   const struct string_array *new,
                   const struct string_array *old)
{
  if (new->n)
    {
      if (!old->n)
        {
          msg (SE, _("%s may only be specified on MSAVE if it was specified "
                     "on the first MSAVE within MATRIX."), keyword);
          return false;
        }
      else if (!string_array_equal_case (old, new))
        {
          msg (SE, _("%s must specify the same variables each time within "
                     "a given MATRIX."), keyword);
          return false;
        }
    }
  return true;
}

static struct matrix_cmd *
matrix_parse_msave (struct matrix_state *s)
{
  struct msave_common common = { .outfile = NULL };
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };

  struct msave_command *msave = &cmd->msave;
  if (lex_token (s->lexer) == T_ID)
    msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
  msave->expr = matrix_parse_exp (s);
  if (!msave->expr)
    return NULL;

  while (lex_match (s->lexer, T_SLASH))
    {
      if (lex_match_id (s->lexer, "TYPE"))
        {
          lex_match (s->lexer, T_EQUALS);

          msave->rowtype = match_rowtype (s->lexer);
          if (!msave->rowtype)
            goto error;
        }
      else if (lex_match_id (s->lexer, "OUTFILE"))
        {
          lex_match (s->lexer, T_EQUALS);

          fh_unref (common.outfile);
          common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
          if (!common.outfile)
            goto error;
        }
      else if (lex_match_id (s->lexer, "VARIABLES"))
        {
          if (!parse_var_names (s->lexer, &common.variables))
            goto error;
        }
      else if (lex_match_id (s->lexer, "FNAMES"))
        {
          if (!parse_var_names (s->lexer, &common.fnames))
            goto error;
        }
      else if (lex_match_id (s->lexer, "SNAMES"))
        {
          if (!parse_var_names (s->lexer, &common.snames))
            goto error;
        }
      else if (lex_match_id (s->lexer, "SPLIT"))
        {
          lex_match (s->lexer, T_EQUALS);

          matrix_expr_destroy (msave->splits);
          msave->splits = matrix_parse_exp (s);
          if (!msave->splits)
            goto error;
        }
      else if (lex_match_id (s->lexer, "FACTOR"))
        {
          lex_match (s->lexer, T_EQUALS);

          matrix_expr_destroy (msave->factors);
          msave->factors = matrix_parse_exp (s);
          if (!msave->factors)
            goto error;
        }
      else
        {
          lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
                               "FNAMES", "SNAMES", "SPLIT", "FACTOR");
          goto error;
        }
    }
  if (!msave->rowtype)
    {
      lex_sbc_missing ("TYPE");
      goto error;
    }
  common.has_splits = msave->splits || common.snames.n;
  common.has_factors = msave->factors || common.fnames.n;

  struct msave_common *c = s->common ? s->common : &common;
  if (c->fnames.n && !msave->factors)
    {
      msg (SE, _("FNAMES requires FACTOR."));
      goto error;
    }
  if (c->snames.n && !msave->splits)
    {
      msg (SE, _("SNAMES requires SPLIT."));
      goto error;
    }
  if (c->has_factors && !common.has_factors)
    {
      msg (SE, _("%s is required because it was present on the first "
                 "MSAVE in this MATRIX command."), "FACTOR");
      goto error;
    }
  if (c->has_splits && !common.has_splits)
    {
      msg (SE, _("%s is required because it was present on the first "
                 "MSAVE in this MATRIX command."), "SPLIT");
      goto error;
    }

  if (!s->common)
    {
      if (!common.outfile)
        {
          lex_sbc_missing ("OUTFILE");
          goto error;
        }
      s->common = xmemdup (&common, sizeof common);
    }
  else
    {
      if (common.outfile)
        {
          bool same = common.outfile == s->common->outfile;
          fh_unref (common.outfile);
          if (!same)
            {
              msg (SE, _("OUTFILE must name the same file on each MSAVE "
                         "within a single MATRIX command."));
              goto error;
            }
        }
      if (!compare_variables ("VARIABLES",
                              &common.variables, &s->common->variables)
          || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
          || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
        goto error;
      msave_common_uninit (&common);
    }
  msave->common = s->common;
  if (!msave->varname_)
    msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
  return cmd;

error:
  msave_common_uninit (&common);
  matrix_cmd_destroy (cmd);
  return NULL;
}

static gsl_vector *
matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
{
  gsl_matrix *m = matrix_expr_evaluate (e);
  if (!m)
    return NULL;

  if (!is_vector (m))
    {
      msg (SE, _("%s expression must evaluate to vector, "
                 "not a %zu×%zu matrix."),
           name, m->size1, m->size2);
      gsl_matrix_free (m);
      return NULL;
    }

  return matrix_to_vector (m);
}

static const char *
msave_add_vars (struct dictionary *d, const struct string_array *vars)
{
  for (size_t i = 0; i < vars->n; i++)
    if (!dict_create_var (d, vars->strings[i], 0))
      return vars->strings[i];
  return NULL;
}

static struct dictionary *
msave_create_dict (const struct msave_common *common)
{
  struct dictionary *dict = dict_create (get_default_encoding ());

  const char *dup_split = msave_add_vars (dict, &common->snames);
  if (dup_split)
    {
      msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
      goto error;
    }

  dict_create_var_assert (dict, "ROWTYPE_", 8);

  const char *dup_factor = msave_add_vars (dict, &common->fnames);
  if (dup_factor)
    {
      msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
      goto error;
    }

  dict_create_var_assert (dict, "VARNAME_", 8);

  const char *dup_var = msave_add_vars (dict, &common->variables);
  if (dup_var)
    {
      msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
      goto error;
    }

  return dict;

error:
  dict_unref (dict);
  return NULL;
}

static void
matrix_cmd_execute_msave (struct msave_command *msave)
{
  struct msave_common *common = msave->common;
  gsl_matrix *m = NULL;
  gsl_vector *factors = NULL;
  gsl_vector *splits = NULL;

  m = matrix_expr_evaluate (msave->expr);
  if (!m)
    goto error;

  if (!common->variables.n)
    for (size_t i = 0; i < m->size2; i++)
      string_array_append_nocopy (&common->variables,
                                  xasprintf ("COL%zu", i + 1));

  if (m->size2 != common->variables.n)
    {
      msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
           m->size2, common->variables.n);
      goto error;
    }

  if (msave->factors)
    {
      factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
      if (!factors)
        goto error;

      if (!common->fnames.n)
        for (size_t i = 0; i < factors->size; i++)
          string_array_append_nocopy (&common->fnames,
                                      xasprintf ("FAC%zu", i + 1));

      if (factors->size != common->fnames.n)
        {
          msg (SE, _("There are %zu factor variables, "
                     "but %zu split values were supplied."),
               common->fnames.n, factors->size);
          goto error;
        }
    }

  if (msave->splits)
    {
      splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
      if (!splits)
        goto error;

      if (!common->fnames.n)
        for (size_t i = 0; i < splits->size; i++)
          string_array_append_nocopy (&common->fnames,
                                      xasprintf ("SPL%zu", i + 1));

      if (splits->size != common->fnames.n)
        {
          msg (SE, _("There are %zu split variables, "
                     "but %zu split values were supplied."),
               common->fnames.n, splits->size);
          goto error;
        }
    }

  if (!common->writer)
    {
      struct dictionary *dict = msave_create_dict (common);
      if (!dict)
        goto error;

      common->writer = any_writer_open (common->outfile, dict);
      if (!common->writer)
        {
          dict_unref (dict);
          goto error;
        }

      common->dict = dict;
    }

  for (size_t y = 0; y < m->size1; y++)
    {
      struct ccase *c = case_create (dict_get_proto (common->dict));
      size_t idx = 0;

      /* Split variables */
      if (splits)
        for (size_t i = 0; i < splits->size; i++)
          case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);

      /* ROWTYPE_. */
      buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
                         msave->rowtype, ' ');

      /* Factors. */
      if (factors)
        for (size_t i = 0; i < factors->size; i++)
          *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);

      /* VARNAME_. */
      buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
                         msave->varname_, ' ');

      /* Continuous variables. */
      for (size_t x = 0; x < m->size2; x++)
        case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
      casewriter_write (common->writer, c);
    }

error:
  gsl_matrix_free (m);
  gsl_vector_free (factors);
  gsl_vector_free (splits);
}
\f
static struct matrix_cmd *
matrix_parse_mget (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };

  struct mget_command *mget = &cmd->mget;

  for (;;)
    {
      if (lex_match_id (s->lexer, "FILE"))
        {
          lex_match (s->lexer, T_EQUALS);

          fh_unref (mget->file);
          mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
          if (!mget->file)
            goto error;
        }
      else if (lex_match_id (s->lexer, "TYPE"))
        {
          lex_match (s->lexer, T_EQUALS);
          while (lex_token (s->lexer) != T_SLASH
                 && lex_token (s->lexer) != T_ENDCMD)
            {
              const char *rowtype = match_rowtype (s->lexer);
              if (!rowtype)
                goto error;

              stringi_set_insert (&mget->rowtypes, rowtype);
            }
        }
      else
        {
          lex_error_expecting (s->lexer, "FILE", "TYPE");
          goto error;
        }
      if (lex_token (s->lexer) == T_ENDCMD)
        break;

      if (!lex_force_match (s->lexer, T_SLASH))
        goto error;
    }
  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static const struct variable *
get_a8_var (const struct dictionary *d, const char *name)
{
  const struct variable *v = dict_lookup_var (d, name);
  if (!v)
    {
      msg (SE, _("Matrix data file lacks %s variable."), name);
      return NULL;
    }
  if (var_get_width (v) != 8)
    {
      msg (SE, _("%s variable in matrix data file must be "
                 "8-byte string, but it has width %d."),
           name, var_get_width (v));
      return NULL;
    }
  return v;
}

static bool
is_rowtype (const union value *v, const char *rowtype)
{
  struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
  ss_rtrim (&vs, ss_cstr (" "));
  return ss_equals_case (vs, ss_cstr (rowtype));
}

static void
matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
                        const struct dictionary *d,
                        const struct variable *rowtype_var,
                        struct matrix_state *s, size_t si, size_t fi,
                        size_t cs, size_t cn)
{
  if (!n_rows)
    return;

  const union value *rowtype_ = case_data (rows[0], rowtype_var);
  const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
                             : is_rowtype (rowtype_, "CORR") ? "CR"
                             : is_rowtype (rowtype_, "MEAN") ? "MN"
                             : is_rowtype (rowtype_, "STDDEV") ? "SD"
                             : is_rowtype (rowtype_, "N") ? "NC"
                             : "CN");

  struct string name = DS_EMPTY_INITIALIZER;
  ds_put_cstr (&name, name_prefix);
  if (fi > 0)
    ds_put_format (&name, "F%zu", fi);
  if (si > 0)
    ds_put_format (&name, "S%zu", si);

  struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
  if (!mv)
    mv = matrix_var_create (s, ds_ss (&name));
  else if (mv->value)
    {
      msg (SW, _("Matrix data file contains variable with existing name %s."),
           ds_cstr (&name));
      goto exit;
    }

  gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
  size_t n_missing = 0;
  for (size_t y = 0; y < n_rows; y++)
    {
      for (size_t x = 0; x < cn; x++)
        {
          struct variable *var = dict_get_var (d, cs + x);
          double value = case_num (rows[y], var);
          if (var_is_num_missing (var, value, MV_ANY))
            {
              n_missing++;
              value = 0.0;
            }
          gsl_matrix_set (m, y, x, value);
        }
    }

  if (n_missing)
    msg (SE, ngettext ("Matrix data file variable %s contains a missing "
                       "value, which was treated as zero.",
                       "Matrix data file variable %s contains %zu missing "
                       "values, which were treated as zero.", n_missing),
         ds_cstr (&name), n_missing);
  mv->value = m;

exit:
  ds_destroy (&name);
  for (size_t y = 0; y < n_rows; y++)
    case_unref (rows[y]);
}

static bool
var_changed (const struct ccase *ca, const struct ccase *cb,
             const struct variable *var)
{
  return !value_equal (case_data (ca, var), case_data (cb, var),
                       var_get_width (var));
}

static bool
vars_changed (const struct ccase *ca, const struct ccase *cb,
              const struct dictionary *d,
              size_t first_var, size_t n_vars)
{
  for (size_t i = 0; i < n_vars; i++)
    {
      const struct variable *v = dict_get_var (d, first_var + i);
      if (var_changed (ca, cb, v))
        return true;
    }
  return false;
}

static void
matrix_cmd_execute_mget (struct mget_command *mget)
{
  struct dictionary *d;
  struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
                                                     &d, NULL);
  if (!r)
    return;

  const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
  const struct variable *varname_ = get_a8_var (d, "VARNAME_");
  if (!rowtype_ || !varname_)
    goto exit;

  if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
    {
      msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
      goto exit;
    }
  if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
    {
      msg (SE, _("Matrix data file contains no continuous variables."));
      goto exit;
    }

  for (size_t i = 0; i < dict_get_var_cnt (d); i++)
    {
      const struct variable *v = dict_get_var (d, i);
      if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
        {
          msg (SE,
               _("Matrix data file contains unexpected string variable %s."),
               var_get_name (v));
          goto exit;
        }
    }

  /* SPLIT variables. */
  size_t ss = 0;
  size_t sn = var_get_dict_index (rowtype_);
  struct ccase *sc = NULL;
  size_t si = 0;

  /* FACTOR variables. */
  size_t fs = var_get_dict_index (rowtype_) + 1;
  size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
  struct ccase *fc = NULL;
  size_t fi = 0;

  /* Continuous variables. */
  size_t cs = var_get_dict_index (varname_) + 1;
  size_t cn = dict_get_var_cnt (d) - cs;
  struct ccase *cc = NULL;

  /* Matrix. */
  struct ccase **rows = NULL;
  size_t allocated_rows = 0;
  size_t n_rows = 0;

  struct ccase *c;
  while ((c = casereader_read (r)) != NULL)
    {
      bool sd = vars_changed (sc, c, d, ss, sn);
      bool fd = sd || vars_changed (fc, c, d, fs, fn);
      bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
      if (sd)
        {
          si++;
          case_unref (sc);
          sc = case_ref (c);
        }
      if (fd)
        {
          fi++;
          case_unref (fc);
          fc = case_ref (c);
        }
      if (md)
        {
          matrix_mget_commit_var (rows, n_rows, d, rowtype_,
                                  mget->state, si, fi, cs, cn);
          n_rows = 0;
          case_unref (cc);
          cc = case_ref (c);
        }

      if (n_rows >= allocated_rows)
        rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
      rows[n_rows++] = c;
    }
  matrix_mget_commit_var (rows, n_rows, d, rowtype_,
                          mget->state, si, fi, cs, cn);
  free (rows);

exit:
  casereader_destroy (r);
}
\f
static bool
matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
{
  if (!lex_force_id (s->lexer))
    return false;

  *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
  if (!*varp)
    *varp = matrix_var_create (s, lex_tokss (s->lexer));
  lex_get (s->lexer);
  return true;
}

static struct matrix_cmd *
matrix_parse_eigen (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_EIGEN,
    .eigen = { .expr = NULL }
  };

  struct eigen_command *eigen = &cmd->eigen;
  if (!lex_force_match (s->lexer, T_LPAREN))
    goto error;
  eigen->expr = matrix_parse_expr (s);
  if (!eigen->expr
      || !lex_force_match (s->lexer, T_COMMA)
      || !matrix_parse_dst_var (s, &eigen->evec)
      || !lex_force_match (s->lexer, T_COMMA)
      || !matrix_parse_dst_var (s, &eigen->eval)
      || !lex_force_match (s->lexer, T_RPAREN))
    goto error;

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_eigen (struct eigen_command *eigen)
{
  gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
  if (!A)
    return;
  if (!is_symmetric (A))
    {
      msg (SE, _("Argument of EIGEN must be symmetric."));
      gsl_matrix_free (A);
      return;
    }

  gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
  gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
  gsl_vector v_eval = to_vector (eval);
  gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
  gsl_eigen_symmv (A, &v_eval, evec, w);
  gsl_eigen_symmv_free (w);

  gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);

  gsl_matrix_free (A);

  gsl_matrix_free (eigen->eval->value);
  eigen->eval->value = eval;

  gsl_matrix_free (eigen->evec->value);
  eigen->evec->value = evec;
}
\f
static struct matrix_cmd *
matrix_parse_setdiag (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_SETDIAG,
    .setdiag = { .dst = NULL }
  };

  struct setdiag_command *setdiag = &cmd->setdiag;
  if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
    goto error;

  setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
  if (!setdiag->dst)
    {
      lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
      goto error;
    }
  lex_get (s->lexer);

  if (!lex_force_match (s->lexer, T_COMMA))
    goto error;

  setdiag->expr = matrix_parse_expr (s);
  if (!setdiag->expr)
    goto error;

  if (!lex_force_match (s->lexer, T_RPAREN))
    goto error;

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
{
  gsl_matrix *dst = setdiag->dst->value;
  if (!dst)
    {
      msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
           setdiag->dst->name);
      return;
    }

  gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
  if (!src)
    return;

  size_t n = MIN (dst->size1, dst->size2);
  if (is_scalar (src))
    {
      double d = to_scalar (src);
      for (size_t i = 0; i < n; i++)
        gsl_matrix_set (dst, i, i, d);
    }
  else if (is_vector (src))
    {
      gsl_vector v = to_vector (src);
      for (size_t i = 0; i < n && i < v.size; i++)
        gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
    }
  else
    msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
               "not a %zu×%zu matrix."),
         src->size1, src->size2);
  gsl_matrix_free (src);
}
\f
static struct matrix_cmd *
matrix_parse_svd (struct matrix_state *s)
{
  struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
  *cmd = (struct matrix_cmd) {
    .type = MCMD_SVD,
    .svd = { .expr = NULL }
  };

  struct svd_command *svd = &cmd->svd;
  if (!lex_force_match (s->lexer, T_LPAREN))
    goto error;
  svd->expr = matrix_parse_expr (s);
  if (!svd->expr
      || !lex_force_match (s->lexer, T_COMMA)
      || !matrix_parse_dst_var (s, &svd->u)
      || !lex_force_match (s->lexer, T_COMMA)
      || !matrix_parse_dst_var (s, &svd->s)
      || !lex_force_match (s->lexer, T_COMMA)
      || !matrix_parse_dst_var (s, &svd->v)
      || !lex_force_match (s->lexer, T_RPAREN))
    goto error;

  return cmd;

error:
  matrix_cmd_destroy (cmd);
  return NULL;
}

static void
matrix_cmd_execute_svd (struct svd_command *svd)
{
  gsl_matrix *m = matrix_expr_evaluate (svd->expr);
  if (!m)
    return;

  if (m->size1 >= m->size2)
    {
      gsl_matrix *A = m;
      gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
      gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
      gsl_vector Sv = gsl_matrix_diagonal (S).vector;
      gsl_vector *work = gsl_vector_alloc (A->size2);
      gsl_linalg_SV_decomp (A, V, &Sv, work);
      gsl_vector_free (work);

      matrix_var_set (svd->u, A);
      matrix_var_set (svd->s, S);
      matrix_var_set (svd->v, V);
    }
  else
    {
      gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
      gsl_matrix_transpose_memcpy (At, m);
      gsl_matrix_free (m);

      gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
      gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
      gsl_vector Stv = gsl_matrix_diagonal (St).vector;
      gsl_vector *work = gsl_vector_alloc (At->size2);
      gsl_linalg_SV_decomp (At, Vt, &Stv, work);
      gsl_vector_free (work);

      matrix_var_set (svd->v, At);
      matrix_var_set (svd->s, St);
      matrix_var_set (svd->u, Vt);
    }
}
\f
static bool
matrix_cmd_execute (struct matrix_cmd *cmd)
{
  switch (cmd->type)
    {
    case MCMD_COMPUTE:
      matrix_cmd_execute_compute (&cmd->compute);
      break;

    case MCMD_PRINT:
      matrix_cmd_execute_print (&cmd->print);
      break;

    case MCMD_DO_IF:
      return matrix_cmd_execute_do_if (&cmd->do_if);

    case MCMD_LOOP:
      matrix_cmd_execute_loop (&cmd->loop);
      break;

    case MCMD_BREAK:
      return false;

    case MCMD_DISPLAY:
      matrix_cmd_execute_display (&cmd->display);
      break;

    case MCMD_RELEASE:
      matrix_cmd_execute_release (&cmd->release);
      break;

    case MCMD_SAVE:
      matrix_cmd_execute_save (&cmd->save);
      break;

    case MCMD_READ:
      matrix_cmd_execute_read (&cmd->read);
      break;

    case MCMD_WRITE:
      matrix_cmd_execute_write (&cmd->write);
      break;

    case MCMD_GET:
      matrix_cmd_execute_get (&cmd->get);
      break;

    case MCMD_MSAVE:
      matrix_cmd_execute_msave (&cmd->msave);
      break;

    case MCMD_MGET:
      matrix_cmd_execute_mget (&cmd->mget);
      break;

    case MCMD_EIGEN:
      matrix_cmd_execute_eigen (&cmd->eigen);
      break;

    case MCMD_SETDIAG:
      matrix_cmd_execute_setdiag (&cmd->setdiag);
      break;

    case MCMD_SVD:
      matrix_cmd_execute_svd (&cmd->svd);
      break;
    }

  return true;
}

static void
matrix_cmds_uninit (struct matrix_cmds *cmds)
{
  for (size_t i = 0; i < cmds->n; i++)
    matrix_cmd_destroy (cmds->commands[i]);
  free (cmds->commands);
}

static void
matrix_cmd_destroy (struct matrix_cmd *cmd)
{
  if (!cmd)
    return;

  switch (cmd->type)
    {
    case MCMD_COMPUTE:
      matrix_lvalue_destroy (cmd->compute.lvalue);
      matrix_expr_destroy (cmd->compute.rvalue);
      break;

    case MCMD_PRINT:
      matrix_expr_destroy (cmd->print.expression);
      free (cmd->print.title);
      print_labels_destroy (cmd->print.rlabels);
      print_labels_destroy (cmd->print.clabels);
      break;

    case MCMD_DO_IF:
      for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
        {
          matrix_expr_destroy (cmd->do_if.clauses[i].condition);
          matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
        }
      free (cmd->do_if.clauses);
      break;

    case MCMD_LOOP:
      matrix_expr_destroy (cmd->loop.start);
      matrix_expr_destroy (cmd->loop.end);
      matrix_expr_destroy (cmd->loop.increment);
      matrix_expr_destroy (cmd->loop.top_condition);
      matrix_expr_destroy (cmd->loop.bottom_condition);
      matrix_cmds_uninit (&cmd->loop.commands);
      break;

    case MCMD_BREAK:
      break;

    case MCMD_DISPLAY:
      break;

    case MCMD_RELEASE:
      free (cmd->release.vars);
      break;

    case MCMD_SAVE:
      matrix_expr_destroy (cmd->save.expression);
      fh_unref (cmd->save.outfile);
      string_array_destroy (cmd->save.variables);
      matrix_expr_destroy (cmd->save.names);
      stringi_set_destroy (&cmd->save.strings);
      break;

    case MCMD_READ:
      matrix_lvalue_destroy (cmd->read.dst);
      matrix_expr_destroy (cmd->read.size);
      break;

    case MCMD_WRITE:
      matrix_expr_destroy (cmd->write.expression);
      break;

    case MCMD_GET:
      matrix_lvalue_destroy (cmd->get.dst);
      fh_unref (cmd->get.file);
      free (cmd->get.encoding);
      string_array_destroy (&cmd->get.variables);
      break;

    case MCMD_MSAVE:
      free (cmd->msave.varname_);
      matrix_expr_destroy (cmd->msave.expr);
      matrix_expr_destroy (cmd->msave.factors);
      matrix_expr_destroy (cmd->msave.splits);
      break;

    case MCMD_MGET:
      fh_unref (cmd->mget.file);
      stringi_set_destroy (&cmd->mget.rowtypes);
      break;

    case MCMD_EIGEN:
      matrix_expr_destroy (cmd->eigen.expr);
      break;

    case MCMD_SETDIAG:
      matrix_expr_destroy (cmd->setdiag.expr);
      break;

    case MCMD_SVD:
      matrix_expr_destroy (cmd->svd.expr);
      break;
    }
  free (cmd);
}

struct matrix_command_name
  {
    const char *name;
    struct matrix_cmd *(*parse) (struct matrix_state *);
  };

static const struct matrix_command_name *
matrix_parse_command_name (struct lexer *lexer)
{
  static const struct matrix_command_name commands[] = {
    { "COMPUTE", matrix_parse_compute },
    { "CALL EIGEN", matrix_parse_eigen },
    { "CALL SETDIAG", matrix_parse_setdiag },
    { "CALL SVD", matrix_parse_svd },
    { "PRINT", matrix_parse_print },
    { "DO IF", matrix_parse_do_if },
    { "LOOP", matrix_parse_loop },
    { "BREAK", matrix_parse_break },
    { "READ", matrix_parse_read },
    { "WRITE", matrix_parse_write },
    { "GET", matrix_parse_get },
    { "SAVE", matrix_parse_save },
    { "MGET", matrix_parse_mget },
    { "MSAVE", matrix_parse_msave },
    { "DISPLAY", matrix_parse_display },
    { "RELEASE", matrix_parse_release },
  };
  static size_t n = sizeof commands / sizeof *commands;

  for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
    if (lex_match_phrase (lexer, c->name))
      return c;
  return NULL;
}

static struct matrix_cmd *
matrix_parse_command (struct matrix_state *s)
{
  size_t nesting_level = SIZE_MAX;

  struct matrix_cmd *c = NULL;
  const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
  if (!cmd)
    lex_error (s->lexer, _("Unknown matrix command."));
  else if (!cmd->parse)
    lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
               cmd->name);
  else
    {
      nesting_level = output_open_group (
        group_item_create_nocopy (utf8_to_title (cmd->name),
                                  utf8_to_title (cmd->name)));
      c = cmd->parse (s);
    }

  if (c)
    lex_end_of_command (s->lexer);
  lex_discard_rest_of_command (s->lexer);
  if (nesting_level != SIZE_MAX)
    output_close_groups (nesting_level);

  return c;
}

int
cmd_matrix (struct lexer *lexer, struct dataset *ds)
{
  if (!lex_force_match (lexer, T_ENDCMD))
    return CMD_FAILURE;

  struct matrix_state state = {
    .session = dataset_session (ds),
    .lexer = lexer,
    .vars = HMAP_INITIALIZER (state.vars),
  };

  for (;;)
    {
      while (lex_match (lexer, T_ENDCMD))
        continue;
      if (lex_token (lexer) == T_STOP)
        {
          msg (SE, _("Unexpected end of input expecting matrix command."));
          break;
        }

      if (lex_match_phrase (lexer, "END MATRIX"))
        break;

      struct matrix_cmd *c = matrix_parse_command (&state);
      if (c)
        {
          matrix_cmd_execute (c);
          matrix_cmd_destroy (c);
        }
    }

  struct matrix_var *var, *next;
  HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
    {
      free (var->name);
      gsl_matrix_free (var->value);
      hmap_delete (&state.vars, &var->hmap_node);
      free (var);
    }
  hmap_destroy (&state.vars);
  fh_unref (state.prev_save_outfile);
  if (state.common)
    {
      dict_unref (state.common->dict);
      casewriter_destroy (state.common->writer);
      free (state.common);
    }
  fh_unref (state.prev_read_file);
  for (size_t i = 0; i < state.n_read_files; i++)
    read_file_destroy (state.read_files[i]);
  free (state.read_files);
  fh_unref (state.prev_write_file);
  for (size_t i = 0; i < state.n_write_files; i++)
    write_file_destroy (state.write_files[i]);
  free (state.write_files);

  return CMD_SUCCESS;
}