MATRIX GET positive test
[pspp] / src / language / stats / matrix.c
index c0c80eda77f3a8005def95de95ef9b24d78c6a58..bcae19971750704affa257ef127ca17693ed9b53 100644 (file)
@@ -59,6 +59,7 @@
 #include "output/pivot-table.h"
 
 #include "gl/c-ctype.h"
+#include "gl/ftoastr.h"
 #include "gl/minmax.h"
 #include "gl/xsize.h"
 
@@ -69,7 +70,7 @@
 struct matrix_var
   {
     struct hmap_node hmap_node;
-    const char *name;
+    char *name;
     gsl_matrix *value;
   };
 
@@ -96,6 +97,14 @@ struct read_file
     char *encoding;
   };
 
+struct write_file
+  {
+    struct file_handle *file;
+    struct dfm_writer *writer;
+    char *encoding;
+    struct u8_line *held;
+  };
+
 struct matrix_state
   {
     struct dataset *dataset;
@@ -104,12 +113,15 @@ struct matrix_state
     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;
   };
 
 static struct matrix_var *
@@ -340,6 +352,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 +464,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;
 }
@@ -809,22 +822,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 +916,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 +951,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;
     }
 
@@ -1200,7 +1215,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;
@@ -1439,8 +1454,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;
@@ -1676,6 +1691,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)
@@ -2238,8 +2315,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 +2329,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 +2369,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 +2389,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 +2443,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++)
@@ -2458,7 +2544,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 +2556,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 +2593,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 +2629,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;
     }
@@ -2621,7 +2707,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 +2720,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;
     }
@@ -2933,7 +3019,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 +3059,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 +3078,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 +3125,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 +3140,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 +3237,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 +3245,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 +3275,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);
 }
 
 \f
@@ -3196,6 +3293,8 @@ struct matrix_cmds
     size_t n;
   };
 
+static void matrix_cmds_uninit (struct matrix_cmds *);
+
 struct matrix_cmd
   {
     enum matrix_cmd_type
@@ -3274,6 +3373,21 @@ struct matrix_cmd
           }
         loop;
 
+        struct display_command
+          {
+            struct matrix_state *state;
+            bool dictionary;
+            bool status;
+          }
+        display;
+
+        struct release_command
+          {
+            struct matrix_var **vars;
+            size_t n_vars;
+          }
+        release;
+
         struct save_command
           {
             struct matrix_expr *expression;
@@ -3292,7 +3406,7 @@ struct matrix_cmd
             int c1, c2;
             enum fmt_type format;
             int w;
-            int d;
+            //int d;
             bool symmetric;
             bool reread;
           }
@@ -3300,33 +3414,21 @@ 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;
@@ -3399,13 +3501,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);
-}
 \f
 static struct matrix_cmd *
 matrix_parse_compute (struct matrix_state *s)
@@ -3777,7 +3874,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 +3895,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
@@ -4535,6 +4635,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 +4681,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 +4727,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 +4771,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 +4856,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 +4872,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 +4883,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 +4892,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 +4902,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 +4911,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 +4936,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 +4958,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 +4970,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 +5005,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 +5027,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 +5045,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 +5057,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 +5076,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 +5131,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 +5147,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 +5182,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 +5222,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 +5236,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 +5272,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);
 }
@@ -5252,11 +5470,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 +5490,10 @@ matrix_parse_get (struct matrix_state *s)
           goto error;
         }
     }
+
+  if (get->user.treatment != MGET_ACCEPT)
+    get->system.treatment = MGET_ERROR;
+
   return cmd;
 
 error:
@@ -5336,6 +5558,22 @@ matrix_cmd_execute_get (struct get_command *get)
         }
     }
 
+  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;
@@ -5636,8 +5874,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;
@@ -6209,8 +6447,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 +6508,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 +6594,114 @@ 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);
+      fh_unref (cmd->save.outfile);
+      string_array_destroy (cmd->save.variables);
+      matrix_expr_destroy (cmd->save.names);
+      stringi_set_destroy (&cmd->save.strings);
+      break;
+
+    case MCMD_READ:
+      matrix_lvalue_destroy (cmd->read.dst);
+      matrix_expr_destroy (cmd->read.size);
+      break;
+
+    case MCMD_WRITE:
+      matrix_expr_destroy (cmd->write.expression);
+      free (cmd->write.format);
+      break;
+
+    case MCMD_GET:
+      matrix_lvalue_destroy (cmd->get.dst);
+      fh_unref (cmd->get.file);
+      free (cmd->get.encoding);
+      string_array_destroy (&cmd->get.variables);
+      break;
+
+    case MCMD_MSAVE:
+      free (cmd->msave.varname_);
+      matrix_expr_destroy (cmd->msave.expr);
+      matrix_expr_destroy (cmd->msave.factors);
+      matrix_expr_destroy (cmd->msave.splits);
+      break;
+
+    case MCMD_MGET:
+      fh_unref (cmd->mget.file);
+      stringi_set_destroy (&cmd->mget.rowtypes);
+      break;
+
+    case MCMD_EIGEN:
+      matrix_expr_destroy (cmd->eigen.expr);
+      break;
+
+    case MCMD_SETDIAG:
+      matrix_expr_destroy (cmd->setdiag.expr);
+      break;
+
+    case MCMD_SVD:
+      matrix_expr_destroy (cmd->svd.expr);
+      break;
+    }
+  free (cmd);
+}
+
 struct matrix_command_name
   {
     const char *name;
@@ -6451,12 +6799,30 @@ cmd_matrix (struct lexer *lexer, struct dataset *ds)
         }
     }
 
+  struct matrix_var *var, *next;
+  HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
+    {
+      free (var->name);
+      gsl_matrix_free (var->value);
+      hmap_delete (&state.vars, &var->hmap_node);
+      free (var);
+    }
+  hmap_destroy (&state.vars);
+  fh_unref (state.prev_save_outfile);
   if (state.common)
     {
       dict_unref (state.common->dict);
       casewriter_destroy (state.common->writer);
       free (state.common);
     }
+  fh_unref (state.prev_read_file);
+  for (size_t i = 0; i < state.n_read_files; i++)
+    read_file_destroy (state.read_files[i]);
+  free (state.read_files);
+  fh_unref (state.prev_write_file);
+  for (size_t i = 0; i < state.n_write_files; i++)
+    write_file_destroy (state.write_files[i]);
+  free (state.write_files);
 
   return CMD_SUCCESS;
 }