X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fmatrix.c;h=b709d65e01e0066639d09cc227494c026d31fbd1;hb=5e6a3c6f8b078875bef5374874f7f004b3077509;hp=c0c80eda77f3a8005def95de95ef9b24d78c6a58;hpb=e80c032386b6524fe760efc6f10bd9adb8a6ac20;p=pspp

diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c
index c0c80eda77..b709d65e01 100644
--- a/src/language/stats/matrix.c
+++ b/src/language/stats/matrix.c
@@ -17,6 +17,7 @@
 #include <config.h>
 
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_cdf.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
@@ -53,12 +54,15 @@
 #include "libpspp/string-array.h"
 #include "libpspp/stringi-set.h"
 #include "libpspp/u8-line.h"
+#include "math/distributions.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/c-strcase.h"
+#include "gl/ftoastr.h"
 #include "gl/minmax.h"
 #include "gl/xsize.h"
 
@@ -69,7 +73,7 @@
 struct matrix_var
   {
     struct hmap_node hmap_node;
-    const char *name;
+    char *name;
     gsl_matrix *value;
   };
 
@@ -80,10 +84,15 @@ struct msave_common
     struct string_array variables;
     struct string_array fnames;
     struct string_array snames;
-    bool has_factors;
-    bool has_splits;
     size_t n_varnames;
 
+    /* Collects and owns factors and splits.  The individual msave_command
+       structs point to these but do not own them. */
+    struct matrix_expr **factors;
+    size_t n_factors, allocated_factors;
+    struct matrix_expr **splits;
+    size_t n_splits, allocated_splits;
+
     /* Execution state. */
     struct dictionary *dict;
     struct casewriter *writer;
@@ -96,6 +105,30 @@ struct read_file
     char *encoding;
   };
 
+struct write_file
+  {
+    struct file_handle *file;
+    struct dfm_writer *writer;
+    char *encoding;
+    struct u8_line *held;
+  };
+
+struct save_file
+  {
+    struct file_handle *file;
+    struct dataset *dataset;
+
+    /* Parameters from parsing the first SAVE command for 'file'. */
+    struct string_array variables;
+    struct matrix_expr *names;
+    struct stringi_set strings;
+
+    /* Results from the first (only) attempt to open this save_file. */
+    bool error;
+    struct casewriter *writer;
+    struct dictionary *dict;
+  };
+
 struct matrix_state
   {
     struct dataset *dataset;
@@ -103,13 +136,19 @@ struct matrix_state
     struct lexer *lexer;
     struct hmap vars;
     bool in_loop;
-    struct file_handle *prev_save_outfile;
-    struct file_handle *prev_write_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;
+
+    struct file_handle *prev_save_file;
+    struct save_file **save_files;
+    size_t n_save_files;
   };
 
 static struct matrix_var *
@@ -141,87 +180,192 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value)
   var->value = value;
 }
 
-#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)
+#define MATRIX_FUNCTIONS                                \
+    F(ABS,      "ABS",      m_m_e, NULL)                \
+    F(ALL,      "ALL",      d_m, NULL)                  \
+    F(ANY,      "ANY",      d_m, NULL)                  \
+    F(ARSIN,    "ARSIN",    m_m_e, "a[-1,1]")         \
+    F(ARTAN,    "ARTAN",    m_m_e, NULL)                \
+    F(BLOCK,    "BLOCK",    m_any, NULL)                \
+    F(CHOL,     "CHOL",     m_m, NULL)                  \
+    F(CMIN,     "CMIN",     m_m, NULL)                  \
+    F(CMAX,     "CMAX",     m_m, NULL)                  \
+    F(COS,      "COS",      m_m_e, NULL)                \
+    F(CSSQ,     "CSSQ",     m_m, NULL)                  \
+    F(CSUM,     "CSUM",     m_m, NULL)                  \
+    F(DESIGN,   "DESIGN",   m_m, NULL)                  \
+    F(DET,      "DET",      d_m, NULL)                  \
+    F(DIAG,     "DIAG",     m_m, NULL)                  \
+    F(EVAL,     "EVAL",     m_m, NULL)                  \
+    F(EXP,      "EXP",      m_m_e, NULL)                \
+    F(GINV,     "GINV",     m_m, NULL)                  \
+    F(GRADE,    "GRADE",    m_m, NULL)                  \
+    F(GSCH,     "GSCH",     m_m, NULL)                  \
+    F(IDENT,    "IDENT",    IDENT, NULL)                \
+    F(INV,      "INV",      m_m, NULL)                  \
+    F(KRONEKER, "KRONEKER", m_mm, NULL)                 \
+    F(LG10,     "LG10",     m_m_e, "a>0")               \
+    F(LN,       "LN",       m_m_e, "a>0")               \
+    F(MAGIC,    "MAGIC",    m_d, NULL)                  \
+    F(MAKE,     "MAKE",     m_ddd, NULL)                \
+    F(MDIAG,    "MDIAG",    m_v, NULL)                  \
+    F(MMAX,     "MMAX",     d_m, NULL)                  \
+    F(MMIN,     "MMIN",     d_m, NULL)                  \
+    F(MOD,      "MOD",      m_md, NULL)                 \
+    F(MSSQ,     "MSSQ",     d_m, NULL)                  \
+    F(MSUM,     "MSUM",     d_m, NULL)                  \
+    F(NCOL,     "NCOL",     d_m, NULL)                  \
+    F(NROW,     "NROW",     d_m, NULL)                  \
+    F(RANK,     "RANK",     d_m, NULL)                  \
+    F(RESHAPE,  "RESHAPE",  m_mdd, NULL)                \
+    F(RMAX,     "RMAX",     m_m, NULL)                  \
+    F(RMIN,     "RMIN",     m_m, NULL)                  \
+    F(RND,      "RND",      m_m_e, NULL)                \
+    F(RNKORDER, "RNKORDER", m_m, NULL)                  \
+    F(RSSQ,     "RSSQ",     m_m, NULL)                  \
+    F(RSUM,     "RSUM",     m_m, NULL)                  \
+    F(SIN,      "SIN",      m_m_e, NULL)                \
+    F(SOLVE,    "SOLVE",    m_mm, NULL)                 \
+    F(SQRT,     "SQRT",     m_m, NULL)                  \
+    F(SSCP,     "SSCP",     m_m, NULL)                  \
+    F(SVAL,     "SVAL",     m_m, NULL)                  \
+    F(SWEEP,    "SWEEP",    m_md, NULL)                 \
+    F(T,        "T",        m_m, NULL)                  \
+    F(TRACE,    "TRACE",    d_m, NULL)                  \
+    F(TRANSPOS, "TRANSPOS", m_m, NULL)                  \
+    F(TRUNC,    "TRUNC",    m_m_e, NULL)                \
+    F(UNIFORM,  "UNIFORM",  m_dd, NULL)                 \
+    F(PDF_BETA, "PDF.BETA", m_mdd_e, "a[0,1] b>0 c>0")  \
+    F(CDF_BETA, "CDF.BETA", m_mdd_e, "a[0,1] b>0 c>0")  \
+    F(IDF_BETA, "IDF.BETA", m_mdd_e, "a[0,1] b>0 c>0")  \
+    F(RV_BETA,  "RV.BETA",  d_dd, "a>0 b>0")            \
+    F(NCDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0")       \
+    F(NPDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0") \
+    /* XXX CDF.BVNOR */ \
+    F(PDF_BVNOR, "PDF.BVNOR", m_mdd_e, "c[-1,1]") \
+    F(CDF_CAUCHY, "CDF.CAUCHY", m_mdd_e, "c>0") \
+    F(IDF_CAUCHY, "IDF.CAUCHY", m_mdd_e, "a(0,1) c>0") \
+    F(PDF_CAUCHY, "PDF.CAUCHY", m_mdd_e, "c>0") \
+    F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
+    F(CDF_CHISQ, "CDF.CHISQ", m_md_e, "a>=0 b>0") \
+    F(IDF_CHISQ, "IDF.CHISQ", m_md_e, "a[0,1) b>0") \
+    F(PDF_CHISQ, "PDF.CHISQ", m_md_e, "a>=0 b>0") \
+    F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
+    F(SIG_CHISQ, "SIG.CHISQ", m_md_e, "a>=0 b>0") \
+    F(CDF_EXP, "CDF.EXP", m_md_e, "a>=0 b>=0") \
+    F(IDF_EXP, "IDF.EXP", m_md_e, "a[0,1) b>0") \
+    F(PDF_EXP, "PDF.EXP", m_md_e, "a>=0 b>0") \
+    F(RV_EXP, "RV.EXP", d_d, "a>0") \
+    F(PDF_XPOWER, "PDF.XPOWER", m_mdd_e, "b>0 c>=0") \
+    F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
+    F(CDF_F, "CDF.F", m_mdd_e, "a>=0 b>0 c>0") \
+    F(IDF_F, "IDF.F", m_mdd_e, "a[0,1) b>0 c>0") \
+    F(PDF_F, "PDF.F", m_mdd_e, "a>=0 b>0 c>0") \
+    F(RV_F, "RV.F", d_dd, "a>0 b>0") \
+    F(SIG_F, "SIG.F", m_mdd_e, "a>=0 b>0 c>0")
+
+struct matrix_function_properties
+  {
+    const char *name;
+    const char *constraints;
+  };
 
+/*
+()
+(a)
+(a > 0)
+(a >= 0)
+(a > 0 && a < 1)
+(a > 0 && a <= 1)
+(a >= 0 && a <= 1)
+(a > 0 && a < 1, b > 0)
+(a >= 0 && a < 1, b > 0)
+(a >= 0 && a <= 1, b > 0)
+(a == 0 || a == 1, b >= 0 && b <= 1)
+(a >= 0 && a < 1, b > 0, c > 0)
+(a >= 0 && a <= 1, b > 0, c > 0)
+(a >= 0 && a < 1, b >= 1, c >= 1)
+(a >= 0 && a <= 1, b, c)
+(a > 0 && a < 1, b, c > 0)
+(a >= 0 && a <= 1, b <= c, c)
+(a > 0 && a == floor (a),
+(a >= 0 && a == floor (a) && a <= b,
+(a >= 0 && a == floor (a) && a <= d,
+(a >= 0 && a == floor (a), b > 0)
+(a > 0 && a == floor (a), b >= 0 && b <= 1)
+(a > 0, b > 0)
+(a > 0, b >= 0)
+(a >= 0, b > 0)
+(a >= 0, b > 0, c)
+(a > 0, b > 0, c > 0)
+(a >= 0, b > 0, c > 0)
+(a >= 0, b > 0, c > 0, d > 0)
+(a >= 0, b > 0, c > 0, d >= 0)
+(a > 0, b >= 1, c >= 1)
+(a >= 1 && a == floor (a), b >= 0 && b <= 1)
+(a >= 1, b > 0 && b <= 1)
+(a >= 1, b == floor (b), c > 0 && c <= 1)
+(a, b)
+(a, b > 0)
+(a, b > 0 && b <= 2)
+(a, b > 0 && b <= 2, c >= -1 && c <= 1)
+(a, b > 0 && b == floor (c), c >= 0 && c <= 1)
+(a, b > 0, c)
+(a, b > 0, c > 0)
+(a, b > 0, c >= 0)
+(a <= b, b)
+(a >= b, b > 0, c > 0)
+(a, b, c)
+(a, b, c > 0)
+(a, b, c >= -1 && c <= 1)
+(a <= c, b <= a, c)
+(a == floor (a), b > 0 && b <= 1)
+b > 0 && b == floor (b),
+b > 0 && b == floor (b) && b <= a,
+c >= 0 && c <= 1)
+c > 0 && c == floor (c) && c <= a)
+c > 0 && c == floor (c) && c <= b,
+d > 0 && d == floor (d) && d <= b)
+*/
+
+enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
+enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
 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_m_e_MIN_ARGS = 1, m_m_e_MAX_ARGS = 1 };
 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
+enum { m_md_e_MIN_ARGS = 2, m_md_e_MAX_ARGS = 2 };
 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
+enum { m_mdd_e_MIN_ARGS = 3, m_mdd_e_MAX_ARGS = 3 };
+enum { m_mddd_e_MIN_ARGS = 4, m_mddd_e_MAX_ARGS = 4 };
 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 double matrix_proto_d_d (double);
+typedef double matrix_proto_d_dd (double, double);
 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 double matrix_proto_m_m_e (double);
 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
+typedef double matrix_proto_m_md_e (double, double);
 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
+typedef double matrix_proto_m_mdd_e (double, double, double);
+typedef double matrix_proto_m_mddd_e (double, double, 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;
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+    static matrix_proto_##PROTO matrix_eval_##ENUM;
 MATRIX_FUNCTIONS
 #undef F
 
@@ -232,7 +376,7 @@ struct matrix_expr
     enum matrix_op
       {
         /* Functions. */
-#define F(NAME, PROTOTYPE) MOP_F_##NAME,
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
         MATRIX_FUNCTIONS
 #undef F
 
@@ -307,7 +451,7 @@ matrix_expr_destroy (struct matrix_expr *e)
 
   switch (e->op)
     {
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
 MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -340,6 +484,7 @@ MATRIX_FUNCTIONS
     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:
@@ -451,7 +596,7 @@ matrix_parse_curly_semi (struct matrix_state *s)
       for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
 
 static bool
-is_vector (gsl_matrix *m)
+is_vector (const gsl_matrix *m)
 {
   return m->size1 <= 1 || m->size2 <= 1;
 }
@@ -465,12 +610,10 @@ to_vector (gsl_matrix *m)
 }
 
 
-static gsl_matrix *
-matrix_eval_ABS (gsl_matrix *m)
+static double
+matrix_eval_ABS (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = fabs (*d);
-  return m;
+  return fabs (d);
 }
 
 static double
@@ -491,20 +634,16 @@ matrix_eval_ANY (gsl_matrix *m)
   return 0.0;
 }
 
-static gsl_matrix *
-matrix_eval_ARSIN (gsl_matrix *m)
+static double
+matrix_eval_ARSIN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = asin (*d);
-  return m;
+  return asin (d);
 }
 
-static gsl_matrix *
-matrix_eval_ARTAN (gsl_matrix *m)
+static double
+matrix_eval_ARTAN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = atan (*d);
-  return m;
+  return atan (d);
 }
 
 static gsl_matrix *
@@ -586,12 +725,10 @@ matrix_eval_CMIN (gsl_matrix *m)
   return matrix_eval_col_extremum (m, true);
 }
 
-static gsl_matrix *
-matrix_eval_COS (gsl_matrix *m)
+static double
+matrix_eval_COS (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = cos (*d);
-  return m;
+  return cos (d);
 }
 
 static gsl_matrix *
@@ -751,12 +888,10 @@ matrix_eval_EVAL (gsl_matrix *m)
   return eval;
 }
 
-static gsl_matrix *
-matrix_eval_EXP (gsl_matrix *m)
+static double
+matrix_eval_EXP (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = exp (*d);
-  return m;
+  return exp (d);
 }
 
 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
@@ -809,22 +944,23 @@ matrix_eval_GINV (gsl_matrix *A)
       gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
 
   /* two dot products to obtain pseudoinverse */
-  tmp_mat = gsl_matrix_alloc (m, n);
-  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat);
+  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_mat, 0., A_pinv);
+      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_mat, U, 0., A_pinv);
+      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);
@@ -902,7 +1038,7 @@ 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)."),
+                 "as rows, but it has dimensions %zu×%zu."),
            v->size1, v->size2);
       return NULL;
     }
@@ -937,9 +1073,10 @@ matrix_eval_GSCH (gsl_matrix *v)
 
   if (ux < v->size1)
     {
-      msg (SE, _("Argument to GSCH with dimensions (%zu,%zu) contains only "
+      msg (SE, _("%zu×%zu argument to GSCH contains only "
                  "%zu linearly independent columns."),
            v->size1, v->size2, ux);
+      gsl_matrix_free (u);
       return NULL;
     }
 
@@ -1000,20 +1137,16 @@ matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
   return k;
 }
 
-static gsl_matrix *
-matrix_eval_LG10 (gsl_matrix *m)
+static double
+matrix_eval_LG10 (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = log10 (*d);
-  return m;
+  return log10 (d);
 }
 
-static gsl_matrix *
-matrix_eval_LN (gsl_matrix *m)
+static double
+matrix_eval_LN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = log (*d);
-  return m;
+  return log (d);
 }
 
 static void
@@ -1200,7 +1333,7 @@ matrix_eval_MAKE (double r, double c, double value)
 static gsl_matrix *
 matrix_eval_MDIAG (gsl_vector *v)
 {
-  gsl_matrix *m = gsl_matrix_alloc (v->size, v->size);
+  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;
@@ -1339,12 +1472,10 @@ matrix_eval_RMIN (gsl_matrix *m)
   return matrix_eval_row_extremum (m, true);
 }
 
-static gsl_matrix *
-matrix_eval_RND (gsl_matrix *m)
+static double
+matrix_eval_RND (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = rint (*d);
-  return m;
+  return rint (d);
 }
 
 struct rank
@@ -1425,12 +1556,10 @@ matrix_eval_RSUM (gsl_matrix *m)
   return matrix_eval_row_sum (m, false);
 }
 
-static gsl_matrix *
-matrix_eval_SIN (gsl_matrix *m)
+static double
+matrix_eval_SIN (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = sin (*d);
-  return m;
+  return sin (d);
 }
 
 static gsl_matrix *
@@ -1439,8 +1568,8 @@ 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)."),
+                 "rows, but the first argument has dimensions %zu×%zu and "
+                 "the second %zu×%zu."),
            m1->size1, m1->size2,
            m2->size1, m2->size2);
       return NULL;
@@ -1588,12 +1717,10 @@ matrix_eval_TRANSPOS (gsl_matrix *m)
     }
 }
 
-static gsl_matrix *
-matrix_eval_TRUNC (gsl_matrix *m)
+static double
+matrix_eval_TRUNC (double d)
 {
-  MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
-    *d = trunc (*d);
-  return m;
+  return trunc (d);
 }
 
 static gsl_matrix *
@@ -1618,6 +1745,168 @@ matrix_eval_UNIFORM (double r_, double c_)
   return m;
 }
 
+static double
+matrix_eval_PDF_BETA (double x, double a, double b)
+{
+  return gsl_ran_beta_pdf (x, a, b);
+}
+
+static double
+matrix_eval_CDF_BETA (double x, double a, double b)
+{
+  return gsl_cdf_beta_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_BETA (double P, double a, double b)
+{
+  return gsl_cdf_beta_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_RV_BETA (double a, double b)
+{
+  return gsl_ran_beta (get_rng (), a, b);
+}
+
+static double
+matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
+{
+  return ncdf_beta (x, a, b, lambda);
+}
+
+static double
+matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
+{
+  return npdf_beta (x, a, b, lambda);
+}
+
+static double
+matrix_eval_PDF_BVNOR (double x0, double x1, double r)
+{
+  return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
+}
+
+static double
+matrix_eval_CDF_CAUCHY (double x, double a, double b)
+{
+  return gsl_cdf_cauchy_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_CAUCHY (double P, double a, double b)
+{
+  return a + b * gsl_cdf_cauchy_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_CAUCHY (double x, double a, double b)
+{
+  return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
+}
+
+static double
+matrix_eval_RV_CAUCHY (double a, double b)
+{
+  return a + b * gsl_ran_cauchy (get_rng (), 1);
+}
+
+static double
+matrix_eval_CDF_CHISQ (double x, double df)
+{
+  return gsl_cdf_chisq_P (x, df);
+}
+
+static double
+matrix_eval_IDF_CHISQ (double P, double df)
+{
+  return gsl_cdf_chisq_Pinv (P, df);
+}
+
+static double
+matrix_eval_PDF_CHISQ (double x, double df)
+{
+  return gsl_ran_chisq_pdf (x, df);
+}
+
+static double
+matrix_eval_RV_CHISQ (double df)
+{
+  return gsl_ran_chisq (get_rng (), df);
+}
+
+static double
+matrix_eval_SIG_CHISQ (double x, double df)
+{
+  return gsl_cdf_chisq_Q (x, df);
+}
+
+static double
+matrix_eval_CDF_EXP (double x, double a)
+{
+  return gsl_cdf_exponential_P (x, 1. / a);
+}
+
+static double
+matrix_eval_IDF_EXP (double P, double a)
+{
+  return gsl_cdf_exponential_Pinv (P, 1. / a);
+}
+
+static double
+matrix_eval_PDF_EXP (double x, double a)
+{
+  return gsl_ran_exponential_pdf (x, 1. / a);
+}
+
+static double
+matrix_eval_RV_EXP (double a)
+{
+  return gsl_ran_exponential (get_rng (), 1. / a);
+}
+
+static double
+matrix_eval_PDF_XPOWER (double x, double a, double b)
+{
+  return gsl_ran_exppow_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_XPOWER (double a, double b)
+{
+  return gsl_ran_exppow (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_F (double x, double df1, double df2)
+{
+  return gsl_cdf_fdist_P (x, df1, df2);
+}
+
+static double
+matrix_eval_IDF_F (double P, double df1, double df2)
+{
+  return idf_fdist (P, df1, df2);
+}
+
+static double
+matrix_eval_RV_F (double df1, double df2)
+{
+  return gsl_ran_fdist (get_rng (), df1, df2);
+}
+
+static double
+matrix_eval_PDF_F (double x, double df1, double df2)
+{
+  return gsl_ran_fdist_pdf (x, df1, df2);
+}
+
+static double
+matrix_eval_SIG_F (double x, double df1, double df2)
+{
+  return gsl_cdf_fdist_Q (x, df1, df2);
+}
+
 struct matrix_function
   {
     const char *name;
@@ -1627,12 +1916,58 @@ struct matrix_function
 
 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
 
+static bool
+word_matches (const char **test, const char **name)
+{
+  size_t test_len = strcspn (*test, ".");
+  size_t name_len = strcspn (*name, ".");
+  if (test_len == name_len)
+    {
+      if (buf_compare_case (*test, *name, test_len))
+        return false;
+    }
+  else if (test_len < 3 || test_len > name_len)
+    return false;
+  else
+    {
+      if (buf_compare_case (*test, *name, test_len))
+        return false;
+    }
+
+  *test += test_len;
+  *name += name_len;
+  if (**test != **name)
+    return false;
+
+  if (**test == '.')
+    {
+      (*test)++;
+      (*name)++;
+    }
+  return true;
+}
+
+/* Returns 0 if TOKEN and FUNC do not match,
+   1 if TOKEN is an acceptable abbreviation for FUNC,
+   2 if TOKEN equals FUNC. */
+static int
+compare_function_names (const char *token_, const char *func_)
+{
+  const char *token = token_;
+  const char *func = func_;
+  while (*token || *func)
+    if (!word_matches (&token, &func))
+      return 0;
+  return !c_strcasecmp (token_, func_) ? 2 : 1;
+}
+
 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 },
+#define F(ENUM, STRING, PROTO, CONSTRAINTS)                             \
+      { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
       MATRIX_FUNCTIONS
 #undef F
     };
@@ -1640,7 +1975,7 @@ matrix_parse_function_name (const char *token)
 
   for (size_t i = 0; i < N_FUNCTIONS; i++)
     {
-      if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
+      if (compare_function_names (token, functions[i].name) > 0)
         return &functions[i];
     }
   return NULL;
@@ -1676,6 +2011,68 @@ read_file_open (struct read_file *rf)
   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)
@@ -2115,7 +2512,7 @@ matrix_op_eval (enum matrix_op op, double a, double b)
     case MOP_OR: return (a > 0) || (b > 0);
     case MOP_XOR: return (a > 0) != (b > 0);
 
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
       MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -2160,7 +2557,7 @@ matrix_op_name (enum matrix_op op)
     case MOP_OR: return "OR";
     case MOP_XOR: return "XOR";
 
-#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
       MATRIX_FUNCTIONS
 #undef F
     case MOP_NEGATE:
@@ -2238,8 +2635,7 @@ matrix_expr_evaluate_elementwise (enum matrix_op op,
   else
     {
       msg (SE, _("Operands to %s must have the same dimensions or one "
-                 "must be a scalar, not matrices with dimensions (%zu,%zu) "
-                 "and (%zu,%zu)."),
+                 "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;
     }
@@ -2253,7 +2649,7 @@ matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
 
   if (a->size2 != b->size1)
     {
-      msg (SE, _("Matrices with dimensions (%zu,%zu) and (%zu,%zu) are "
+      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;
@@ -2293,14 +2689,14 @@ matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
   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)."),
+                 "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)."),
+                 "right-hand side, not a matrix with dimensions %zu×%zu."),
            b->size1, b->size2);
       return NULL;
     }
@@ -2313,27 +2709,37 @@ matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
     }
   long int bl = bf;
 
-  gsl_matrix *tmp = gsl_matrix_alloc (x->size1, x->size2);
-  gsl_matrix *y = gsl_matrix_alloc (x->size1, x->size2);
+  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, &tmp);
-        square_matrix (&x, &tmp);
+        mul_matrix (&y, x, y, &t);
+        square_matrix (&x, &t);
       }
     else
-      square_matrix (&x, &tmp);
+      square_matrix (&x, &t);
 
-  mul_matrix (&y, x, y, &tmp);
+  mul_matrix (&y, x, y, &t);
   if (bf < 0)
     invert_matrix (y);
 
-  if (tmp != x_)
-    gsl_matrix_free (tmp);
+  /* 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;
 }
 
@@ -2357,8 +2763,8 @@ matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
       return NULL;
     }
 
-  long int n = (end > start && by > 0 ? (end - start + by) / by
-                : end < start && by < 0 ? (start - end - by) / -by
+  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++)
@@ -2441,6 +2847,7 @@ matrix_to_vector (gsl_matrix *m)
   assert (!v.owner);
   v.owner = 1;
   m->owner = 0;
+  gsl_matrix_free (m);
   return xmemdup (&v, sizeof v);
 }
 
@@ -2458,7 +2865,7 @@ struct index_vector
 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
 
 static bool
-matrix_normalize_index_vector (gsl_matrix *m, size_t size,
+matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
                                enum index_type index_type, size_t other_size,
                                struct index_vector *iv)
 {
@@ -2470,26 +2877,26 @@ matrix_normalize_index_vector (gsl_matrix *m, size_t size,
             {
             case IV_VECTOR:
               msg (SE, _("Vector index must be scalar or vector, not a "
-                         "matrix with dimensions (%zu,%zu)."),
+                         "%zu×%zu matrix."),
                    m->size1, m->size2);
               break;
 
             case IV_ROW:
               msg (SE, _("Matrix row index must be scalar or vector, not a "
-                         "matrix with dimensions (%zu,%zu)."),
+                         "%zu×%zu matrix."),
                    m->size1, m->size2);
               break;
 
             case IV_COLUMN:
               msg (SE, _("Matrix column index must be scalar or vector, not a "
-                         "matrix with dimensions (%zu,%zu)."),
+                         "%zu×%zu matrix."),
                    m->size1, m->size2);
               break;
             }
           return false;
         }
 
-      gsl_vector v = to_vector (m);
+      gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
       *iv = (struct index_vector) {
         .indexes = xnmalloc (v.size, sizeof *iv->indexes),
         .n = v.size,
@@ -2507,14 +2914,14 @@ matrix_normalize_index_vector (gsl_matrix *m, size_t size,
                   break;
 
                 case IV_ROW:
-                  msg (SE, _("%g is not a valid row index for a matrix "
-                             "with dimensions (%zu,%zu)."),
+                  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 matrix "
-                             "with dimensions (%zu,%zu)."),
+                  msg (SE, _("%g is not a valid column index for "
+                             "a %zu×%zu matrix."),
                        index, other_size, size);
                   break;
                 }
@@ -2543,8 +2950,8 @@ matrix_expr_evaluate_vec_all (gsl_matrix *sm)
 {
   if (!is_vector (sm))
     {
-      msg (SE, _("Vector index operator must be applied to vector, "
-                 "not a matrix with dimensions (%zu,%zu)."),
+      msg (SE, _("Vector index operator may not be applied to "
+                 "a %zu×%zu matrix."),
            sm->size1, sm->size2);
       return NULL;
     }
@@ -2608,10 +3015,10 @@ matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
   return dm;
 }
 
-#define F(NAME, PROTOTYPE)                              \
-  static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
-    const char *name, gsl_matrix *[], size_t,           \
-    matrix_proto_##PROTOTYPE *);
+#define F(ENUM, STRING, PROTO, CONSTRAINTS)                             \
+  static gsl_matrix *matrix_expr_evaluate_##PROTO (                     \
+    const struct matrix_function_properties *, gsl_matrix *[], size_t,  \
+    matrix_proto_##PROTO *);
 MATRIX_FUNCTIONS
 #undef F
 
@@ -2621,7 +3028,7 @@ 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, "
-                 "but it has dimensions (%zu,%zu)."),
+                 "not a %zu×%zu matrix."),
            name, index + 1, subs[index]->size1, subs[index]->size2);
       return false;
     }
@@ -2634,7 +3041,7 @@ 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, "
-                 "but it has dimensions (%zu,%zu)."),
+                 "not a %zu×%zu matrix."),
            name, index + 1, subs[index]->size1, subs[index]->size2);
       return false;
     }
@@ -2653,41 +3060,194 @@ to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
   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)
+static int
+parse_constraint_value (const char **constraintsp)
 {
-  assert (n_subs == 1);
+  char *tail;
+  long retval = strtol (*constraintsp, &tail, 10);
+  *constraintsp = tail;
+  return retval;
+}
 
-  double d;
-  return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
+static bool
+check_constraints (const struct matrix_function_properties *props,
+                   gsl_matrix *args[], size_t n_args)
+{
+  const char *constraints = props->constraints;
+  if (!constraints)
+    return true;
+
+  size_t arg_index = SIZE_MAX;
+  while (*constraints)
+    {
+      if (*constraints >= 'a' && *constraints <= 'd')
+        {
+          arg_index = *constraints++ - 'a';
+          assert (arg_index < n_args);
+        }
+      else if (*constraints == '[' || *constraints == '(')
+        {
+          assert (arg_index < n_args);
+          bool open_lower = *constraints++ == '(';
+          int minimum = parse_constraint_value (&constraints);
+          assert (*constraints == ',');
+          constraints++;
+          int maximum = parse_constraint_value (&constraints);
+          assert (*constraints == ']' || *constraints == ')');
+          bool open_upper = *constraints++ == ')';
+
+          MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+            if ((open_lower ? *d <= minimum : *d < minimum)
+                || (open_upper ? *d >= maximum : *d > maximum))
+              {
+                if (!is_scalar (args[arg_index]))
+                  msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
+                             "function %s has value %g, which is outside "
+                             "the valid range %c%d,%d%c."),
+                       y + 1, x + 1, arg_index + 1, props->name, *d,
+                       open_lower ? '(' : '[',
+                       minimum, maximum,
+                       open_upper ? ')' : ']');
+                else
+                  msg (ME, _("Argument %zu to matrix function %s has value %g, "
+                             "which is outside the valid range %c%d,%d%c."),
+                       arg_index + 1, props->name, *d,
+                       open_lower ? '(' : '[',
+                       minimum, maximum,
+                       open_upper ? ')' : ']');
+                return false;
+              }
+        }
+      else if (*constraints == '>' || *constraints == '<')
+        {
+          bool greater = *constraints++ == '>';
+          bool equal;
+          if (*constraints == '=')
+            {
+              equal = true;
+              constraints++;
+            }
+          else
+            equal = false;
+          int comparand = parse_constraint_value (&constraints);
+
+          assert (arg_index < n_args);
+          MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+            if (greater
+                ? (equal ? !(*d >= comparand) : !(*d > comparand))
+                : (equal ? !(*d <= comparand) : !(*d < comparand)))
+              {
+                struct string s = DS_EMPTY_INITIALIZER;
+                if (!is_scalar (args[arg_index]))
+                  ds_put_format (&s, _("Row %zu, column %zu of argument %zu "
+                                       "to matrix function %s has "
+                                       "invalid value %g."),
+                                 y + 1, x + 1, arg_index + 1, props->name, *d);
+                else
+                  ds_put_format (&s, _("Argument %zu to matrix function %s "
+                                       "has invalid value %g."),
+                                 arg_index + 1, props->name, *d);
+
+                ds_put_cstr (&s, "  ");
+                if (greater && equal)
+                  ds_put_format (&s, _("This argument must be greater than or "
+                                       "equal to %d."), comparand);
+                else if (greater && !equal)
+                  ds_put_format (&s, _("This argument must be greater than %d."),
+                                 comparand);
+                else if (equal)
+                  ds_put_format (&s, _("This argument must be less than or "
+                                       "equal to %d."), comparand);
+                else
+                  ds_put_format (&s, _("This argument must be less than %d."),
+                                 comparand);
+                msg (ME, ds_cstr (&s));
+                ds_destroy (&s);
+                return false;
+              }
+        }
+      else
+        {
+          assert (*constraints == ' ');
+          constraints++;
+        }
+    }
+  return true;
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_dd (const char *name,
+matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
                            gsl_matrix *subs[], size_t n_subs,
-                           matrix_proto_m_dd *f)
+                           matrix_proto_d_d *f)
+{
+  assert (n_subs == 1);
+
+  double d;
+  if (!to_scalar_args (props->name, subs, n_subs, &d))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f (d));
+  return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
+                           gsl_matrix *subs[], size_t n_subs,
+                           matrix_proto_d_dd *f)
+{
+  assert (n_subs == 2);
+
+  double d[2];
+  if (!to_scalar_args (props->name, subs, n_subs, d))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  gsl_matrix *m = gsl_matrix_alloc (1, 1);
+  gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
+  return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
+                          gsl_matrix *subs[], size_t n_subs,
+                          matrix_proto_m_d *f)
+{
+  assert (n_subs == 1);
+
+  double d;
+  return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
+                           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;
+  return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_ddd (const char *name,
+matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
                            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;
+  return to_scalar_args (props->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,
+matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_m *f)
 {
@@ -2696,29 +3256,102 @@ matrix_expr_evaluate_m_m (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_md (const char *name UNUSED,
+matrix_expr_evaluate_m_m_e (const struct matrix_function_properties *props,
+                            gsl_matrix *subs[], size_t n_subs,
+                            matrix_proto_m_m_e *f)
+{
+  assert (n_subs == 1);
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+      *a = f (*a);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_md *f)
 {
   assert (n_subs == 2);
-  if (!check_scalar_arg (name, subs, 1))
+  if (!check_scalar_arg (props->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,
+matrix_expr_evaluate_m_md_e (const struct matrix_function_properties *props,
+                             gsl_matrix *subs[], size_t n_subs,
+                             matrix_proto_m_md_e *f)
+{
+  assert (n_subs == 2);
+  if (!check_scalar_arg (props->name, subs, 1))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props 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))
+  if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->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,
+matrix_expr_evaluate_m_mdd_e (const struct matrix_function_properties *props,
+                              gsl_matrix *subs[], size_t n_subs,
+                              matrix_proto_m_mdd_e *f)
+{
+  assert (n_subs == 3);
+  if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  double c = to_scalar (subs[2]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b, c);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mddd_e (const struct matrix_function_properties *props,
+                               gsl_matrix *subs[], size_t n_subs,
+                               matrix_proto_m_mddd_e *f)
+{
+  assert (n_subs == 4);
+  for (size_t i = 1; i < 4; i++)
+    if (!check_scalar_arg (props->name, subs, i))
+    return NULL;
+
+  if (!check_constraints (props, subs, n_subs))
+    return NULL;
+
+  double b = to_scalar (subs[1]);
+  double c = to_scalar (subs[2]);
+  double d = to_scalar (subs[3]);
+  MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+    *a = f (*a, b, c, d);
+  return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
                            gsl_matrix *subs[], size_t n_subs,
                            matrix_proto_m_mm *f)
 {
@@ -2727,19 +3360,19 @@ matrix_expr_evaluate_m_mm (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_v (const char *name,
+matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_v *f)
 {
   assert (n_subs == 1);
-  if (!check_vector_arg (name, subs, 0))
+  if (!check_vector_arg (props->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,
+matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_d_m *f)
 {
@@ -2751,7 +3384,7 @@ matrix_expr_evaluate_d_m (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_m_any (const char *name UNUSED,
+matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
                           gsl_matrix *subs[], size_t n_subs,
                           matrix_proto_m_any *f)
 {
@@ -2759,14 +3392,14 @@ matrix_expr_evaluate_m_any (const char *name UNUSED,
 }
 
 static gsl_matrix *
-matrix_expr_evaluate_IDENT (const char *name,
+matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
                             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))
+  if (!to_scalar_args (props->name, subs, n_subs, d))
     return NULL;
 
   return f (d[0], d[n_subs - 1]);
@@ -2825,11 +3458,16 @@ matrix_expr_evaluate (const struct matrix_expr *e)
   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); \
+#define F(ENUM, STRING, PROTO, CONSTRAINTS)                             \
+      case MOP_F_##ENUM:                                                \
+        {                                                               \
+          static const struct matrix_function_properties props = {      \
+            .name = STRING,                                             \
+            .constraints = CONSTRAINTS,                                 \
+          };                                                            \
+          result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
+                                                 matrix_eval_##ENUM);   \
+        }                                                               \
       break;
       MATRIX_FUNCTIONS
 #undef F
@@ -2933,7 +3571,7 @@ matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
   if (!is_scalar (m))
     {
       msg (SE, _("Expression for %s must evaluate to scalar, "
-                 "not a matrix with dimensions (%zu,%zu)."),
+                 "not a %zu×%zu matrix."),
            context, m->size1, m->size2);
       gsl_matrix_free (m);
       return false;
@@ -2973,8 +3611,11 @@ 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]);
+    {
+      for (size_t i = 0; i < lvalue->n_indexes; i++)
+        matrix_expr_destroy (lvalue->indexes[i]);
+      free (lvalue);
+    }
 }
 
 static struct matrix_lvalue *
@@ -2989,6 +3630,7 @@ matrix_lvalue_parse (struct matrix_state *s)
       if (!lvalue->var)
         {
           msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
+          free (lvalue);
           return NULL;
         }
 
@@ -3035,7 +3677,10 @@ matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
   else
     m = NULL;
 
-  return matrix_normalize_index_vector (m, size, index_type, other_size, iv);
+  bool ok = matrix_normalize_index_vector (m, size, index_type,
+                                           other_size, iv);
+  gsl_matrix_free (m);
+  return ok;
 }
 
 static bool
@@ -3047,15 +3692,15 @@ matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
   /* Convert source matrix 'sm' to source vector 'sv'. */
   if (!is_vector (sm))
     {
-      msg (SE, _("Can't assign matrix with dimensions (%zu,%zu) to subvector."),
+      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 vector with %zu elements "
-                 "to subvector with %zu."), sv.size, iv->n);
+      msg (SE, _("Can't assign %zu-element vector "
+                 "to %zu-element subvector."), sv.size, iv->n);
       return false;
     }
 
@@ -3144,7 +3789,7 @@ matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
     }
   else if (dm->size1 == 0 || dm->size2 == 0)
     {
-      msg (SE, _("Cannot index matrix with dimensions (%zu,%zu)."),
+      msg (SE, _("Cannot index %zu×%zu matrix."),
            dm->size1, dm->size2);
       return false;
     }
@@ -3152,9 +3797,8 @@ matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
     {
       if (!is_vector (dm))
         {
-          msg (SE, _("Can't use vector indexing on matrix %s with "
-                     "dimensions (%zu,%zu)."),
-               lvalue->var->name, dm->size1, dm->size2);
+          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],
@@ -3183,8 +3827,13 @@ static bool
 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
 {
   struct index_vector iv0, iv1;
-  return (matrix_lvalue_evaluate (lvalue, &iv0, &iv1)
-          && matrix_lvalue_assign (lvalue, &iv0, &iv1, sm));
+  if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
+    {
+      gsl_matrix_free (sm);
+      return false;
+    }
+
+  return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
 }
 
 
@@ -3196,6 +3845,8 @@ struct matrix_cmds
     size_t n;
   };
 
+static void matrix_cmds_uninit (struct matrix_cmds *);
+
 struct matrix_cmd
   {
     enum matrix_cmd_type
@@ -3274,13 +3925,23 @@ struct matrix_cmd
           }
         loop;
 
+        struct display_command
+          {
+            struct matrix_state *state;
+          }
+        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;
+            struct save_file *sf;
           }
         save;
 
@@ -3292,7 +3953,6 @@ struct matrix_cmd
             int c1, c2;
             enum fmt_type format;
             int w;
-            int d;
             bool symmetric;
             bool reread;
           }
@@ -3300,39 +3960,29 @@ struct matrix_cmd
 
         struct write_command
           {
+            struct write_file *wf;
             struct matrix_expr *expression;
-            struct file_handle *outfile;
-            char *encoding;
             int c1, c2;
-            enum fmt_type format;
-            int w;
-            int d;
+
+            /* If this is nonnull, WRITE uses this format.
+
+               If this is NULL, WRITE uses free-field format with as many
+               digits of precision as needed. */
+            struct fmt_spec *format;
+
             bool triangular;
             bool hold;
           }
         write;
 
-        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 get_command
           {
             struct matrix_lvalue *dst;
+            struct dataset *dataset;
             struct file_handle *file;
             char *encoding;
-            struct string_array variables;
+            struct var_syntax *vars;
+            size_t n_vars;
             struct matrix_var *names;
 
             /* Treatment of missing values. */
@@ -3355,11 +4005,10 @@ struct matrix_cmd
         struct msave_command
           {
             struct msave_common *common;
-            char *varname_;
             struct matrix_expr *expr;
             const char *rowtype;
-            struct matrix_expr *factors;
-            struct matrix_expr *splits;
+            const struct matrix_expr *factors;
+            const struct matrix_expr *splits;
           }
          msave;
 
@@ -3367,6 +4016,7 @@ struct matrix_cmd
           {
             struct matrix_state *state;
             struct file_handle *file;
+            char *encoding;
             struct stringi_set rowtypes;
           }
         mget;
@@ -3399,13 +4049,8 @@ struct matrix_cmd
 
 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 *);
 
-static void
-matrix_cmd_destroy (struct matrix_cmd *cmd)
-{
-  /* XXX */
-  free (cmd);
-}
 
 static struct matrix_cmd *
 matrix_parse_compute (struct matrix_state *s)
@@ -3777,7 +4422,9 @@ matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
     }
 
   string_array_destroy (rlabels);
+  free (rlabels);
   string_array_destroy (clabels);
+  free (clabels);
 }
 
 static void
@@ -3796,6 +4443,7 @@ create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
   if (!labels)
     d->hide_all_labels = true;
   string_array_destroy (labels);
+  free (labels);
 }
 
 static void
@@ -4130,32 +4778,18 @@ matrix_parse_break (struct matrix_state *s)
 static struct matrix_cmd *
 matrix_parse_display (struct matrix_state *s)
 {
-  bool dictionary = false;
-  bool status = false;
-  while (lex_token (s->lexer) == T_ID)
+  while (lex_token (s->lexer) != T_ENDCMD)
     {
-      if (lex_match_id (s->lexer, "DICTIONARY"))
-        dictionary = true;
-      else if (lex_match_id (s->lexer, "STATUS"))
-        status = true;
-      else
+      if (!lex_match_id (s->lexer, "DICTIONARY")
+          && !lex_match_id (s->lexer, "STATUS"))
         {
           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,
-    }
-  };
+  *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
   return cmd;
 }
 
@@ -4175,9 +4809,11 @@ 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)"));
+  struct pivot_dimension *attr_dimension
+    = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
+  pivot_category_create_group (attr_dimension->root, N_("Dimension"),
+                               N_("Rows"), N_("Columns"));
+  pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
 
   const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
   size_t n_vars = 0;
@@ -4202,6 +4838,7 @@ matrix_cmd_execute_display (struct display_command *cmd)
       for (size_t j = 0; j < sizeof values / sizeof *values; j++)
         pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
     }
+  free (vars);
   pivot_table_submit (table);
 }
 
@@ -4244,18 +4881,199 @@ matrix_cmd_execute_release (struct release_command *cmd)
     }
 }
 
+static struct save_file *
+save_file_create (struct matrix_state *s, struct file_handle *fh,
+                  struct string_array *variables,
+                  struct matrix_expr *names,
+                  struct stringi_set *strings)
+{
+  for (size_t i = 0; i < s->n_save_files; i++)
+    {
+      struct save_file *sf = s->save_files[i];
+      if (fh_equal (sf->file, fh))
+        {
+          fh_unref (fh);
+
+          string_array_destroy (variables);
+          matrix_expr_destroy (names);
+          stringi_set_destroy (strings);
+
+          return sf;
+        }
+    }
+
+  struct save_file *sf = xmalloc (sizeof *sf);
+  *sf = (struct save_file) {
+    .file = fh,
+    .dataset = s->dataset,
+    .variables = *variables,
+    .names = names,
+    .strings = STRINGI_SET_INITIALIZER (sf->strings),
+  };
+
+  stringi_set_swap (&sf->strings, strings);
+  stringi_set_destroy (strings);
+
+  s->save_files = xrealloc (s->save_files,
+                             (s->n_save_files + 1) * sizeof *s->save_files);
+  s->save_files[s->n_save_files++] = sf;
+  return sf;
+}
+
+static struct casewriter *
+save_file_open (struct save_file *sf, gsl_matrix *m)
+{
+  if (sf->writer || sf->error)
+    {
+      if (sf->writer)
+        {
+          size_t n_variables = caseproto_get_n_widths (
+            casewriter_get_proto (sf->writer));
+          if (m->size2 != n_variables)
+            {
+              msg (ME, _("The first SAVE to %s within this matrix program "
+                         "had %zu columns, so a %zu×%zu matrix cannot be "
+                         "saved to it."),
+                   sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
+                   n_variables, m->size1, m->size2);
+              return NULL;
+            }
+        }
+      return sf->writer;
+    }
+
+  bool ok = true;
+  struct dictionary *dict = dict_create (get_default_encoding ());
+
+  /* Fill 'names' with user-specified names if there were any, either from
+     sf->variables or sf->names. */
+  struct string_array names = { .n = 0 };
+  if (sf->variables.n)
+    string_array_clone (&names, &sf->variables);
+  else if (sf->names)
+    {
+      gsl_matrix *nm = matrix_expr_evaluate (sf->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;
+            }
+        }
+      gsl_matrix_free (nm);
+    }
+
+  struct stringi_set strings;
+  stringi_set_clone (&strings, &sf->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 (ME, _("Duplicate variable name %s in SAVE statement."),
+               name);
+          ok = false;
+        }
+    }
+
+  if (!stringi_set_is_empty (&strings))
+    {
+      size_t count = stringi_set_count (&strings);
+      const char *example = stringi_set_node_get_string (
+        stringi_set_first (&strings));
+
+      if (count == 1)
+        msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
+                   "variable %s."), example);
+      else
+        msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
+                           "unknown variable: %s.",
+                           "The SAVE command STRINGS subcommand specifies %zu "
+                           "unknown variables, including %s.",
+                           count),
+             count, example);
+      ok = false;
+    }
+  stringi_set_destroy (&strings);
+  string_array_destroy (&names);
+
+  if (!ok)
+    {
+      dict_unref (dict);
+      sf->error = true;
+      return NULL;
+    }
+
+  if (sf->file == fh_inline_file ())
+    sf->writer = autopaging_writer_create (dict_get_proto (dict));
+  else
+    sf->writer = any_writer_open (sf->file, dict);
+  if (sf->writer)
+    sf->dict = dict;
+  else
+    {
+      dict_unref (dict);
+      sf->error = true;
+    }
+  return sf->writer;
+}
+
+static void
+save_file_destroy (struct save_file *sf)
+{
+  if (sf)
+    {
+      if (sf->file == fh_inline_file () && sf->writer && sf->dict)
+        {
+          dataset_set_dict (sf->dataset, sf->dict);
+          dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
+        }
+      else
+        {
+          casewriter_destroy (sf->writer);
+          dict_unref (sf->dict);
+        }
+      fh_unref (sf->file);
+      string_array_destroy (&sf->variables);
+      matrix_expr_destroy (sf->names);
+      stringi_set_destroy (&sf->strings);
+      free (sf);
+    }
+}
+
 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)
-    }
+    .save = { .expression = NULL },
   };
 
+  struct file_handle *fh = NULL;
   struct save_command *save = &cmd->save;
+
+  struct string_array variables = STRING_ARRAY_INITIALIZER;
+  struct matrix_expr *names = NULL;
+  struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
+
   save->expression = matrix_parse_exp (s);
   if (!save->expression)
     goto error;
@@ -4266,11 +5084,11 @@ matrix_parse_save (struct matrix_state *s)
         {
           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)
+          fh_unref (fh);
+          fh = (lex_match (s->lexer, T_ASTERISK)
+                ? fh_inline_file ()
+                : fh_parse (s->lexer, FH_REF_FILE, s->session));
+          if (!fh)
             goto error;
         }
       else if (lex_match_id (s->lexer, "VARIABLES"))
@@ -4286,10 +5104,8 @@ matrix_parse_save (struct matrix_state *s)
           if (!ok)
             goto error;
 
-          string_array_destroy (save->variables);
-          if (!save->variables)
-            save->variables = xmalloc (sizeof *save->variables);
-          *save->variables = (struct string_array) {
+          string_array_clear (&variables);
+          variables = (struct string_array) {
             .strings = names,
             .n = n,
             .allocated = n,
@@ -4298,9 +5114,9 @@ matrix_parse_save (struct matrix_state *s)
       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)
+          matrix_expr_destroy (names);
+          names = matrix_parse_exp (s);
+          if (!names)
             goto error;
         }
       else if (lex_match_id (s->lexer, "STRINGS"))
@@ -4308,7 +5124,7 @@ matrix_parse_save (struct matrix_state *s)
           lex_match (s->lexer, T_EQUALS);
           while (lex_token (s->lexer) == T_ID)
             {
-              stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
+              stringi_set_insert (&strings, lex_tokcstr (s->lexer));
               lex_get (s->lexer);
               if (!lex_match (s->lexer, T_COMMA))
                 break;
@@ -4322,29 +5138,34 @@ matrix_parse_save (struct matrix_state *s)
         }
     }
 
-  if (!save->outfile)
+  if (!fh)
     {
-      if (s->prev_save_outfile)
-        save->outfile = fh_ref (s->prev_save_outfile);
+      if (s->prev_save_file)
+        fh = fh_ref (s->prev_save_file);
       else
         {
           lex_sbc_missing ("OUTFILE");
           goto error;
         }
     }
-  fh_unref (s->prev_save_outfile);
-  s->prev_save_outfile = fh_ref (save->outfile);
+  fh_unref (s->prev_save_file);
+  s->prev_save_file = fh_ref (fh);
 
-  if (save->variables && save->names)
+  if (variables.n && names)
     {
       msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
-      matrix_expr_destroy (save->names);
-      save->names = NULL;
+      matrix_expr_destroy (names);
+      names = NULL;
     }
 
+  save->sf = save_file_create (s, fh, &variables, names, &strings);
   return cmd;
 
 error:
+  string_array_destroy (&variables);
+  matrix_expr_destroy (names);
+  stringi_set_destroy (&strings);
+  fh_unref (fh);
   matrix_cmd_destroy (cmd);
   return NULL;
 }
@@ -4352,90 +5173,18 @@ error:
 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)
+  struct casewriter *writer = save_file_open (save->sf, m);
+  if (!writer)
     {
-      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;
-            }
-        }
+      gsl_matrix_free (m);
+      return;
     }
 
-  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);
+  const struct caseproto *proto = casewriter_get_proto (writer);
   for (size_t y = 0; y < m->size1; y++)
     {
       struct ccase *c = case_create (proto);
@@ -4451,10 +5200,7 @@ matrix_cmd_execute_save (const struct save_command *save)
         }
       casewriter_write (writer, c);
     }
-  casewriter_destroy (writer);
-
   gsl_matrix_free (m);
-  dict_unref (dict);
 }
 
 static struct matrix_cmd *
@@ -4535,6 +5281,7 @@ matrix_parse_read (struct matrix_state *s)
       else if (lex_match_id (s->lexer, "SIZE"))
         {
           lex_match (s->lexer, T_EQUALS);
+          matrix_expr_destroy (read->size);
           read->size = matrix_parse_exp (s);
           if (!read->size)
             goto error;
@@ -4580,14 +5327,15 @@ matrix_parse_read (struct matrix_state *s)
                 }
               lex_get (s->lexer);
             }
-          else if (!fmt_from_name (p, &read->format))
+          else if (fmt_from_name (p, &read->format))
+            lex_get (s->lexer);
+          else
             {
               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
@@ -4625,6 +5373,7 @@ matrix_parse_read (struct matrix_state *s)
   s->prev_read_file = fh_ref (fh);
 
   read->rf = read_file_create (s, fh);
+  fh = NULL;
   if (encoding)
     {
       free (read->rf->encoding);
@@ -4668,33 +5417,84 @@ matrix_parse_read (struct matrix_state *s)
   return cmd;
 
 error:
+  fh_unref (fh);
   matrix_cmd_destroy (cmd);
   free (encoding);
   return NULL;
 }
 
+static void
+parse_error (const struct dfm_reader *reader, enum fmt_type format,
+             struct substring data, size_t y, size_t x,
+             int first_column, int last_column, char *error)
+{
+  int line_number = dfm_get_line_number (reader);
+  struct msg_location *location = xmalloc (sizeof *location);
+  *location = (struct msg_location) {
+    .file_name = xstrdup (dfm_get_file_name (reader)),
+    .first_line = line_number,
+    .last_line = line_number + 1,
+    .first_column = first_column,
+    .last_column = last_column,
+  };
+  struct msg *m = xmalloc (sizeof *m);
+  *m = (struct msg) {
+    .category = MSG_C_DATA,
+    .severity = MSG_S_WARNING,
+    .location = location,
+    .text = xasprintf (_("Error reading \"%.*s\" as format %s "
+                         "for matrix row %zu, column %zu: %s"),
+                       (int) data.length, data.string, fmt_name (format),
+                       y + 1, x + 1, error),
+  };
+  msg_emit (m);
+
+  free (error);
+}
+
 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)
+                       gsl_matrix *m, struct substring p, size_t y, size_t x,
+                       const char *line_start)
 {
   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 */
+  char *error;
+  double f;
+  if (fmt_is_numeric (read->format))
+    {
+      union value v;
+      error = data_in (p, input_encoding, read->format,
+                       settings_get_fmt_settings (), &v, 0, NULL);
+      if (!error && v.f == SYSMIS)
+        error = xstrdup (_("Matrix data may not contain missing value."));
+      f = v.f;
+    }
+    else
+      {
+        uint8_t s[sizeof (double)];
+        union value v = { .s = s };
+        error = data_in (p, input_encoding, read->format,
+                         settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
+        memcpy (&f, s, sizeof f);
+      }
+
   if (error)
-    msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
+    {
+      int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
+      int c2 = c1 + ss_utf8_count_columns (p) - 1;
+      parse_error (reader, read->format, p, y, x, c1, c2, error);
+    }
   else
     {
-      gsl_matrix_set (m, y, x, v.f);
+      gsl_matrix_set (m, y, x, f);
       if (read->symmetric && x != y)
-        gsl_matrix_set (m, x, y, v.f);
+        gsl_matrix_set (m, x, y, f);
     }
 }
 
 static bool
 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
-                  struct substring *line)
+                  struct substring *line, const char **startp)
 {
   if (dfm_eof (reader))
     {
@@ -4702,8 +5502,10 @@ matrix_read_line (struct read_command *read, struct dfm_reader *reader,
       return false;
     }
   dfm_expand_tabs (reader);
-  *line = ss_substr (dfm_get_record (reader),
-                     read->c1 - 1, read->c2 - read->c1);
+  struct substring record = dfm_get_record (reader);
+  /* XXX need to recode record into UTF-8 */
+  *startp = record.string;
+  *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
   return true;
 }
 
@@ -4716,6 +5518,7 @@ matrix_read (struct read_command *read, struct dfm_reader *reader,
       size_t nx = read->symmetric ? y + 1 : m->size2;
 
       struct substring line = ss_empty ();
+      const char *line_start = line.string;
       for (size_t x = 0; x < nx; x++)
         {
           struct substring p;
@@ -4726,7 +5529,7 @@ matrix_read (struct read_command *read, struct dfm_reader *reader,
                   ss_ltrim (&line, ss_cstr (" ,"));
                   if (!ss_is_empty (line))
                     break;
-                  if (!matrix_read_line (read, reader, &line))
+                  if (!matrix_read_line (read, reader, &line, &line_start))
                     return;
                   dfm_forward_record (reader);
                 }
@@ -4735,7 +5538,7 @@ matrix_read (struct read_command *read, struct dfm_reader *reader,
             }
           else
             {
-              if (!matrix_read_line (read, reader, &line))
+              if (!matrix_read_line (read, reader, &line, &line_start))
                 return;
               size_t fields_per_line = (read->c2 - read->c1) / read->w;
               int f = x % fields_per_line;
@@ -4745,7 +5548,7 @@ matrix_read (struct read_command *read, struct dfm_reader *reader,
               p = ss_substr (line, read->w * f, read->w);
             }
 
-          matrix_read_set_field (read, reader, m, p, y, x);
+          matrix_read_set_field (read, reader, m, p, y, x, line_start);
         }
 
       if (read->w)
@@ -4754,8 +5557,11 @@ matrix_read (struct read_command *read, struct dfm_reader *reader,
         {
           ss_ltrim (&line, ss_cstr (" ,"));
           if (!ss_is_empty (line))
-            msg (SW, _("Trailing garbage on line \"%.*s\""),
-                 (int) line.length, line.string);
+            {
+              /* XXX */
+              msg (SW, _("Trailing garbage on line \"%.*s\""),
+                   (int) line.length, line.string);
+            }
         }
     }
 }
@@ -4776,7 +5582,8 @@ matrix_cmd_execute_read (struct read_command *read)
 
       if (!is_vector (m))
         {
-          msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
+          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);
@@ -4797,7 +5604,9 @@ matrix_cmd_execute_read (struct read_command *read)
         }
       else
         {
-          msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
+          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);
@@ -4807,7 +5616,9 @@ matrix_cmd_execute_read (struct read_command *read)
 
       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]);
+          msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
+                     "are outside valid range."),
+               d[0], d[1]);
           free (iv0.indexes);
           free (iv1.indexes);
           return;
@@ -4840,8 +5651,8 @@ matrix_cmd_execute_read (struct read_command *read)
         {
           if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
             {
-              msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
-                         "(%zu,%zu)."),
+              msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
+                         "differ from submatrix dimensions %zu×%zu."),
                    size[0], size[1],
                    submatrix_size[0], submatrix_size[1]);
               free (iv0.indexes);
@@ -4862,7 +5673,7 @@ matrix_cmd_execute_read (struct read_command *read)
 
   if (read->symmetric && size[0] != size[1])
     {
-      msg (SE, _("Cannot read matrix with non-square dimensions (%zu,%zu) "
+      msg (SE, _("Cannot read non-square %zu×%zu matrix "
                  "using READ with MODE=SYMMETRIC."),
            size[0], size[1]);
       free (iv0.indexes);
@@ -4880,9 +5691,10 @@ 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)
@@ -4891,16 +5703,17 @@ matrix_parse_write (struct matrix_state *s)
   int by = 0;
   int repetitions = 0;
   int record_width = 0;
-  bool seen_format = false;
+  enum fmt_type format = FMT_F;
+  bool has_format = false;
   while (lex_match (s->lexer, T_SLASH))
     {
       if (lex_match_id (s->lexer, "OUTFILE"))
         {
           lex_match (s->lexer, T_EQUALS);
 
-          fh_unref (write->outfile);
-          write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
-          if (!write->outfile)
+          fh_unref (fh);
+          fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
+          if (!fh)
             goto error;
         }
       else if (lex_match_id (s->lexer, "ENCODING"))
@@ -4909,8 +5722,8 @@ matrix_parse_write (struct matrix_state *s)
 	  if (!lex_force_string (s->lexer))
 	    goto error;
 
-          free (write->encoding);
-          write->encoding = ss_xstrdup (lex_tokss (s->lexer));
+          free (encoding);
+          encoding = ss_xstrdup (lex_tokss (s->lexer));
 
 	  lex_get (s->lexer);
 	}
@@ -4964,12 +5777,11 @@ matrix_parse_write (struct matrix_state *s)
         write->hold = true;
       else if (lex_match_id (s->lexer, "FORMAT"))
         {
-          if (seen_format)
+          if (has_format || write->format)
             {
               lex_sbc_only_once ("FORMAT");
               goto error;
             }
-          seen_format = true;
 
           lex_match (s->lexer, T_EQUALS);
 
@@ -4981,21 +5793,25 @@ matrix_parse_write (struct matrix_state *s)
             {
               repetitions = atoi (p);
               p += strspn (p, "0123456789");
-              if (!fmt_from_name (p, &write->format))
+              if (!fmt_from_name (p, &format))
                 {
                   lex_error (s->lexer, _("Unknown format %s."), p);
                   goto error;
                 }
+              has_format = true;
               lex_get (s->lexer);
             }
-          else if (!fmt_from_name (p, &write->format))
+          else if (fmt_from_name (p, &format))
             {
-              struct fmt_spec format;
-              if (!parse_format_specifier (s->lexer, &format))
+              has_format = true;
+              lex_get (s->lexer);
+            }
+          else
+            {
+              struct fmt_spec spec;
+              if (!parse_format_specifier (s->lexer, &spec))
                 goto error;
-              write->format = format.type;
-              write->w = format.w;
-              write->d = format.d;
+              write->format = xmemdup (&spec, sizeof spec);
             }
         }
       else
@@ -5012,18 +5828,27 @@ matrix_parse_write (struct matrix_state *s)
       goto error;
     }
 
-  if (!write->outfile)
+  if (!fh)
     {
-      if (s->prev_write_outfile)
-        write->outfile = fh_ref (s->prev_write_outfile);
+      if (s->prev_write_file)
+        fh = fh_ref (s->prev_write_file);
       else
         {
           lex_sbc_missing ("OUTFILE");
           goto error;
         }
     }
-  fh_unref (s->prev_write_outfile);
-  s->prev_write_outfile = fh_ref (write->outfile);
+  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:
 
@@ -5043,7 +5868,7 @@ matrix_parse_write (struct matrix_state *s)
       goto error;
     }
   int w = (repetitions ? record_width / repetitions
-           : write->w ? write->w
+           : write->format ? write->format->w
            : by);
   if (by && w != by)
     {
@@ -5057,10 +5882,28 @@ matrix_parse_write (struct matrix_state *s)
              w, by);
       goto error;
     }
-  write->w = w;
+  if (w && !write->format)
+    {
+      write->format = xmalloc (sizeof *write->format);
+      *write->format = (struct fmt_spec) { .type = format, .w = w };
+
+      if (!fmt_check_output (write->format))
+        goto error;
+    };
+
+  if (write->format && fmt_var_width (write->format) > sizeof (double))
+    {
+      char s[FMT_STRING_LEN_MAX + 1];
+      fmt_to_string (write->format, s);
+      msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
+           s, sizeof (double));
+      goto error;
+    }
+
   return cmd;
 
 error:
+  fh_unref (fh);
   matrix_cmd_destroy (cmd);
   return NULL;
 }
@@ -5075,56 +5918,77 @@ matrix_cmd_execute_write (struct write_command *write)
   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)."),
+                 "the matrix to be written has dimensions %zu×%zu."),
            m->size1, m->size2);
       gsl_matrix_free (m);
       return;
     }
 
-  struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
-  if (!writer)
+  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 = U8_LINE_INITIALIZER;
+  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));
+          char *s;
+          double f = gsl_matrix_get (m, y, x);
+          if (write->format)
+            {
+              union value v;
+              if (fmt_is_numeric (write->format->type))
+                v.f = f;
+              else
+                v.s = (uint8_t *) &f;
+              s = data_out (&v, NULL, write->format, settings);
+            }
+          else
+            {
+              s = xmalloc (DBL_BUFSIZE_BOUND);
+              if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
+                  >= DBL_BUFSIZE_BOUND)
+                abort ();
+            }
           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);
+              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);
+          u8_line_put (line, x0, x0 + width, s, len);
           free (s);
 
-          x0 += write->w ? write->w : width + 1;
+          x0 += write->format ? write->format->w : width + 1;
         }
 
-      dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
-      u8_line_clear (&line);
+      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);
     }
-  u8_line_destroy (&line);
-  dfm_close_writer (writer);
+  if (!write->hold)
+    {
+      u8_line_destroy (line);
+      free (line);
+      line = NULL;
+    }
+  write->wf->held = line;
 
   gsl_matrix_free (m);
 }
@@ -5136,6 +6000,7 @@ matrix_parse_get (struct matrix_state *s)
   *cmd = (struct matrix_cmd) {
     .type = MCMD_GET,
     .get = {
+      .dataset = s->dataset,
       .user = { .treatment = MGET_ERROR },
       .system = { .treatment = MGET_ERROR },
     }
@@ -5150,25 +6015,20 @@ matrix_parse_get (struct matrix_state *s)
     {
       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;
+          if (lex_match (s->lexer, T_ASTERISK))
+            get->file = NULL;
+          else
+            {
+              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;
@@ -5182,40 +6042,14 @@ matrix_parse_get (struct matrix_state *s)
         {
           lex_match (s->lexer, T_EQUALS);
 
-          struct dictionary *dict = NULL;
-          if (!get->file)
+          if (get->n_vars)
             {
-              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);
+              lex_sbc_only_once ("VARIABLES");
               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);
+          if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
+            goto error;
         }
       else if (lex_match_id (s->lexer, "NAMES"))
         {
@@ -5252,11 +6086,11 @@ matrix_parse_get (struct matrix_state *s)
         {
 	  lex_match (s->lexer, T_EQUALS);
           if (lex_match_id (s->lexer, "OMIT"))
-            get->user.treatment = MGET_OMIT;
+            get->system.treatment = MGET_OMIT;
           else if (lex_is_number (s->lexer))
             {
-              get->user.treatment = MGET_RECODE;
-              get->user.substitute = lex_number (s->lexer);
+              get->system.treatment = MGET_RECODE;
+              get->system.substitute = lex_number (s->lexer);
               lex_get (s->lexer);
             }
           else
@@ -5272,6 +6106,10 @@ matrix_parse_get (struct matrix_state *s)
           goto error;
         }
     }
+
+  if (get->user.treatment != MGET_ACCEPT)
+    get->system.treatment = MGET_ERROR;
+
   return cmd;
 
 error:
@@ -5280,62 +6118,52 @@ error:
 }
 
 static void
-matrix_cmd_execute_get (struct get_command *get)
+matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
+                          const struct dictionary *dict)
 {
-  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);
+  struct variable **vars;
   size_t n_vars = 0;
 
-  if (get->variables.n)
+  if (get->n_vars)
     {
-      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;
-        }
+      if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
+                                &vars, &n_vars, PV_NUMERIC))
+        return;
     }
   else
     {
-      for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
+      n_vars = dict_get_var_cnt (dict);
+      vars = xnmalloc (n_vars, sizeof *vars);
+      for (size_t i = 0; i < n_vars; i++)
         {
-          const struct variable *var = dict_get_var (dict, i);
+          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;
+          vars[i] = var;
         }
     }
 
+  if (get->names)
+    {
+      gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
+      for (size_t i = 0; i < n_vars; i++)
+        {
+          char s[sizeof (double)];
+          double f;
+          buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
+          memcpy (&f, s, sizeof f);
+          gsl_matrix_set (names, i, 0, f);
+        }
+
+      gsl_matrix_free (get->names->value);
+      get->names->value = names;
+    }
+
   size_t n_rows = 0;
   gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
   long long int casenum = 1;
@@ -5394,7 +6222,6 @@ matrix_cmd_execute_get (struct get_command *get)
       if (keep)
         n_rows++;
     }
-  casereader_destroy (reader);
   if (!error)
     {
       m->size1 = n_rows;
@@ -5402,9 +6229,56 @@ matrix_cmd_execute_get (struct get_command *get)
     }
   else
     gsl_matrix_free (m);
-  dict_unref (dict);
   free (vars);
 }
+
+static bool
+matrix_open_casereader (const char *command_name,
+                        struct file_handle *file, const char *encoding,
+                        struct dataset *dataset,
+                        struct casereader **readerp, struct dictionary **dictp)
+{
+  if (file)
+    {
+       *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
+       return *readerp != NULL;
+    }
+  else
+    {
+      if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
+        {
+          msg (ME, _("The %s command cannot read an empty active file."),
+               command_name);
+          return false;
+        }
+      *readerp = proc_open (dataset);
+      *dictp = dict_ref (dataset_dict (dataset));
+      return true;
+    }
+}
+
+static void
+matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
+                         struct casereader *reader, struct dictionary *dict)
+{
+  dict_unref (dict);
+  casereader_destroy (reader);
+  if (!file)
+    proc_commit (dataset);
+}
+
+static void
+matrix_cmd_execute_get (struct get_command *get)
+{
+  struct casereader *r;
+  struct dictionary *d;
+  if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
+                              &r, &d))
+    {
+      matrix_cmd_execute_get__ (get, r, d);
+      matrix_close_casereader (get->file, get->dataset, r, d);
+    }
+}
 
 static const char *
 match_rowtype (struct lexer *lexer)
@@ -5430,8 +6304,6 @@ parse_var_names (struct lexer *lexer, struct string_array *sa)
   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,
@@ -5440,6 +6312,17 @@ parse_var_names (struct lexer *lexer, struct string_array *sa)
 
   if (ok)
     {
+      for (size_t i = 0; i < n_names; i++)
+        if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
+            || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
+          {
+            msg (SE, _("Variable name %s is reserved."), names[i]);
+            for (size_t j = 0; j < n_names; j++)
+              free (names[i]);
+            free (names);
+            return false;
+          }
+
       string_array_clear (sa);
       sa->strings = names;
       sa->n = sa->allocated = n_names;
@@ -5448,7 +6331,7 @@ parse_var_names (struct lexer *lexer, struct string_array *sa)
 }
 
 static void
-msave_common_uninit (struct msave_common *common)
+msave_common_destroy (struct msave_common *common)
 {
   if (common)
     {
@@ -5456,6 +6339,19 @@ msave_common_uninit (struct msave_common *common)
       string_array_destroy (&common->variables);
       string_array_destroy (&common->fnames);
       string_array_destroy (&common->snames);
+
+      for (size_t i = 0; i < common->n_factors; i++)
+        matrix_expr_destroy (common->factors[i]);
+      free (common->factors);
+
+      for (size_t i = 0; i < common->n_splits; i++)
+        matrix_expr_destroy (common->splits[i]);
+      free (common->splits);
+
+      dict_unref (common->dict);
+      casewriter_destroy (common->writer);
+
+      free (common);
     }
 }
 
@@ -5485,16 +6381,19 @@ compare_variables (const char *keyword,
 static struct matrix_cmd *
 matrix_parse_msave (struct matrix_state *s)
 {
-  struct msave_common common = { .outfile = NULL };
+  struct msave_common *common = xmalloc (sizeof *common);
+  *common = (struct msave_common) { .outfile = NULL };
+
   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
   *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
 
+  struct matrix_expr *splits = NULL;
+  struct matrix_expr *factors = 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;
+    goto error;
 
   while (lex_match (s->lexer, T_SLASH))
     {
@@ -5510,42 +6409,42 @@ matrix_parse_msave (struct matrix_state *s)
         {
           lex_match (s->lexer, T_EQUALS);
 
-          fh_unref (common.outfile);
-          common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
-          if (!common.outfile)
+          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))
+          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))
+          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))
+          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)
+          matrix_expr_destroy (splits);
+          splits = matrix_parse_exp (s);
+          if (!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)
+          matrix_expr_destroy (factors);
+          factors = matrix_parse_exp (s);
+          if (!factors)
             goto error;
         }
       else
@@ -5560,49 +6459,31 @@ matrix_parse_msave (struct matrix_state *s)
       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)
+      if (common->fnames.n && !factors)
+        {
+          msg (SE, _("FNAMES requires FACTOR."));
+          goto error;
+        }
+      if (common->snames.n && !splits)
+        {
+          msg (SE, _("SNAMES requires SPLIT."));
+          goto error;
+        }
+      if (!common->outfile)
         {
           lex_sbc_missing ("OUTFILE");
           goto error;
         }
-      s->common = xmemdup (&common, sizeof common);
+      s->common = common;
     }
   else
     {
-      if (common.outfile)
+      if (common->outfile)
         {
-          bool same = common.outfile == s->common->outfile;
-          fh_unref (common.outfile);
-          if (!same)
+          if (!fh_equal (common->outfile, s->common->outfile))
             {
               msg (SE, _("OUTFILE must name the same file on each MSAVE "
                          "within a single MATRIX command."));
@@ -5610,25 +6491,47 @@ matrix_parse_msave (struct matrix_state *s)
             }
         }
       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))
+                              &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_destroy (common);
     }
   msave->common = s->common;
-  if (!msave->varname_)
-    msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
+
+  struct msave_common *c = s->common;
+  if (factors)
+    {
+      if (c->n_factors >= c->allocated_factors)
+        c->factors = x2nrealloc (c->factors, &c->allocated_factors,
+                                 sizeof *c->factors);
+      c->factors[c->n_factors++] = factors;
+    }
+  if (c->n_factors > 0)
+    msave->factors = c->factors[c->n_factors - 1];
+
+  if (splits)
+    {
+      if (c->n_splits >= c->allocated_splits)
+        c->splits = x2nrealloc (c->splits, &c->allocated_splits,
+                                sizeof *c->splits);
+      c->splits[c->n_splits++] = splits;
+    }
+  if (c->n_splits > 0)
+    msave->splits = c->splits[c->n_splits - 1];
+
   return cmd;
 
 error:
-  msave_common_uninit (&common);
+  matrix_expr_destroy (splits);
+  matrix_expr_destroy (factors);
+  msave_common_destroy (common);
   matrix_cmd_destroy (cmd);
   return NULL;
 }
 
 static gsl_vector *
-matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
+matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
 {
   gsl_matrix *m = matrix_expr_evaluate (e);
   if (!m)
@@ -5636,8 +6539,8 @@ matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
 
   if (!is_vector (m))
     {
-      msg (SE, _("%s expression must evaluate to vector, not a matrix with "
-                 "dimensions (%zu,%zu)."),
+      msg (SE, _("%s expression must evaluate to vector, "
+                 "not a %zu×%zu matrix."),
            name, m->size1, m->size2);
       gsl_matrix_free (m);
       return NULL;
@@ -5663,8 +6566,9 @@ msave_create_dict (const struct msave_common *common)
   const char *dup_split = msave_add_vars (dict, &common->snames);
   if (dup_split)
     {
-      msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
-      goto error;
+      /* Should not be possible because the parser ensures that the names are
+         unique. */
+      NOT_REACHED ();
     }
 
   dict_create_var_assert (dict, "ROWTYPE_", 8);
@@ -5708,10 +6612,10 @@ matrix_cmd_execute_msave (struct msave_command *msave)
     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)
+  else if (m->size2 != common->variables.n)
     {
-      msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
+      msg (SE,
+           _("Matrix on MSAVE has %zu columns but there are %zu variables."),
            m->size2, common->variables.n);
       goto error;
     }
@@ -5726,11 +6630,10 @@ matrix_cmd_execute_msave (struct msave_command *msave)
         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)
+      else if (factors->size != common->fnames.n)
         {
           msg (SE, _("There are %zu factor variables, "
-                     "but %zu split values were supplied."),
+                     "but %zu factor values were supplied."),
                common->fnames.n, factors->size);
           goto error;
         }
@@ -5742,16 +6645,15 @@ matrix_cmd_execute_msave (struct msave_command *msave)
       if (!splits)
         goto error;
 
-      if (!common->fnames.n)
+      if (!common->snames.n)
         for (size_t i = 0; i < splits->size; i++)
-          string_array_append_nocopy (&common->fnames,
+          string_array_append_nocopy (&common->snames,
                                       xasprintf ("SPL%zu", i + 1));
-
-      if (splits->size != common->fnames.n)
+      else if (splits->size != common->snames.n)
         {
           msg (SE, _("There are %zu split variables, "
                      "but %zu split values were supplied."),
-               common->fnames.n, splits->size);
+               common->snames.n, splits->size);
           goto error;
         }
     }
@@ -5772,6 +6674,8 @@ matrix_cmd_execute_msave (struct msave_command *msave)
       common->dict = dict;
     }
 
+  bool matrix = (!strcmp (msave->rowtype, "COV")
+                 || !strcmp (msave->rowtype, "CORR"));
   for (size_t y = 0; y < m->size1; y++)
     {
       struct ccase *c = case_create (dict_get_proto (common->dict));
@@ -5792,8 +6696,11 @@ matrix_cmd_execute_msave (struct msave_command *msave)
           *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
 
       /* VARNAME_. */
+      const char *varname_ = (matrix && y < common->variables.n
+                              ? common->variables.strings[y]
+                              : "");
       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
-                         msave->varname_, ' ');
+                         varname_, ' ');
 
       /* Continuous variables. */
       for (size_t x = 0; x < m->size2; x++)
@@ -5811,11 +6718,18 @@ 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 } };
+  *cmd = (struct matrix_cmd) {
+    .type = MCMD_MGET,
+    .mget = {
+      .state = s,
+      .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
+    },
+  };
 
   struct mget_command *mget = &cmd->mget;
 
-  for (;;)
+  lex_match (s->lexer, T_SLASH);
+  while (lex_token (s->lexer) != T_ENDCMD)
     {
       if (lex_match_id (s->lexer, "FILE"))
         {
@@ -5826,6 +6740,17 @@ matrix_parse_mget (struct matrix_state *s)
           if (!mget->file)
             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 (mget->encoding);
+          mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
+
+	  lex_get (s->lexer);
+	}
       else if (lex_match_id (s->lexer, "TYPE"))
         {
           lex_match (s->lexer, T_EQUALS);
@@ -5844,11 +6769,7 @@ matrix_parse_mget (struct matrix_state *s)
           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;
+      lex_match (s->lexer, T_SLASH);
     }
   return cmd;
 
@@ -5877,34 +6798,82 @@ get_a8_var (const struct dictionary *d, const char *name)
 }
 
 static bool
-is_rowtype (const union value *v, const char *rowtype)
+var_changed (const struct ccase *ca, const struct ccase *cb,
+             const struct variable *var)
 {
-  struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
-  ss_rtrim (&vs, ss_cstr (" "));
-  return ss_equals_case (vs, ss_cstr (rowtype));
+  return (ca && cb
+          ? !value_equal (case_data (ca, var), case_data (cb, var),
+                          var_get_width (var))
+          : ca || cb);
+}
+
+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 bool
+vars_all_missing (const struct ccase *c, const struct dictionary *d,
+                  size_t first_var, size_t n_vars)
+{
+  for (size_t i = 0; i < n_vars; i++)
+    if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
+      return false;
+  return true;
 }
 
 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)
+                        const struct stringi_set *accepted_rowtypes,
+                        struct matrix_state *s,
+                        size_t ss, size_t sn, size_t si,
+                        size_t fs, size_t fn, size_t fi,
+                        size_t cs, size_t cn,
+                        struct pivot_table *pt,
+                        struct pivot_dimension *var_dimension)
 {
   if (!n_rows)
-    return;
+    goto exit;
+
+  /* Is this a matrix for pooled data, either where there are no factor
+     variables or the factor variables are missing? */
+  bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
 
-  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 substring rowtype = case_ss (rows[0], rowtype_var);
+  ss_rtrim (&rowtype, ss_cstr (" "));
+  if (!stringi_set_is_empty (accepted_rowtypes)
+      && !stringi_set_contains_len (accepted_rowtypes,
+                                    rowtype.string, rowtype.length))
+    goto exit;
+
+  const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
+                        : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
+                        : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
+                        : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
+                        : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
+                        : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
+                        : NULL);
+  if (!prefix)
+    {
+      msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
+           (int) rowtype.length, rowtype.string);
+      goto exit;
+    }
 
   struct string name = DS_EMPTY_INITIALIZER;
-  ds_put_cstr (&name, name_prefix);
-  if (fi > 0)
+  ds_put_cstr (&name, prefix);
+  if (!pooled)
     ds_put_format (&name, "F%zu", fi);
   if (si > 0)
     ds_put_format (&name, "S%zu", si);
@@ -5916,7 +6885,7 @@ matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
     {
       msg (SW, _("Matrix data file contains variable with existing name %s."),
            ds_cstr (&name));
-      goto exit;
+      goto exit_free_name;
     }
 
   gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
@@ -5936,6 +6905,28 @@ matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
         }
     }
 
+  int var_index = pivot_category_create_leaf (
+    var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
+  double values[] = { n_rows, cn };
+  for (size_t j = 0; j < sn; j++)
+    {
+      struct variable *var = dict_get_var (d, ss + j);
+      const union value *value = case_data (rows[0], var);
+      pivot_table_put2 (pt, j, var_index,
+                        pivot_value_new_var_value (var, value));
+    }
+  for (size_t j = 0; j < fn; j++)
+    {
+      struct variable *var = dict_get_var (d, fs + j);
+      const union value sysmis = { .f = SYSMIS };
+      const union value *value = pooled ? &sysmis : case_data (rows[0], var);
+      pivot_table_put2 (pt, j + sn, var_index,
+                        pivot_value_new_var_value (var, value));
+    }
+  for (size_t j = 0; j < sizeof values / sizeof *values; j++)
+    pivot_table_put2 (pt, j + sn + fn, var_index,
+                      pivot_value_new_integer (values[j]));
+
   if (n_missing)
     msg (SE, ngettext ("Matrix data file variable %s contains a missing "
                        "value, which was treated as zero.",
@@ -5944,57 +6935,32 @@ matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
          ds_cstr (&name), n_missing);
   mv->value = m;
 
-exit:
+exit_free_name:
   ds_destroy (&name);
+
+exit:
   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)
+matrix_cmd_execute_mget__ (struct mget_command *mget,
+                           struct casereader *r, const struct dictionary *d)
 {
-  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;
+    return;
 
   if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
     {
       msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
-      goto exit;
+      return;
     }
   if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
     {
       msg (SE, _("Matrix data file contains no continuous variables."));
-      goto exit;
+      return;
     }
 
   for (size_t i = 0; i < dict_get_var_cnt (d); i++)
@@ -6005,7 +6971,7 @@ matrix_cmd_execute_mget (struct mget_command *mget)
           msg (SE,
                _("Matrix data file contains unexpected string variable %s."),
                var_get_name (v));
-          goto exit;
+          return;
         }
     }
 
@@ -6026,6 +6992,32 @@ matrix_cmd_execute_mget (struct mget_command *mget)
   size_t cn = dict_get_var_cnt (d) - cs;
   struct ccase *cc = NULL;
 
+  /* Pivot table. */
+  struct pivot_table *pt = pivot_table_create (
+    N_("Matrix Variables Created by MGET"));
+  struct pivot_dimension *attr_dimension = pivot_dimension_create (
+    pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
+  struct pivot_dimension *var_dimension = pivot_dimension_create (
+    pt, PIVOT_AXIS_ROW, N_("Variable"));
+  if (sn > 0)
+    {
+      struct pivot_category *splits = pivot_category_create_group (
+        attr_dimension->root, N_("Split Values"));
+      for (size_t i = 0; i < sn; i++)
+        pivot_category_create_leaf (splits, pivot_value_new_variable (
+                                      dict_get_var (d, ss + i)));
+    }
+  if (fn > 0)
+    {
+      struct pivot_category *factors = pivot_category_create_group (
+        attr_dimension->root, N_("Factors"));
+      for (size_t i = 0; i < fn; i++)
+        pivot_category_create_leaf (factors, pivot_value_new_variable (
+                                      dict_get_var (d, fs + i)));
+    }
+  pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
+                                N_("Rows"), N_("Columns"));
+
   /* Matrix. */
   struct ccase **rows = NULL;
   size_t allocated_rows = 0;
@@ -6034,25 +7026,30 @@ matrix_cmd_execute_mget (struct mget_command *mget)
   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)
+      bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
+
+      enum
         {
-          fi++;
-          case_unref (fc);
-          fc = case_ref (c);
+          SPLITS_CHANGED,
+          FACTORS_CHANGED,
+          ROWTYPE_CHANGED,
+          NOTHING_CHANGED
         }
-      if (md)
+      change
+        = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
+           : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
+           : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
+           : NOTHING_CHANGED);
+
+      if (change != NOTHING_CHANGED)
         {
-          matrix_mget_commit_var (rows, n_rows, d, rowtype_,
-                                  mget->state, si, fi, cs, cn);
+          matrix_mget_commit_var (rows, n_rows, d,
+                                  rowtype_, &mget->rowtypes,
+                                  mget->state,
+                                  ss, sn, si,
+                                  fs, fn, fi,
+                                  cs, cn,
+                                  pt, var_dimension);
           n_rows = 0;
           case_unref (cc);
           cc = case_ref (c);
@@ -6061,13 +7058,61 @@ matrix_cmd_execute_mget (struct mget_command *mget)
       if (n_rows >= allocated_rows)
         rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
       rows[n_rows++] = c;
+
+      if (change == SPLITS_CHANGED)
+        {
+          si++;
+          case_unref (sc);
+          sc = case_ref (c);
+
+          /* Reset the factor number, if there are factors. */
+          if (fn)
+            {
+              fi = 0;
+              if (row_has_factors)
+                fi++;
+              case_unref (fc);
+              fc = case_ref (c);
+            }
+        }
+      else if (change == FACTORS_CHANGED)
+        {
+          if (row_has_factors)
+            fi++;
+          case_unref (fc);
+          fc = case_ref (c);
+        }
     }
-  matrix_mget_commit_var (rows, n_rows, d, rowtype_,
-                          mget->state, si, fi, cs, cn);
+  matrix_mget_commit_var (rows, n_rows, d,
+                          rowtype_, &mget->rowtypes,
+                          mget->state,
+                          ss, sn, si,
+                          fs, fn, fi,
+                          cs, cn,
+                          pt, var_dimension);
   free (rows);
 
-exit:
-  casereader_destroy (r);
+  case_unref (sc);
+  case_unref (fc);
+  case_unref (cc);
+
+  if (var_dimension->n_leaves)
+    pivot_table_submit (pt);
+  else
+    pivot_table_unref (pt);
+}
+
+static void
+matrix_cmd_execute_mget (struct mget_command *mget)
+{
+  struct casereader *r;
+  struct dictionary *d;
+  if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
+                              mget->state->dataset, &r, &d))
+    {
+      matrix_cmd_execute_mget__ (mget, r, d);
+      matrix_close_casereader (mget->file, mget->state->dataset, r, d);
+    }
 }
 
 static bool
@@ -6209,8 +7254,8 @@ matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
         gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
     }
   else
-    msg (SE, _("SETDIAG argument 2 must be a scalar or a vector but it "
-               "has dimensions (%zu,%zu)."),
+    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);
 }
@@ -6270,6 +7315,8 @@ matrix_cmd_execute_svd (struct svd_command *svd)
     {
       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;
@@ -6354,6 +7401,107 @@ matrix_cmd_execute (struct matrix_cmd *cmd)
   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);
+      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);
+      free (cmd->write.format);
+      break;
+
+    case MCMD_GET:
+      matrix_lvalue_destroy (cmd->get.dst);
+      fh_unref (cmd->get.file);
+      free (cmd->get.encoding);
+      var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
+      break;
+
+    case MCMD_MSAVE:
+      matrix_expr_destroy (cmd->msave.expr);
+      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;
@@ -6425,6 +7573,7 @@ cmd_matrix (struct lexer *lexer, struct dataset *ds)
     return CMD_FAILURE;
 
   struct matrix_state state = {
+    .dataset = ds,
     .session = dataset_session (ds),
     .lexer = lexer,
     .vars = HMAP_INITIALIZER (state.vars),
@@ -6451,12 +7600,28 @@ cmd_matrix (struct lexer *lexer, struct dataset *ds)
         }
     }
 
-  if (state.common)
+  struct matrix_var *var, *next;
+  HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
     {
-      dict_unref (state.common->dict);
-      casewriter_destroy (state.common->writer);
-      free (state.common);
-    }
+      free (var->name);
+      gsl_matrix_free (var->value);
+      hmap_delete (&state.vars, &var->hmap_node);
+      free (var);
+    }
+  hmap_destroy (&state.vars);
+  msave_common_destroy (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);
+  fh_unref (state.prev_save_file);
+  for (size_t i = 0; i < state.n_save_files; i++)
+    save_file_destroy (state.save_files[i]);
+  free (state.save_files);
 
   return CMD_SUCCESS;
 }