c59b3fd2db789e18b8fc66cc9da40039639e1a6c
[pspp] / src / language / stats / matrix.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2021 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_eigen.h>
21 #include <gsl/gsl_linalg.h>
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_permutation.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_vector.h>
26 #include <limits.h>
27 #include <math.h>
28 #include <uniwidth.h>
29
30 #include "data/any-reader.h"
31 #include "data/any-writer.h"
32 #include "data/casereader.h"
33 #include "data/casewriter.h"
34 #include "data/data-in.h"
35 #include "data/data-out.h"
36 #include "data/dataset.h"
37 #include "data/dictionary.h"
38 #include "data/file-handle-def.h"
39 #include "language/command.h"
40 #include "language/data-io/data-reader.h"
41 #include "language/data-io/data-writer.h"
42 #include "language/data-io/file-handle.h"
43 #include "language/lexer/format-parser.h"
44 #include "language/lexer/lexer.h"
45 #include "language/lexer/variable-parser.h"
46 #include "libpspp/array.h"
47 #include "libpspp/assertion.h"
48 #include "libpspp/compiler.h"
49 #include "libpspp/hmap.h"
50 #include "libpspp/i18n.h"
51 #include "libpspp/misc.h"
52 #include "libpspp/str.h"
53 #include "libpspp/string-array.h"
54 #include "libpspp/stringi-set.h"
55 #include "libpspp/u8-line.h"
56 #include "math/random.h"
57 #include "output/driver.h"
58 #include "output/output-item.h"
59 #include "output/pivot-table.h"
60
61 #include "gl/c-ctype.h"
62 #include "gl/minmax.h"
63 #include "gl/xsize.h"
64
65 #include "gettext.h"
66 #define _(msgid) gettext (msgid)
67 #define N_(msgid) (msgid)
68
69 struct matrix_var
70   {
71     struct hmap_node hmap_node;
72     char *name;
73     gsl_matrix *value;
74   };
75
76 struct msave_common
77   {
78     /* Configuration. */
79     struct file_handle *outfile;
80     struct string_array variables;
81     struct string_array fnames;
82     struct string_array snames;
83     bool has_factors;
84     bool has_splits;
85     size_t n_varnames;
86
87     /* Execution state. */
88     struct dictionary *dict;
89     struct casewriter *writer;
90   };
91
92 struct read_file
93   {
94     struct file_handle *file;
95     struct dfm_reader *reader;
96     char *encoding;
97   };
98
99 struct matrix_state
100   {
101     struct dataset *dataset;
102     struct session *session;
103     struct lexer *lexer;
104     struct hmap vars;
105     bool in_loop;
106     struct file_handle *prev_save_outfile;
107     struct file_handle *prev_write_outfile;
108     struct msave_common *common;
109
110     struct file_handle *prev_read_file;
111     struct read_file **read_files;
112     size_t n_read_files;
113   };
114
115 static struct matrix_var *
116 matrix_var_lookup (struct matrix_state *s, struct substring name)
117 {
118   struct matrix_var *var;
119
120   HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
121                            utf8_hash_case_substring (name, 0), &s->vars)
122     if (!utf8_sscasecmp (ss_cstr (var->name), name))
123       return var;
124
125   return NULL;
126 }
127
128 static struct matrix_var *
129 matrix_var_create (struct matrix_state *s, struct substring name)
130 {
131   struct matrix_var *var = xmalloc (sizeof *var);
132   *var = (struct matrix_var) { .name = ss_xstrdup (name) };
133   hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
134   return var;
135 }
136
137 static void
138 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
139 {
140   gsl_matrix_free (var->value);
141   var->value = value;
142 }
143 \f
144 #define MATRIX_FUNCTIONS \
145     F(ABS, m_m) \
146     F(ALL, d_m) \
147     F(ANY, d_m) \
148     F(ARSIN, m_m) \
149     F(ARTAN, m_m) \
150     F(BLOCK, m_any) \
151     F(CHOL, m_m) \
152     F(CMIN, m_m) \
153     F(CMAX, m_m) \
154     F(COS, m_m) \
155     F(CSSQ, m_m) \
156     F(CSUM, m_m) \
157     F(DESIGN, m_m) \
158     F(DET, d_m) \
159     F(DIAG, m_m) \
160     F(EVAL, m_m) \
161     F(EXP, m_m) \
162     F(GINV, m_m) \
163     F(GRADE, m_m) \
164     F(GSCH, m_m) \
165     F(IDENT, IDENT) \
166     F(INV, m_m) \
167     F(KRONEKER, m_mm) \
168     F(LG10, m_m) \
169     F(LN, m_m) \
170     F(MAGIC, m_d) \
171     F(MAKE, m_ddd) \
172     F(MDIAG, m_v) \
173     F(MMAX, d_m) \
174     F(MMIN, d_m) \
175     F(MOD, m_md) \
176     F(MSSQ, d_m) \
177     F(MSUM, d_m) \
178     F(NCOL, d_m) \
179     F(NROW, d_m) \
180     F(RANK, d_m) \
181     F(RESHAPE, m_mdd) \
182     F(RMAX, m_m) \
183     F(RMIN, m_m) \
184     F(RND, m_m) \
185     F(RNKORDER, m_m) \
186     F(RSSQ, m_m) \
187     F(RSUM, m_m) \
188     F(SIN, m_m) \
189     F(SOLVE, m_mm) \
190     F(SQRT, m_m) \
191     F(SSCP, m_m) \
192     F(SVAL, m_m) \
193     F(SWEEP, m_md) \
194     F(T, m_m) \
195     F(TRACE, d_m) \
196     F(TRANSPOS, m_m) \
197     F(TRUNC, m_m) \
198     F(UNIFORM, m_dd)
199
200 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
201 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
202 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
203 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
204 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
205 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
206 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
207 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
208 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
209 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
210 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
211
212 typedef gsl_matrix *matrix_proto_m_d (double);
213 typedef gsl_matrix *matrix_proto_m_dd (double, double);
214 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
215 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
216 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
217 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
218 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
219 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
220 typedef double matrix_proto_d_m (gsl_matrix *);
221 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
222 typedef gsl_matrix *matrix_proto_IDENT (double, double);
223
224 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
225 MATRIX_FUNCTIONS
226 #undef F
227 \f
228 /* Matrix expressions. */
229
230 struct matrix_expr
231   {
232     enum matrix_op
233       {
234         /* Functions. */
235 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
236         MATRIX_FUNCTIONS
237 #undef F
238
239         /* Elementwise and scalar arithmetic. */
240         MOP_NEGATE,             /* unary - */
241         MOP_ADD_ELEMS,          /* + */
242         MOP_SUB_ELEMS,          /* - */
243         MOP_MUL_ELEMS,          /* &* */
244         MOP_DIV_ELEMS,          /* / and &/ */
245         MOP_EXP_ELEMS,          /* &** */
246         MOP_SEQ,                /* a:b */
247         MOP_SEQ_BY,             /* a:b:c */
248
249         /* Matrix arithmetic. */
250         MOP_MUL_MAT,            /* * */
251         MOP_EXP_MAT,            /* ** */
252
253         /* Relational. */
254         MOP_GT,                 /* > */
255         MOP_GE,                 /* >= */
256         MOP_LT,                 /* < */
257         MOP_LE,                 /* <= */
258         MOP_EQ,                 /* = */
259         MOP_NE,                 /* <> */
260
261         /* Logical. */
262         MOP_NOT,                /* NOT */
263         MOP_AND,                /* AND */
264         MOP_OR,                 /* OR */
265         MOP_XOR,                /* XOR */
266
267         /* {}. */
268         MOP_PASTE_HORZ,         /* a, b, c, ... */
269         MOP_PASTE_VERT,         /* a; b; c; ... */
270         MOP_EMPTY,              /* {} */
271
272         /* Sub-matrices. */
273         MOP_VEC_INDEX,          /* x(y) */
274         MOP_VEC_ALL,            /* x(:) */
275         MOP_MAT_INDEX,          /* x(y,z) */
276         MOP_ROW_INDEX,          /* x(y,:) */
277         MOP_COL_INDEX,          /* x(:,z) */
278
279         /* Literals. */
280         MOP_NUMBER,
281         MOP_VARIABLE,
282
283         /* Oddball stuff. */
284         MOP_EOF,                /* EOF('file') */
285       }
286     op;
287
288     union
289       {
290         struct
291           {
292             struct matrix_expr **subs;
293             size_t n_subs;
294           };
295
296         double number;
297         struct matrix_var *variable;
298         struct read_file *eof;
299       };
300   };
301
302 static void
303 matrix_expr_destroy (struct matrix_expr *e)
304 {
305   if (!e)
306     return;
307
308   switch (e->op)
309     {
310 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
311 MATRIX_FUNCTIONS
312 #undef F
313     case MOP_NEGATE:
314     case MOP_ADD_ELEMS:
315     case MOP_SUB_ELEMS:
316     case MOP_MUL_ELEMS:
317     case MOP_DIV_ELEMS:
318     case MOP_EXP_ELEMS:
319     case MOP_SEQ:
320     case MOP_SEQ_BY:
321     case MOP_MUL_MAT:
322     case MOP_EXP_MAT:
323     case MOP_GT:
324     case MOP_GE:
325     case MOP_LT:
326     case MOP_LE:
327     case MOP_EQ:
328     case MOP_NE:
329     case MOP_NOT:
330     case MOP_AND:
331     case MOP_OR:
332     case MOP_XOR:
333     case MOP_EMPTY:
334     case MOP_PASTE_HORZ:
335     case MOP_PASTE_VERT:
336     case MOP_VEC_INDEX:
337     case MOP_VEC_ALL:
338     case MOP_MAT_INDEX:
339     case MOP_ROW_INDEX:
340     case MOP_COL_INDEX:
341       for (size_t i = 0; i < e->n_subs; i++)
342         matrix_expr_destroy (e->subs[i]);
343       free (e->subs);
344       break;
345
346     case MOP_NUMBER:
347     case MOP_VARIABLE:
348     case MOP_EOF:
349       break;
350     }
351   free (e);
352 }
353
354 static struct matrix_expr *
355 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
356                          size_t n_subs)
357 {
358   struct matrix_expr *e = xmalloc (sizeof *e);
359   *e = (struct matrix_expr) {
360     .op = op,
361     .subs = xmemdup (subs, n_subs * sizeof *subs),
362     .n_subs = n_subs
363   };
364   return e;
365 }
366
367 static struct matrix_expr *
368 matrix_expr_create_0 (enum matrix_op op)
369 {
370   struct matrix_expr *sub;
371   return matrix_expr_create_subs (op, &sub, 0);
372 }
373
374 static struct matrix_expr *
375 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
376 {
377   return matrix_expr_create_subs (op, &sub, 1);
378 }
379
380 static struct matrix_expr *
381 matrix_expr_create_2 (enum matrix_op op,
382                       struct matrix_expr *sub0, struct matrix_expr *sub1)
383 {
384   struct matrix_expr *subs[] = { sub0, sub1 };
385   return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
386 }
387
388 static struct matrix_expr *
389 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
390                       struct matrix_expr *sub1, struct matrix_expr *sub2)
391 {
392   struct matrix_expr *subs[] = { sub0, sub1, sub2 };
393   return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
394 }
395
396 static struct matrix_expr *
397 matrix_expr_create_number (double number)
398 {
399   struct matrix_expr *e = xmalloc (sizeof *e);
400   *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
401   return e;
402 }
403
404 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
405
406 static struct matrix_expr *
407 matrix_parse_curly_comma (struct matrix_state *s)
408 {
409   struct matrix_expr *lhs = matrix_parse_or_xor (s);
410   if (!lhs)
411     return NULL;
412
413   while (lex_match (s->lexer, T_COMMA))
414     {
415       struct matrix_expr *rhs = matrix_parse_or_xor (s);
416       if (!rhs)
417         {
418           matrix_expr_destroy (lhs);
419           return NULL;
420         }
421       lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
422     }
423   return lhs;
424 }
425
426 static struct matrix_expr *
427 matrix_parse_curly_semi (struct matrix_state *s)
428 {
429   if (lex_token (s->lexer) == T_RCURLY)
430     return matrix_expr_create_0 (MOP_EMPTY);
431
432   struct matrix_expr *lhs = matrix_parse_curly_comma (s);
433   if (!lhs)
434     return NULL;
435
436   while (lex_match (s->lexer, T_SEMICOLON))
437     {
438       struct matrix_expr *rhs = matrix_parse_curly_comma (s);
439       if (!rhs)
440         {
441           matrix_expr_destroy (lhs);
442           return NULL;
443         }
444       lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
445     }
446   return lhs;
447 }
448
449 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M)                     \
450   for (size_t Y = 0; Y < (M)->size1; Y++)                       \
451     for (size_t X = 0; X < (M)->size2; X++)                     \
452       for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
453
454 static bool
455 is_vector (const gsl_matrix *m)
456 {
457   return m->size1 <= 1 || m->size2 <= 1;
458 }
459
460 static gsl_vector
461 to_vector (gsl_matrix *m)
462 {
463   return (m->size1 == 1
464           ? gsl_matrix_row (m, 0).vector
465           : gsl_matrix_column (m, 0).vector);
466 }
467
468
469 static gsl_matrix *
470 matrix_eval_ABS (gsl_matrix *m)
471 {
472   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
473     *d = fabs (*d);
474   return m;
475 }
476
477 static double
478 matrix_eval_ALL (gsl_matrix *m)
479 {
480   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
481     if (*d == 0.0)
482       return 0.0;
483   return 1.0;
484 }
485
486 static double
487 matrix_eval_ANY (gsl_matrix *m)
488 {
489   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
490     if (*d != 0.0)
491       return !.0;
492   return 0.0;
493 }
494
495 static gsl_matrix *
496 matrix_eval_ARSIN (gsl_matrix *m)
497 {
498   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
499     *d = asin (*d);
500   return m;
501 }
502
503 static gsl_matrix *
504 matrix_eval_ARTAN (gsl_matrix *m)
505 {
506   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
507     *d = atan (*d);
508   return m;
509 }
510
511 static gsl_matrix *
512 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
513 {
514   size_t r = 0;
515   size_t c = 0;
516   for (size_t i = 0; i < n; i++)
517     {
518       r += m[i]->size1;
519       c += m[i]->size2;
520     }
521   gsl_matrix *block = gsl_matrix_calloc (r, c);
522   r = c = 0;
523   for (size_t i = 0; i < n; i++)
524     {
525       for (size_t y = 0; y < m[i]->size1; y++)
526         for (size_t x = 0; x < m[i]->size2; x++)
527           gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
528       r += m[i]->size1;
529       c += m[i]->size2;
530     }
531   return block;
532 }
533
534 static gsl_matrix *
535 matrix_eval_CHOL (gsl_matrix *m)
536 {
537   if (!gsl_linalg_cholesky_decomp1 (m))
538     {
539       for (size_t y = 0; y < m->size1; y++)
540         for (size_t x = y + 1; x < m->size2; x++)
541           gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
542
543       for (size_t y = 0; y < m->size1; y++)
544         for (size_t x = 0; x < y; x++)
545           gsl_matrix_set (m, y, x, 0);
546       return m;
547     }
548   else
549     {
550       msg (SE, _("Input to CHOL function is not positive-definite."));
551       return NULL;
552     }
553 }
554
555 static gsl_matrix *
556 matrix_eval_col_extremum (gsl_matrix *m, bool min)
557 {
558   if (m->size1 <= 1)
559     return m;
560   else if (!m->size2)
561     return gsl_matrix_alloc (1, 0);
562
563   gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
564   for (size_t x = 0; x < m->size2; x++)
565     {
566       double ext = gsl_matrix_get (m, 0, x);
567       for (size_t y = 1; y < m->size1; y++)
568         {
569           double value = gsl_matrix_get (m, y, x);
570           if (min ? value < ext : value > ext)
571             ext = value;
572         }
573       gsl_matrix_set (cext, 0, x, ext);
574     }
575   return cext;
576 }
577
578 static gsl_matrix *
579 matrix_eval_CMAX (gsl_matrix *m)
580 {
581   return matrix_eval_col_extremum (m, false);
582 }
583
584 static gsl_matrix *
585 matrix_eval_CMIN (gsl_matrix *m)
586 {
587   return matrix_eval_col_extremum (m, true);
588 }
589
590 static gsl_matrix *
591 matrix_eval_COS (gsl_matrix *m)
592 {
593   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
594     *d = cos (*d);
595   return m;
596 }
597
598 static gsl_matrix *
599 matrix_eval_col_sum (gsl_matrix *m, bool square)
600 {
601   if (m->size1 == 0)
602     return m;
603   else if (!m->size2)
604     return gsl_matrix_alloc (1, 0);
605
606   gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
607   for (size_t x = 0; x < m->size2; x++)
608     {
609       double sum = 0;
610       for (size_t y = 0; y < m->size1; y++)
611         {
612           double d = gsl_matrix_get (m, y, x);
613           sum += square ? pow2 (d) : d;
614         }
615       gsl_matrix_set (result, 0, x, sum);
616     }
617   return result;
618 }
619
620 static gsl_matrix *
621 matrix_eval_CSSQ (gsl_matrix *m)
622 {
623   return matrix_eval_col_sum (m, true);
624 }
625
626 static gsl_matrix *
627 matrix_eval_CSUM (gsl_matrix *m)
628 {
629   return matrix_eval_col_sum (m, false);
630 }
631
632 static int
633 compare_double_3way (const void *a_, const void *b_)
634 {
635   const double *a = a_;
636   const double *b = b_;
637   return *a < *b ? -1 : *a > *b;
638 }
639
640 static gsl_matrix *
641 matrix_eval_DESIGN (gsl_matrix *m)
642 {
643   double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
644   gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
645   gsl_matrix_transpose_memcpy (&m2, m);
646
647   for (size_t y = 0; y < m2.size1; y++)
648     qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
649
650   size_t *n = xcalloc (m2.size1, sizeof *n);
651   size_t n_total = 0;
652   for (size_t i = 0; i < m2.size1; i++)
653     {
654       double *row = tmp + m2.size2 * i;
655       for (size_t j = 0; j < m2.size2; )
656         {
657           size_t k;
658           for (k = j + 1; k < m2.size2; k++)
659             if (row[j] != row[k])
660               break;
661           row[n[i]++] = row[j];
662           j = k;
663         }
664
665       if (n[i] <= 1)
666         msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
667       else
668         n_total += n[i];
669     }
670
671   gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
672   size_t x = 0;
673   for (size_t i = 0; i < m->size2; i++)
674     {
675       if (n[i] <= 1)
676         continue;
677
678       const double *unique = tmp + m2.size2 * i;
679       for (size_t j = 0; j < n[i]; j++, x++)
680         {
681           double value = unique[j];
682           for (size_t y = 0; y < m->size1; y++)
683             gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
684         }
685     }
686
687   free (n);
688   free (tmp);
689
690   return result;
691 }
692
693 static double
694 matrix_eval_DET (gsl_matrix *m)
695 {
696   gsl_permutation *p = gsl_permutation_alloc (m->size1);
697   int signum;
698   gsl_linalg_LU_decomp (m, p, &signum);
699   gsl_permutation_free (p);
700   return gsl_linalg_LU_det (m, signum);
701 }
702
703 static gsl_matrix *
704 matrix_eval_DIAG (gsl_matrix *m)
705 {
706   gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
707   for (size_t i = 0; i < diag->size1; i++)
708     gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
709   return diag;
710 }
711
712 static bool
713 is_symmetric (const gsl_matrix *m)
714 {
715   if (m->size1 != m->size2)
716     return false;
717
718   for (size_t y = 0; y < m->size1; y++)
719     for (size_t x = 0; x < y; x++)
720       if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
721         return false;
722
723   return true;
724 }
725
726 static int
727 compare_double_desc (const void *a_, const void *b_)
728 {
729   const double *a = a_;
730   const double *b = b_;
731   return *a > *b ? -1 : *a < *b;
732 }
733
734 static gsl_matrix *
735 matrix_eval_EVAL (gsl_matrix *m)
736 {
737   if (!is_symmetric (m))
738     {
739       msg (SE, _("Argument of EVAL must be symmetric."));
740       return NULL;
741     }
742
743   gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
744   gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
745   gsl_vector v_eval = to_vector (eval);
746   gsl_eigen_symm (m, &v_eval, w);
747   gsl_eigen_symm_free (w);
748
749   assert (v_eval.stride == 1);
750   qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
751
752   return eval;
753 }
754
755 static gsl_matrix *
756 matrix_eval_EXP (gsl_matrix *m)
757 {
758   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
759     *d = exp (*d);
760   return m;
761 }
762
763 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
764    marked as:
765
766    Charl Linssen <charl@itfromb.it>
767    Feb 2016
768    PUBLIC DOMAIN */
769 static gsl_matrix *
770 matrix_eval_GINV (gsl_matrix *A)
771 {
772   size_t n = A->size1;
773   size_t m = A->size2;
774   bool swap = m > n;
775   gsl_matrix *tmp_mat = NULL;
776   if (swap)
777     {
778       /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
779       tmp_mat = gsl_matrix_alloc (m, n);
780       gsl_matrix_transpose_memcpy (tmp_mat, A);
781       A = tmp_mat;
782       size_t i = m;
783       m = n;
784       n = i;
785     }
786
787   /* Do SVD. */
788   gsl_matrix *V = gsl_matrix_alloc (m, m);
789   gsl_vector *u = gsl_vector_alloc (m);
790
791   gsl_vector *tmp_vec = gsl_vector_alloc (m);
792   gsl_linalg_SV_decomp (A, V, u, tmp_vec);
793   gsl_vector_free (tmp_vec);
794
795   /* Compute Σ⁻¹. */
796   gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
797   gsl_matrix_set_zero (Sigma_pinv);
798   double cutoff = 1e-15 * gsl_vector_max (u);
799
800   for (size_t i = 0; i < m; ++i)
801     {
802       double x = gsl_vector_get (u, i);
803       gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
804     }
805
806   /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
807   gsl_matrix *U = gsl_matrix_calloc (n, n);
808   for (size_t i = 0; i < n; i++)
809     for (size_t j = 0; j < m; j++)
810       gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
811
812   /* two dot products to obtain pseudoinverse */
813   gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
814   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
815
816   gsl_matrix *A_pinv;
817   if (swap)
818     {
819       A_pinv = gsl_matrix_alloc (n, m);
820       gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
821     }
822   else
823     {
824       A_pinv = gsl_matrix_alloc (m, n);
825       gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
826     }
827
828   gsl_matrix_free (tmp_mat);
829   gsl_matrix_free (tmp_mat2);
830   gsl_matrix_free (U);
831   gsl_matrix_free (Sigma_pinv);
832   gsl_vector_free (u);
833   gsl_matrix_free (V);
834
835   return A_pinv;
836 }
837
838 struct grade
839   {
840     size_t y, x;
841     double value;
842   };
843
844 static int
845 grade_compare_3way (const void *a_, const void *b_)
846 {
847   const struct grade *a = a_;
848   const struct grade *b = b_;
849
850   return (a->value < b->value ? -1
851           : a->value > b->value ? 1
852           : a->y < b->y ? -1
853           : a->y > b->y ? 1
854           : a->x < b->x ? -1
855           : a->x > b->x);
856 }
857
858 static gsl_matrix *
859 matrix_eval_GRADE (gsl_matrix *m)
860 {
861   size_t n = m->size1 * m->size2;
862   struct grade *grades = xmalloc (n * sizeof *grades);
863
864   size_t i = 0;
865   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
866     grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
867   qsort (grades, n, sizeof *grades, grade_compare_3way);
868
869   for (size_t i = 0; i < n; i++)
870     gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
871
872   free (grades);
873
874   return m;
875 }
876
877 static double
878 dot (gsl_vector *a, gsl_vector *b)
879 {
880   double result = 0.0;
881   for (size_t i = 0; i < a->size; i++)
882     result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
883   return result;
884 }
885
886 static double
887 norm2 (gsl_vector *v)
888 {
889   double result = 0.0;
890   for (size_t i = 0; i < v->size; i++)
891     result += pow2 (gsl_vector_get (v, i));
892   return result;
893 }
894
895 static double
896 norm (gsl_vector *v)
897 {
898   return sqrt (norm2 (v));
899 }
900
901 static gsl_matrix *
902 matrix_eval_GSCH (gsl_matrix *v)
903 {
904   if (v->size2 < v->size1)
905     {
906       msg (SE, _("GSCH requires its argument to have at least as many columns "
907                  "as rows, but it has dimensions %zu×%zu."),
908            v->size1, v->size2);
909       return NULL;
910     }
911   if (!v->size1 || !v->size2)
912     return v;
913
914   gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
915   size_t ux = 0;
916   for (size_t vx = 0; vx < v->size2; vx++)
917     {
918       gsl_vector u_i = gsl_matrix_column (u, ux).vector;
919       gsl_vector v_i = gsl_matrix_column (v, vx).vector;
920
921       gsl_vector_memcpy (&u_i, &v_i);
922       for (size_t j = 0; j < ux; j++)
923         {
924           gsl_vector u_j = gsl_matrix_column (u, j).vector;
925           double scale = dot (&u_j, &u_i) / norm2 (&u_j);
926           for (size_t k = 0; k < u_i.size; k++)
927             gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
928                                       - scale * gsl_vector_get (&u_j, k)));
929         }
930
931       double len = norm (&u_i);
932       if (len > 1e-15)
933         {
934           gsl_vector_scale (&u_i, 1.0 / len);
935           if (++ux >= v->size1)
936             break;
937         }
938     }
939
940   if (ux < v->size1)
941     {
942       msg (SE, _("%zu×%zu argument to GSCH contains only "
943                  "%zu linearly independent columns."),
944            v->size1, v->size2, ux);
945       gsl_matrix_free (u);
946       return NULL;
947     }
948
949   u->size2 = v->size1;
950   return u;
951 }
952
953 static gsl_matrix *
954 matrix_eval_IDENT (double s1, double s2)
955 {
956   if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
957     {
958       msg (SE, _("Arguments to IDENT must be non-negative integers."));
959       return NULL;
960     }
961
962   gsl_matrix *m = gsl_matrix_alloc (s1, s2);
963   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
964     *d = x == y;
965   return m;
966 }
967
968 static void
969 invert_matrix (gsl_matrix *x)
970 {
971   gsl_permutation *p = gsl_permutation_alloc (x->size1);
972   int signum;
973   gsl_linalg_LU_decomp (x, p, &signum);
974   gsl_linalg_LU_invx (x, p);
975   gsl_permutation_free (p);
976 }
977
978 static gsl_matrix *
979 matrix_eval_INV (gsl_matrix *m)
980 {
981   invert_matrix (m);
982   return m;
983 }
984
985 static gsl_matrix *
986 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
987 {
988   gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
989                                     a->size2 * b->size2);
990   size_t y = 0;
991   for (size_t ar = 0; ar < a->size1; ar++)
992     for (size_t br = 0; br < b->size1; br++, y++)
993       {
994         size_t x = 0;
995         for (size_t ac = 0; ac < a->size2; ac++)
996           for (size_t bc = 0; bc < b->size2; bc++, x++)
997             {
998               double av = gsl_matrix_get (a, ar, ac);
999               double bv = gsl_matrix_get (b, br, bc);
1000               gsl_matrix_set (k, y, x, av * bv);
1001             }
1002       }
1003   return k;
1004 }
1005
1006 static gsl_matrix *
1007 matrix_eval_LG10 (gsl_matrix *m)
1008 {
1009   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1010     *d = log10 (*d);
1011   return m;
1012 }
1013
1014 static gsl_matrix *
1015 matrix_eval_LN (gsl_matrix *m)
1016 {
1017   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1018     *d = log (*d);
1019   return m;
1020 }
1021
1022 static void
1023 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1024 {
1025   /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1026   size_t y = 0;
1027   size_t x = n / 2;
1028   for (size_t i = 1; i <= n * n; i++)
1029     {
1030       gsl_matrix_set (m, y, x, i);
1031
1032       size_t y1 = !y ? n - 1 : y - 1;
1033       size_t x1 = x + 1 >= n ? 0 : x + 1;
1034       if (gsl_matrix_get (m, y1, x1) == 0)
1035         {
1036           y = y1;
1037           x = x1;
1038         }
1039       else
1040         y = y + 1 >= n ? 0 : y + 1;
1041     }
1042 }
1043
1044 static void
1045 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1046 {
1047   double a = gsl_matrix_get (m, y1, x1);
1048   double b = gsl_matrix_get (m, y2, x2);
1049   gsl_matrix_set (m, y1, x1, b);
1050   gsl_matrix_set (m, y2, x2, a);
1051 }
1052
1053 static void
1054 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1055 {
1056   size_t x, y;
1057
1058   /* A. Umar, "On the Construction of Even Order Magic Squares",
1059      https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1060   x = y = 0;
1061   for (size_t i = 1; i <= n * n / 2; i++)
1062     {
1063       gsl_matrix_set (m, y, x, i);
1064       if (++y >= n)
1065         {
1066           y = 0;
1067           x++;
1068         }
1069     }
1070
1071   x = n - 1;
1072   y = 0;
1073   for (size_t i = n * n; i > n * n / 2; i--)
1074     {
1075       gsl_matrix_set (m, y, x, i);
1076       if (++y >= n)
1077         {
1078           y = 0;
1079           x--;
1080         }
1081     }
1082
1083   for (size_t y = 0; y < n; y++)
1084     for (size_t x = 0; x < n / 2; x++)
1085       {
1086         unsigned int d = gsl_matrix_get (m, y, x);
1087         if (d % 2 != (y < n / 2))
1088           magic_exchange (m, y, x, y, n - x - 1);
1089       }
1090
1091   size_t y1 = n / 2;
1092   size_t y2 = n - 1;
1093   size_t x1 = n / 2 - 1;
1094   size_t x2 = n / 2;
1095   magic_exchange (m, y1, x1, y2, x1);
1096   magic_exchange (m, y1, x2, y2, x2);
1097 }
1098
1099 static void
1100 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1101 {
1102   /* A. Umar, "On the Construction of Even Order Magic Squares",
1103      https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1104   size_t x, y;
1105
1106   x = y = 0;
1107   for (size_t i = 1; ; i++)
1108     {
1109       gsl_matrix_set (m, y, x, i);
1110       if (++y == n / 2 - 1)
1111         y += 2;
1112       else if (y >= n)
1113         {
1114           y = 0;
1115           if (++x >= n / 2)
1116             break;
1117         }
1118     }
1119
1120   x = n - 1;
1121   y = 0;
1122   for (size_t i = n * n; ; i--)
1123     {
1124       gsl_matrix_set (m, y, x, i);
1125       if (++y == n / 2 - 1)
1126         y += 2;
1127       else if (y >= n)
1128         {
1129           y = 0;
1130           if (--x < n / 2)
1131             break;
1132         }
1133     }
1134   for (size_t y = 0; y < n; y++)
1135     if (y != n / 2 - 1 && y != n / 2)
1136       for (size_t x = 0; x < n / 2; x++)
1137         {
1138           unsigned int d = gsl_matrix_get (m, y, x);
1139           if (d % 2 != (y < n / 2))
1140             magic_exchange (m, y, x, y, n - x - 1);
1141         }
1142
1143   size_t a0 = (n * n - 2 * n) / 2 + 1;
1144   for (size_t i = 0; i < n / 2; i++)
1145     {
1146       size_t a = a0 + i;
1147       gsl_matrix_set (m, n / 2, i, a);
1148       gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1149     }
1150   for (size_t i = 0; i < n / 2; i++)
1151     {
1152       size_t a = a0 + i + n / 2;
1153       gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1154       gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1155     }
1156   for (size_t x = 1; x < n / 2; x += 2)
1157     magic_exchange (m, n / 2, x, n / 2 - 1, x);
1158   for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1159     magic_exchange (m, n / 2, x, n / 2 - 1, x);
1160   size_t x1 = n / 2 - 2;
1161   size_t x2 = n / 2 + 1;
1162   size_t y1 = n / 2 - 2;
1163   size_t y2 = n / 2 + 1;
1164   magic_exchange (m, y1, x1, y2, x1);
1165   magic_exchange (m, y1, x2, y2, x2);
1166 }
1167
1168 static gsl_matrix *
1169 matrix_eval_MAGIC (double n_)
1170 {
1171   if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1172     {
1173       msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1174       return NULL;
1175     }
1176   size_t n = n_;
1177
1178   gsl_matrix *m = gsl_matrix_calloc (n, n);
1179   if (n % 2)
1180     matrix_eval_MAGIC_odd (m, n);
1181   else if (n % 4)
1182     matrix_eval_MAGIC_singly_even (m, n);
1183   else
1184     matrix_eval_MAGIC_doubly_even (m, n);
1185   return m;
1186 }
1187
1188 static gsl_matrix *
1189 matrix_eval_MAKE (double r, double c, double value)
1190 {
1191   if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1192     {
1193       msg (SE, _("Size arguments to MAKE must be integers."));
1194       return NULL;
1195     }
1196
1197   gsl_matrix *m = gsl_matrix_alloc (r, c);
1198   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1199     *d = value;
1200   return m;
1201 }
1202
1203 static gsl_matrix *
1204 matrix_eval_MDIAG (gsl_vector *v)
1205 {
1206   gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1207   gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1208   gsl_vector_memcpy (&diagonal, v);
1209   return m;
1210 }
1211
1212 static double
1213 matrix_eval_MMAX (gsl_matrix *m)
1214 {
1215   return gsl_matrix_max (m);
1216 }
1217
1218 static double
1219 matrix_eval_MMIN (gsl_matrix *m)
1220 {
1221   return gsl_matrix_min (m);
1222 }
1223
1224 static gsl_matrix *
1225 matrix_eval_MOD (gsl_matrix *m, double divisor)
1226 {
1227   if (divisor == 0.0)
1228     {
1229       msg (SE, _("Divisor argument to MOD function must be nonzero."));
1230       return NULL;
1231     }
1232
1233   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1234     *d = fmod (*d, divisor);
1235   return m;
1236 }
1237
1238 static double
1239 matrix_eval_MSSQ (gsl_matrix *m)
1240 {
1241   double mssq = 0.0;
1242   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1243     mssq += *d * *d;
1244   return mssq;
1245 }
1246
1247 static double
1248 matrix_eval_MSUM (gsl_matrix *m)
1249 {
1250   double msum = 0.0;
1251   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1252     msum += *d;
1253   return msum;
1254 }
1255
1256 static double
1257 matrix_eval_NCOL (gsl_matrix *m)
1258 {
1259   return m->size2;
1260 }
1261
1262 static double
1263 matrix_eval_NROW (gsl_matrix *m)
1264 {
1265   return m->size1;
1266 }
1267
1268 static double
1269 matrix_eval_RANK (gsl_matrix *m)
1270 {
1271   gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1272   gsl_linalg_QR_decomp (m, tau);
1273   gsl_vector_free (tau);
1274
1275   return gsl_linalg_QRPT_rank (m, -1);
1276 }
1277
1278 static gsl_matrix *
1279 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1280 {
1281   if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1282     {
1283       msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1284       return NULL;
1285     }
1286   size_t r = r_;
1287   size_t c = c_;
1288   if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1289     {
1290       msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1291                  "product of matrix dimensions."));
1292       return NULL;
1293     }
1294
1295   gsl_matrix *dst = gsl_matrix_alloc (r, c);
1296   size_t y1 = 0;
1297   size_t x1 = 0;
1298   MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1299     {
1300       gsl_matrix_set (dst, y1, x1, *d);
1301       if (++x1 >= c)
1302         {
1303           x1 = 0;
1304           y1++;
1305         }
1306     }
1307   return dst;
1308 }
1309
1310 static gsl_matrix *
1311 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1312 {
1313   if (m->size2 <= 1)
1314     return m;
1315   else if (!m->size1)
1316     return gsl_matrix_alloc (0, 1);
1317
1318   gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1319   for (size_t y = 0; y < m->size1; y++)
1320     {
1321       double ext = gsl_matrix_get (m, y, 0);
1322       for (size_t x = 1; x < m->size2; x++)
1323         {
1324           double value = gsl_matrix_get (m, y, x);
1325           if (min ? value < ext : value > ext)
1326             ext = value;
1327         }
1328       gsl_matrix_set (rext, y, 0, ext);
1329     }
1330   return rext;
1331 }
1332
1333 static gsl_matrix *
1334 matrix_eval_RMAX (gsl_matrix *m)
1335 {
1336   return matrix_eval_row_extremum (m, false);
1337 }
1338
1339 static gsl_matrix *
1340 matrix_eval_RMIN (gsl_matrix *m)
1341 {
1342   return matrix_eval_row_extremum (m, true);
1343 }
1344
1345 static gsl_matrix *
1346 matrix_eval_RND (gsl_matrix *m)
1347 {
1348   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1349     *d = rint (*d);
1350   return m;
1351 }
1352
1353 struct rank
1354   {
1355     size_t y, x;
1356     double value;
1357   };
1358
1359 static int
1360 rank_compare_3way (const void *a_, const void *b_)
1361 {
1362   const struct rank *a = a_;
1363   const struct rank *b = b_;
1364
1365   return a->value < b->value ? -1 : a->value > b->value;
1366 }
1367
1368 static gsl_matrix *
1369 matrix_eval_RNKORDER (gsl_matrix *m)
1370 {
1371   size_t n = m->size1 * m->size2;
1372   struct rank *ranks = xmalloc (n * sizeof *ranks);
1373   size_t i = 0;
1374   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1375     ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1376   qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1377
1378   for (size_t i = 0; i < n; )
1379     {
1380       size_t j;
1381       for (j = i + 1; j < n; j++)
1382         if (ranks[i].value != ranks[j].value)
1383           break;
1384
1385       double rank = (i + j + 1.0) / 2.0;
1386       for (size_t k = i; k < j; k++)
1387         gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1388
1389       i = j;
1390     }
1391
1392   free (ranks);
1393
1394   return m;
1395 }
1396
1397 static gsl_matrix *
1398 matrix_eval_row_sum (gsl_matrix *m, bool square)
1399 {
1400   if (m->size1 == 0)
1401     return m;
1402   else if (!m->size1)
1403     return gsl_matrix_alloc (0, 1);
1404
1405   gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1406   for (size_t y = 0; y < m->size1; y++)
1407     {
1408       double sum = 0;
1409       for (size_t x = 0; x < m->size2; x++)
1410         {
1411           double d = gsl_matrix_get (m, y, x);
1412           sum += square ? pow2 (d) : d;
1413         }
1414       gsl_matrix_set (result, y, 0, sum);
1415     }
1416   return result;
1417 }
1418
1419 static gsl_matrix *
1420 matrix_eval_RSSQ (gsl_matrix *m)
1421 {
1422   return matrix_eval_row_sum (m, true);
1423 }
1424
1425 static gsl_matrix *
1426 matrix_eval_RSUM (gsl_matrix *m)
1427 {
1428   return matrix_eval_row_sum (m, false);
1429 }
1430
1431 static gsl_matrix *
1432 matrix_eval_SIN (gsl_matrix *m)
1433 {
1434   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1435     *d = sin (*d);
1436   return m;
1437 }
1438
1439 static gsl_matrix *
1440 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1441 {
1442   if (m1->size1 != m2->size1)
1443     {
1444       msg (SE, _("SOLVE requires its arguments to have the same number of "
1445                  "rows, but the first argument has dimensions %zu×%zu and "
1446                  "the second %zu×%zu."),
1447            m1->size1, m1->size2,
1448            m2->size1, m2->size2);
1449       return NULL;
1450     }
1451
1452   gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1453   gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1454   int signum;
1455   gsl_linalg_LU_decomp (m1, p, &signum);
1456   for (size_t i = 0; i < m2->size2; i++)
1457     {
1458       gsl_vector bi = gsl_matrix_column (m2, i).vector;
1459       gsl_vector xi = gsl_matrix_column (x, i).vector;
1460       gsl_linalg_LU_solve (m1, p, &bi, &xi);
1461     }
1462   gsl_permutation_free (p);
1463   return x;
1464 }
1465
1466 static gsl_matrix *
1467 matrix_eval_SQRT (gsl_matrix *m)
1468 {
1469   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1470     {
1471       if (*d < 0)
1472         {
1473           msg (SE, _("Argument to SQRT must be nonnegative."));
1474           return NULL;
1475         }
1476       *d = sqrt (*d);
1477     }
1478   return m;
1479 }
1480
1481 static gsl_matrix *
1482 matrix_eval_SSCP (gsl_matrix *m)
1483 {
1484   gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1485   gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1486   return sscp;
1487 }
1488
1489 static gsl_matrix *
1490 matrix_eval_SVAL (gsl_matrix *m)
1491 {
1492   gsl_matrix *tmp_mat = NULL;
1493   if (m->size2 > m->size1)
1494     {
1495       tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1496       gsl_matrix_transpose_memcpy (tmp_mat, m);
1497       m = tmp_mat;
1498     }
1499
1500   /* Do SVD. */
1501   gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1502   gsl_vector *S = gsl_vector_alloc (m->size2);
1503   gsl_vector *work = gsl_vector_alloc (m->size2);
1504   gsl_linalg_SV_decomp (m, V, S, work);
1505
1506   gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1507   for (size_t i = 0; i < m->size2; i++)
1508     gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1509
1510   gsl_matrix_free (V);
1511   gsl_vector_free (S);
1512   gsl_vector_free (work);
1513   gsl_matrix_free (tmp_mat);
1514
1515   return vals;
1516 }
1517
1518 static gsl_matrix *
1519 matrix_eval_SWEEP (gsl_matrix *m, double d)
1520 {
1521   if (d < 1 || d > SIZE_MAX)
1522     {
1523       msg (SE, _("Scalar argument to SWEEP must be integer."));
1524       return NULL;
1525     }
1526   size_t k = d - 1;
1527   if (k >= MIN (m->size1, m->size2))
1528     {
1529       msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1530                  "equal to the smaller of the matrix argument's rows and "
1531                  "columns."));
1532       return NULL;
1533     }
1534
1535   double m_kk = gsl_matrix_get (m, k, k);
1536   if (fabs (m_kk) > 1e-19)
1537     {
1538       gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1539       MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1540         {
1541           double m_ij = gsl_matrix_get (m, i, j);
1542           double m_ik = gsl_matrix_get (m, i, k);
1543           double m_kj = gsl_matrix_get (m, k, j);
1544           *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1545                    : i != k ? -m_ik
1546                    : j != k ? m_kj
1547                    : 1.0) / m_kk;
1548         }
1549       return a;
1550     }
1551   else
1552     {
1553       for (size_t i = 0; i < m->size1; i++)
1554         {
1555           gsl_matrix_set (m, i, k, 0);
1556           gsl_matrix_set (m, k, i, 0);
1557         }
1558       return m;
1559     }
1560 }
1561
1562 static double
1563 matrix_eval_TRACE (gsl_matrix *m)
1564 {
1565   double sum = 0;
1566   size_t n = MIN (m->size1, m->size2);
1567   for (size_t i = 0; i < n; i++)
1568     sum += gsl_matrix_get (m, i, i);
1569   return sum;
1570 }
1571
1572 static gsl_matrix *
1573 matrix_eval_T (gsl_matrix *m)
1574 {
1575   return matrix_eval_TRANSPOS (m);
1576 }
1577
1578 static gsl_matrix *
1579 matrix_eval_TRANSPOS (gsl_matrix *m)
1580 {
1581   if (m->size1 == m->size2)
1582     {
1583       gsl_matrix_transpose (m);
1584       return m;
1585     }
1586   else
1587     {
1588       gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1589       gsl_matrix_transpose_memcpy (t, m);
1590       return t;
1591     }
1592 }
1593
1594 static gsl_matrix *
1595 matrix_eval_TRUNC (gsl_matrix *m)
1596 {
1597   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1598     *d = trunc (*d);
1599   return m;
1600 }
1601
1602 static gsl_matrix *
1603 matrix_eval_UNIFORM (double r_, double c_)
1604 {
1605   if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1606     {
1607       msg (SE, _("Arguments to UNIFORM must be integers."));
1608       return NULL;
1609     }
1610   size_t r = r_;
1611   size_t c = c_;
1612   if (size_overflow_p (xtimes (r, xmax (c, 1))))
1613     {
1614       msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1615       return NULL;
1616     }
1617
1618   gsl_matrix *m = gsl_matrix_alloc (r, c);
1619   MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1620     *d = gsl_ran_flat (get_rng (), 0, 1);
1621   return m;
1622 }
1623
1624 struct matrix_function
1625   {
1626     const char *name;
1627     enum matrix_op op;
1628     size_t min_args, max_args;
1629   };
1630
1631 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1632
1633 static const struct matrix_function *
1634 matrix_parse_function_name (const char *token)
1635 {
1636   static const struct matrix_function functions[] =
1637     {
1638 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1639       MATRIX_FUNCTIONS
1640 #undef F
1641     };
1642   enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1643
1644   for (size_t i = 0; i < N_FUNCTIONS; i++)
1645     {
1646       if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1647         return &functions[i];
1648     }
1649   return NULL;
1650 }
1651
1652 static struct read_file *
1653 read_file_create (struct matrix_state *s, struct file_handle *fh)
1654 {
1655   for (size_t i = 0; i < s->n_read_files; i++)
1656     {
1657       struct read_file *rf = s->read_files[i];
1658       if (rf->file == fh)
1659         {
1660           fh_unref (fh);
1661           return rf;
1662         }
1663     }
1664
1665   struct read_file *rf = xmalloc (sizeof *rf);
1666   *rf = (struct read_file) { .file = fh };
1667
1668   s->read_files = xrealloc (s->read_files,
1669                             (s->n_read_files + 1) * sizeof *s->read_files);
1670   s->read_files[s->n_read_files++] = rf;
1671   return rf;
1672 }
1673
1674 static struct dfm_reader *
1675 read_file_open (struct read_file *rf)
1676 {
1677   if (!rf->reader)
1678     rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1679   return rf->reader;
1680 }
1681
1682 static void
1683 read_file_destroy (struct read_file *rf)
1684 {
1685   if (rf)
1686     {
1687       fh_unref (rf->file);
1688       dfm_close_reader (rf->reader);
1689       free (rf->encoding);
1690       free (rf);
1691     }
1692 }
1693
1694 static bool
1695 matrix_parse_function (struct matrix_state *s, const char *token,
1696                        struct matrix_expr **exprp)
1697 {
1698   *exprp = NULL;
1699   if (lex_next_token (s->lexer, 1) != T_LPAREN)
1700     return false;
1701
1702   if (lex_match_id (s->lexer, "EOF"))
1703     {
1704       lex_get (s->lexer);
1705       struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1706       if (!fh)
1707         return true;
1708
1709       if (!lex_force_match (s->lexer, T_RPAREN))
1710         {
1711           fh_unref (fh);
1712           return true;
1713         }
1714
1715       struct read_file *rf = read_file_create (s, fh);
1716
1717       struct matrix_expr *e = xmalloc (sizeof *e);
1718       *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1719       *exprp = e;
1720       return true;
1721     }
1722
1723   const struct matrix_function *f = matrix_parse_function_name (token);
1724   if (!f)
1725     return false;
1726
1727   lex_get_n (s->lexer, 2);
1728
1729   struct matrix_expr *e = xmalloc (sizeof *e);
1730   *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1731
1732   size_t allocated_subs = 0;
1733   do
1734     {
1735       struct matrix_expr *sub = matrix_parse_expr (s);
1736       if (!sub)
1737         goto error;
1738
1739       if (e->n_subs >= allocated_subs)
1740         e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1741       e->subs[e->n_subs++] = sub;
1742     }
1743   while (lex_match (s->lexer, T_COMMA));
1744   if (!lex_force_match (s->lexer, T_RPAREN))
1745     goto error;
1746
1747   *exprp = e;
1748   return true;
1749
1750 error:
1751   matrix_expr_destroy (e);
1752   return true;
1753 }
1754
1755 static struct matrix_expr *
1756 matrix_parse_primary (struct matrix_state *s)
1757 {
1758   if (lex_is_number (s->lexer))
1759     {
1760       double number = lex_number (s->lexer);
1761       lex_get (s->lexer);
1762       return matrix_expr_create_number (number);
1763     }
1764   else if (lex_is_string (s->lexer))
1765     {
1766       char string[sizeof (double)];
1767       buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1768       lex_get (s->lexer);
1769
1770       double number;
1771       memcpy (&number, string, sizeof number);
1772       return matrix_expr_create_number (number);
1773     }
1774   else if (lex_match (s->lexer, T_LPAREN))
1775     {
1776       struct matrix_expr *e = matrix_parse_or_xor (s);
1777       if (!e || !lex_force_match (s->lexer, T_RPAREN))
1778         {
1779           matrix_expr_destroy (e);
1780           return NULL;
1781         }
1782       return e;
1783     }
1784   else if (lex_match (s->lexer, T_LCURLY))
1785     {
1786       struct matrix_expr *e = matrix_parse_curly_semi (s);
1787       if (!e || !lex_force_match (s->lexer, T_RCURLY))
1788         {
1789           matrix_expr_destroy (e);
1790           return NULL;
1791         }
1792       return e;
1793     }
1794   else if (lex_token (s->lexer) == T_ID)
1795     {
1796       struct matrix_expr *retval;
1797       if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1798         return retval;
1799
1800       struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1801       if (!var)
1802         {
1803           lex_error (s->lexer, _("Unknown variable %s."),
1804                      lex_tokcstr (s->lexer));
1805           return NULL;
1806         }
1807       lex_get (s->lexer);
1808
1809       struct matrix_expr *e = xmalloc (sizeof *e);
1810       *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1811       return e;
1812     }
1813   else if (lex_token (s->lexer) == T_ALL)
1814     {
1815       struct matrix_expr *retval;
1816       if (matrix_parse_function (s, "ALL", &retval))
1817         return retval;
1818     }
1819
1820   lex_error (s->lexer, NULL);
1821   return NULL;
1822 }
1823
1824 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1825
1826 static bool
1827 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1828 {
1829   if (lex_match (s->lexer, T_COLON))
1830     {
1831       *indexp = NULL;
1832       return true;
1833     }
1834   else
1835     {
1836       *indexp = matrix_parse_expr (s);
1837       return *indexp != NULL;
1838     }
1839 }
1840
1841 static struct matrix_expr *
1842 matrix_parse_postfix (struct matrix_state *s)
1843 {
1844   struct matrix_expr *lhs = matrix_parse_primary (s);
1845   if (!lhs || !lex_match (s->lexer, T_LPAREN))
1846     return lhs;
1847
1848   struct matrix_expr *i0;
1849   if (!matrix_parse_index_expr (s, &i0))
1850     {
1851       matrix_expr_destroy (lhs);
1852       return NULL;
1853     }
1854   if (lex_match (s->lexer, T_RPAREN))
1855     return (i0
1856             ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1857             : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1858   else if (lex_match (s->lexer, T_COMMA))
1859     {
1860       struct matrix_expr *i1;
1861       if (!matrix_parse_index_expr (s, &i1)
1862           || !lex_force_match (s->lexer, T_RPAREN))
1863         {
1864           matrix_expr_destroy (lhs);
1865           matrix_expr_destroy (i0);
1866           matrix_expr_destroy (i1);
1867           return NULL;
1868         }
1869       return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1870               : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1871               : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1872               : lhs);
1873     }
1874   else
1875     {
1876       lex_error_expecting (s->lexer, "`)'", "`,'");
1877       return NULL;
1878     }
1879 }
1880
1881 static struct matrix_expr *
1882 matrix_parse_unary (struct matrix_state *s)
1883 {
1884   if (lex_match (s->lexer, T_DASH))
1885     {
1886       struct matrix_expr *lhs = matrix_parse_unary (s);
1887       return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1888     }
1889   else if (lex_match (s->lexer, T_PLUS))
1890     return matrix_parse_unary (s);
1891   else
1892     return matrix_parse_postfix (s);
1893 }
1894
1895 static struct matrix_expr *
1896 matrix_parse_seq (struct matrix_state *s)
1897 {
1898   struct matrix_expr *start = matrix_parse_unary (s);
1899   if (!start || !lex_match (s->lexer, T_COLON))
1900     return start;
1901
1902   struct matrix_expr *end = matrix_parse_unary (s);
1903   if (!end)
1904     {
1905       matrix_expr_destroy (start);
1906       return NULL;
1907     }
1908
1909   if (lex_match (s->lexer, T_COLON))
1910     {
1911       struct matrix_expr *increment = matrix_parse_unary (s);
1912       if (!increment)
1913         {
1914           matrix_expr_destroy (start);
1915           matrix_expr_destroy (end);
1916           return NULL;
1917         }
1918       return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1919     }
1920   else
1921     return matrix_expr_create_2 (MOP_SEQ, start, end);
1922 }
1923
1924 static struct matrix_expr *
1925 matrix_parse_exp (struct matrix_state *s)
1926 {
1927   struct matrix_expr *lhs = matrix_parse_seq (s);
1928   if (!lhs)
1929     return NULL;
1930
1931   for (;;)
1932     {
1933       enum matrix_op op;
1934       if (lex_match (s->lexer, T_EXP))
1935         op = MOP_EXP_MAT;
1936       else if (lex_match_phrase (s->lexer, "&**"))
1937         op = MOP_EXP_ELEMS;
1938       else
1939         return lhs;
1940
1941       struct matrix_expr *rhs = matrix_parse_seq (s);
1942       if (!rhs)
1943         {
1944           matrix_expr_destroy (lhs);
1945           return NULL;
1946         }
1947       lhs = matrix_expr_create_2 (op, lhs, rhs);
1948     }
1949 }
1950
1951 static struct matrix_expr *
1952 matrix_parse_mul_div (struct matrix_state *s)
1953 {
1954   struct matrix_expr *lhs = matrix_parse_exp (s);
1955   if (!lhs)
1956     return NULL;
1957
1958   for (;;)
1959     {
1960       enum matrix_op op;
1961       if (lex_match (s->lexer, T_ASTERISK))
1962         op = MOP_MUL_MAT;
1963       else if (lex_match (s->lexer, T_SLASH))
1964         op = MOP_DIV_ELEMS;
1965       else if (lex_match_phrase (s->lexer, "&*"))
1966         op = MOP_MUL_ELEMS;
1967       else if (lex_match_phrase (s->lexer, "&/"))
1968         op = MOP_DIV_ELEMS;
1969       else
1970         return lhs;
1971
1972       struct matrix_expr *rhs = matrix_parse_exp (s);
1973       if (!rhs)
1974         {
1975           matrix_expr_destroy (lhs);
1976           return NULL;
1977         }
1978       lhs = matrix_expr_create_2 (op, lhs, rhs);
1979     }
1980 }
1981
1982 static struct matrix_expr *
1983 matrix_parse_add_sub (struct matrix_state *s)
1984 {
1985   struct matrix_expr *lhs = matrix_parse_mul_div (s);
1986   if (!lhs)
1987     return NULL;
1988
1989   for (;;)
1990     {
1991       enum matrix_op op;
1992       if (lex_match (s->lexer, T_PLUS))
1993         op = MOP_ADD_ELEMS;
1994       else if (lex_match (s->lexer, T_DASH))
1995         op = MOP_SUB_ELEMS;
1996       else if (lex_token (s->lexer) == T_NEG_NUM)
1997         op = MOP_ADD_ELEMS;
1998       else
1999         return lhs;
2000
2001       struct matrix_expr *rhs = matrix_parse_mul_div (s);
2002       if (!rhs)
2003         {
2004           matrix_expr_destroy (lhs);
2005           return NULL;
2006         }
2007       lhs = matrix_expr_create_2 (op, lhs, rhs);
2008     }
2009 }
2010
2011 static struct matrix_expr *
2012 matrix_parse_relational (struct matrix_state *s)
2013 {
2014   struct matrix_expr *lhs = matrix_parse_add_sub (s);
2015   if (!lhs)
2016     return NULL;
2017
2018   for (;;)
2019     {
2020       enum matrix_op op;
2021       if (lex_match (s->lexer, T_GT))
2022         op = MOP_GT;
2023       else if (lex_match (s->lexer, T_GE))
2024         op = MOP_GE;
2025       else if (lex_match (s->lexer, T_LT))
2026         op = MOP_LT;
2027       else if (lex_match (s->lexer, T_LE))
2028         op = MOP_LE;
2029       else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2030         op = MOP_EQ;
2031       else if (lex_match (s->lexer, T_NE))
2032         op = MOP_NE;
2033       else
2034         return lhs;
2035
2036       struct matrix_expr *rhs = matrix_parse_add_sub (s);
2037       if (!rhs)
2038         {
2039           matrix_expr_destroy (lhs);
2040           return NULL;
2041         }
2042       lhs = matrix_expr_create_2 (op, lhs, rhs);
2043     }
2044 }
2045
2046 static struct matrix_expr *
2047 matrix_parse_not (struct matrix_state *s)
2048 {
2049   if (lex_match (s->lexer, T_NOT))
2050     {
2051       struct matrix_expr *lhs = matrix_parse_not (s);
2052       return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2053     }
2054   else
2055     return matrix_parse_relational (s);
2056 }
2057
2058 static struct matrix_expr *
2059 matrix_parse_and (struct matrix_state *s)
2060 {
2061   struct matrix_expr *lhs = matrix_parse_not (s);
2062   if (!lhs)
2063     return NULL;
2064
2065   while (lex_match (s->lexer, T_AND))
2066     {
2067       struct matrix_expr *rhs = matrix_parse_not (s);
2068       if (!rhs)
2069         {
2070           matrix_expr_destroy (lhs);
2071           return NULL;
2072         }
2073       lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2074     }
2075   return lhs;
2076 }
2077
2078 static struct matrix_expr *
2079 matrix_parse_or_xor (struct matrix_state *s)
2080 {
2081   struct matrix_expr *lhs = matrix_parse_and (s);
2082   if (!lhs)
2083     return NULL;
2084
2085   for (;;)
2086     {
2087       enum matrix_op op;
2088       if (lex_match (s->lexer, T_OR))
2089         op = MOP_OR;
2090       else if (lex_match_id (s->lexer, "XOR"))
2091         op = MOP_XOR;
2092       else
2093         return lhs;
2094
2095       struct matrix_expr *rhs = matrix_parse_and (s);
2096       if (!rhs)
2097         {
2098           matrix_expr_destroy (lhs);
2099           return NULL;
2100         }
2101       lhs = matrix_expr_create_2 (op, lhs, rhs);
2102     }
2103 }
2104
2105 static struct matrix_expr *
2106 matrix_parse_expr (struct matrix_state *s)
2107 {
2108   return matrix_parse_or_xor (s);
2109 }
2110 \f
2111 /* Expression evaluation. */
2112
2113 static double
2114 matrix_op_eval (enum matrix_op op, double a, double b)
2115 {
2116   switch (op)
2117     {
2118     case MOP_ADD_ELEMS: return a + b;
2119     case MOP_SUB_ELEMS: return a - b;
2120     case MOP_MUL_ELEMS: return a * b;
2121     case MOP_DIV_ELEMS: return a / b;
2122     case MOP_EXP_ELEMS: return pow (a, b);
2123     case MOP_GT: return a > b;
2124     case MOP_GE: return a >= b;
2125     case MOP_LT: return a < b;
2126     case MOP_LE: return a <= b;
2127     case MOP_EQ: return a == b;
2128     case MOP_NE: return a != b;
2129     case MOP_AND: return (a > 0) && (b > 0);
2130     case MOP_OR: return (a > 0) || (b > 0);
2131     case MOP_XOR: return (a > 0) != (b > 0);
2132
2133 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2134       MATRIX_FUNCTIONS
2135 #undef F
2136     case MOP_NEGATE:
2137     case MOP_SEQ:
2138     case MOP_SEQ_BY:
2139     case MOP_MUL_MAT:
2140     case MOP_EXP_MAT:
2141     case MOP_NOT:
2142     case MOP_PASTE_HORZ:
2143     case MOP_PASTE_VERT:
2144     case MOP_EMPTY:
2145     case MOP_VEC_INDEX:
2146     case MOP_VEC_ALL:
2147     case MOP_MAT_INDEX:
2148     case MOP_ROW_INDEX:
2149     case MOP_COL_INDEX:
2150     case MOP_NUMBER:
2151     case MOP_VARIABLE:
2152     case MOP_EOF:
2153       NOT_REACHED ();
2154     }
2155   NOT_REACHED ();
2156 }
2157
2158 static const char *
2159 matrix_op_name (enum matrix_op op)
2160 {
2161   switch (op)
2162     {
2163     case MOP_ADD_ELEMS: return "+";
2164     case MOP_SUB_ELEMS: return "-";
2165     case MOP_MUL_ELEMS: return "&*";
2166     case MOP_DIV_ELEMS: return "&/";
2167     case MOP_EXP_ELEMS: return "&**";
2168     case MOP_GT: return ">";
2169     case MOP_GE: return ">=";
2170     case MOP_LT: return "<";
2171     case MOP_LE: return "<=";
2172     case MOP_EQ: return "=";
2173     case MOP_NE: return "<>";
2174     case MOP_AND: return "AND";
2175     case MOP_OR: return "OR";
2176     case MOP_XOR: return "XOR";
2177
2178 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2179       MATRIX_FUNCTIONS
2180 #undef F
2181     case MOP_NEGATE:
2182     case MOP_SEQ:
2183     case MOP_SEQ_BY:
2184     case MOP_MUL_MAT:
2185     case MOP_EXP_MAT:
2186     case MOP_NOT:
2187     case MOP_PASTE_HORZ:
2188     case MOP_PASTE_VERT:
2189     case MOP_EMPTY:
2190     case MOP_VEC_INDEX:
2191     case MOP_VEC_ALL:
2192     case MOP_MAT_INDEX:
2193     case MOP_ROW_INDEX:
2194     case MOP_COL_INDEX:
2195     case MOP_NUMBER:
2196     case MOP_VARIABLE:
2197     case MOP_EOF:
2198       NOT_REACHED ();
2199     }
2200   NOT_REACHED ();
2201 }
2202
2203 static bool
2204 is_scalar (const gsl_matrix *m)
2205 {
2206   return m->size1 == 1 && m->size2 == 1;
2207 }
2208
2209 static double
2210 to_scalar (const gsl_matrix *m)
2211 {
2212   assert (is_scalar (m));
2213   return gsl_matrix_get (m, 0, 0);
2214 }
2215
2216 static gsl_matrix *
2217 matrix_expr_evaluate_elementwise (enum matrix_op op,
2218                                   gsl_matrix *a, gsl_matrix *b)
2219 {
2220   if (is_scalar (b))
2221     {
2222       double be = to_scalar (b);
2223       for (size_t r = 0; r < a->size1; r++)
2224         for (size_t c = 0; c < a->size2; c++)
2225           {
2226             double *ae = gsl_matrix_ptr (a, r, c);
2227             *ae = matrix_op_eval (op, *ae, be);
2228           }
2229       return a;
2230     }
2231   else if (is_scalar (a))
2232     {
2233       double ae = to_scalar (a);
2234       for (size_t r = 0; r < b->size1; r++)
2235         for (size_t c = 0; c < b->size2; c++)
2236           {
2237             double *be = gsl_matrix_ptr (b, r, c);
2238             *be = matrix_op_eval (op, ae, *be);
2239           }
2240       return b;
2241     }
2242   else if (a->size1 == b->size1 && a->size2 == b->size2)
2243     {
2244       for (size_t r = 0; r < a->size1; r++)
2245         for (size_t c = 0; c < a->size2; c++)
2246           {
2247             double *ae = gsl_matrix_ptr (a, r, c);
2248             double be = gsl_matrix_get (b, r, c);
2249             *ae = matrix_op_eval (op, *ae, be);
2250           }
2251       return a;
2252     }
2253   else
2254     {
2255       msg (SE, _("Operands to %s must have the same dimensions or one "
2256                  "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2257            matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2258       return NULL;
2259     }
2260 }
2261
2262 static gsl_matrix *
2263 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2264 {
2265   if (is_scalar (a) || is_scalar (b))
2266     return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2267
2268   if (a->size2 != b->size1)
2269     {
2270       msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2271                  "not conformable for multiplication."),
2272            a->size1, a->size2, b->size1, b->size2);
2273       return NULL;
2274     }
2275
2276   gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2277   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2278   return c;
2279 }
2280
2281 static void
2282 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2283 {
2284   gsl_matrix *tmp = *a;
2285   *a = *b;
2286   *b = tmp;
2287 }
2288
2289 static void
2290 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2291             gsl_matrix **tmp)
2292 {
2293   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2294   swap_matrix (z, tmp);
2295 }
2296
2297 static void
2298 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2299 {
2300   mul_matrix (x, *x, *x, tmp);
2301 }
2302
2303 static gsl_matrix *
2304 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2305 {
2306   gsl_matrix *x = x_;
2307   if (x->size1 != x->size2)
2308     {
2309       msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2310                  "the left-hand size, not one with dimensions %zu×%zu."),
2311            x->size1, x->size2);
2312       return NULL;
2313     }
2314   if (!is_scalar (b))
2315     {
2316       msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2317                  "right-hand side, not a matrix with dimensions %zu×%zu."),
2318            b->size1, b->size2);
2319       return NULL;
2320     }
2321   double bf = to_scalar (b);
2322   if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2323     {
2324       msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2325                  "or outside the valid range."), bf);
2326       return NULL;
2327     }
2328   long int bl = bf;
2329
2330   gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2331   gsl_matrix *y = y_;
2332   gsl_matrix_set_identity (y);
2333   if (bl == 0)
2334     return y;
2335
2336   gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2337   gsl_matrix *t = t_;
2338   for (unsigned long int n = labs (bl); n > 1; n /= 2)
2339     if (n & 1)
2340       {
2341         mul_matrix (&y, x, y, &t);
2342         square_matrix (&x, &t);
2343       }
2344     else
2345       square_matrix (&x, &t);
2346
2347   mul_matrix (&y, x, y, &t);
2348   if (bf < 0)
2349     invert_matrix (y);
2350
2351   /* Garbage collection.
2352
2353      There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2354      a permutation of them.  We are returning one of them; that one must not be
2355      destroyed.  We must not destroy 'x_' because the caller owns it. */
2356   if (y != y_)
2357     gsl_matrix_free (y_);
2358   if (y != t_)
2359     gsl_matrix_free (t_);
2360
2361   return y;
2362 }
2363
2364 static gsl_matrix *
2365 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2366                           gsl_matrix *by_)
2367 {
2368   if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2369     {
2370       msg (SE, _("All operands of : operator must be scalars."));
2371       return NULL;
2372     }
2373
2374   long int start = to_scalar (start_);
2375   long int end = to_scalar (end_);
2376   long int by = by_ ? to_scalar (by_) : 1;
2377
2378   if (!by)
2379     {
2380       msg (SE, _("The increment operand to : must be nonzero."));
2381       return NULL;
2382     }
2383
2384   long int n = (end >= start && by > 0 ? (end - start + by) / by
2385                 : end < start && by < 0 ? (start - end - by) / -by
2386                 : 0);
2387   gsl_matrix *m = gsl_matrix_alloc (1, n);
2388   for (long int i = 0; i < n; i++)
2389     gsl_matrix_set (m, 0, i, start + i * by);
2390   return m;
2391 }
2392
2393 static gsl_matrix *
2394 matrix_expr_evaluate_not (gsl_matrix *a)
2395 {
2396   for (size_t r = 0; r < a->size1; r++)
2397     for (size_t c = 0; c < a->size2; c++)
2398       {
2399         double *ae = gsl_matrix_ptr (a, r, c);
2400         *ae = !(*ae > 0);
2401       }
2402   return a;
2403 }
2404
2405 static gsl_matrix *
2406 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2407 {
2408   if (a->size1 != b->size1)
2409     {
2410       if (!a->size1 || !a->size2)
2411         return b;
2412       else if (!b->size1 || !b->size2)
2413         return a;
2414
2415       msg (SE, _("All columns in a matrix must have the same number of rows, "
2416                  "but this tries to paste matrices with %zu and %zu rows."),
2417            a->size1, b->size1);
2418       return NULL;
2419     }
2420
2421   gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2422   for (size_t y = 0; y < a->size1; y++)
2423     {
2424       for (size_t x = 0; x < a->size2; x++)
2425         gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2426       for (size_t x = 0; x < b->size2; x++)
2427         gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2428     }
2429   return c;
2430 }
2431
2432 static gsl_matrix *
2433 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2434 {
2435   if (a->size2 != b->size2)
2436     {
2437       if (!a->size1 || !a->size2)
2438         return b;
2439       else if (!b->size1 || !b->size2)
2440         return a;
2441
2442       msg (SE, _("All rows in a matrix must have the same number of columns, "
2443                  "but this tries to stack matrices with %zu and %zu columns."),
2444            a->size2, b->size2);
2445       return NULL;
2446     }
2447
2448   gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2449   for (size_t x = 0; x < a->size2; x++)
2450     {
2451       for (size_t y = 0; y < a->size1; y++)
2452         gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2453       for (size_t y = 0; y < b->size1; y++)
2454         gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2455     }
2456   return c;
2457 }
2458
2459 static gsl_vector *
2460 matrix_to_vector (gsl_matrix *m)
2461 {
2462   assert (m->owner);
2463   gsl_vector v = to_vector (m);
2464   assert (v.block == m->block || !v.block);
2465   assert (!v.owner);
2466   v.owner = 1;
2467   m->owner = 0;
2468   return xmemdup (&v, sizeof v);
2469 }
2470
2471 enum index_type {
2472   IV_ROW,
2473   IV_COLUMN,
2474   IV_VECTOR
2475 };
2476
2477 struct index_vector
2478   {
2479     size_t *indexes;
2480     size_t n;
2481   };
2482 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2483
2484 static bool
2485 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2486                                enum index_type index_type, size_t other_size,
2487                                struct index_vector *iv)
2488 {
2489   if (m)
2490     {
2491       if (!is_vector (m))
2492         {
2493           switch (index_type)
2494             {
2495             case IV_VECTOR:
2496               msg (SE, _("Vector index must be scalar or vector, not a "
2497                          "%zu×%zu matrix."),
2498                    m->size1, m->size2);
2499               break;
2500
2501             case IV_ROW:
2502               msg (SE, _("Matrix row index must be scalar or vector, not a "
2503                          "%zu×%zu matrix."),
2504                    m->size1, m->size2);
2505               break;
2506
2507             case IV_COLUMN:
2508               msg (SE, _("Matrix column index must be scalar or vector, not a "
2509                          "%zu×%zu matrix."),
2510                    m->size1, m->size2);
2511               break;
2512             }
2513           return false;
2514         }
2515
2516       gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2517       *iv = (struct index_vector) {
2518         .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2519         .n = v.size,
2520       };
2521       for (size_t i = 0; i < v.size; i++)
2522         {
2523           double index = gsl_vector_get (&v, i);
2524           if (index < 1 || index >= size + 1)
2525             {
2526               switch (index_type)
2527                 {
2528                 case IV_VECTOR:
2529                   msg (SE, _("Index %g is out of range for vector "
2530                              "with %zu elements."), index, size);
2531                   break;
2532
2533                 case IV_ROW:
2534                   msg (SE, _("%g is not a valid row index for "
2535                              "a %zu×%zu matrix."),
2536                        index, size, other_size);
2537                   break;
2538
2539                 case IV_COLUMN:
2540                   msg (SE, _("%g is not a valid column index for "
2541                              "a %zu×%zu matrix."),
2542                        index, other_size, size);
2543                   break;
2544                 }
2545
2546               free (iv->indexes);
2547               return false;
2548             }
2549           iv->indexes[i] = index - 1;
2550         }
2551       return true;
2552     }
2553   else
2554     {
2555       *iv = (struct index_vector) {
2556         .indexes = xnmalloc (size, sizeof *iv->indexes),
2557         .n = size,
2558       };
2559       for (size_t i = 0; i < size; i++)
2560         iv->indexes[i] = i;
2561       return true;
2562     }
2563 }
2564
2565 static gsl_matrix *
2566 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2567 {
2568   if (!is_vector (sm))
2569     {
2570       msg (SE, _("Vector index operator may not be applied to "
2571                  "a %zu×%zu matrix."),
2572            sm->size1, sm->size2);
2573       return NULL;
2574     }
2575
2576   return sm;
2577 }
2578
2579 static gsl_matrix *
2580 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2581 {
2582   if (!matrix_expr_evaluate_vec_all (sm))
2583     return NULL;
2584
2585   gsl_vector sv = to_vector (sm);
2586   struct index_vector iv;
2587   if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2588     return NULL;
2589
2590   gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2591                                      sm->size1 == 1 ? iv.n : 1);
2592   gsl_vector dv = to_vector (dm);
2593   for (size_t dx = 0; dx < iv.n; dx++)
2594     {
2595       size_t sx = iv.indexes[dx];
2596       gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2597     }
2598   free (iv.indexes);
2599
2600   return dm;
2601 }
2602
2603 static gsl_matrix *
2604 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2605                                 gsl_matrix *im1)
2606 {
2607   struct index_vector iv0;
2608   if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2609     return NULL;
2610
2611   struct index_vector iv1;
2612   if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2613                                       &iv1))
2614     {
2615       free (iv0.indexes);
2616       return NULL;
2617     }
2618
2619   gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2620   for (size_t dy = 0; dy < iv0.n; dy++)
2621     {
2622       size_t sy = iv0.indexes[dy];
2623
2624       for (size_t dx = 0; dx < iv1.n; dx++)
2625         {
2626           size_t sx = iv1.indexes[dx];
2627           gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2628         }
2629     }
2630   free (iv0.indexes);
2631   free (iv1.indexes);
2632   return dm;
2633 }
2634
2635 #define F(NAME, PROTOTYPE)                              \
2636   static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2637     const char *name, gsl_matrix *[], size_t,           \
2638     matrix_proto_##PROTOTYPE *);
2639 MATRIX_FUNCTIONS
2640 #undef F
2641
2642 static bool
2643 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2644 {
2645   if (!is_scalar (subs[index]))
2646     {
2647       msg (SE, _("Function %s argument %zu must be a scalar, "
2648                  "not a %zu×%zu matrix."),
2649            name, index + 1, subs[index]->size1, subs[index]->size2);
2650       return false;
2651     }
2652   return true;
2653 }
2654
2655 static bool
2656 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2657 {
2658   if (!is_vector (subs[index]))
2659     {
2660       msg (SE, _("Function %s argument %zu must be a vector, "
2661                  "not a %zu×%zu matrix."),
2662            name, index + 1, subs[index]->size1, subs[index]->size2);
2663       return false;
2664     }
2665   return true;
2666 }
2667
2668 static bool
2669 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2670 {
2671   for (size_t i = 0; i < n_subs; i++)
2672     {
2673       if (!check_scalar_arg (name, subs, i))
2674         return false;
2675       d[i] = to_scalar (subs[i]);
2676     }
2677   return true;
2678 }
2679
2680 static gsl_matrix *
2681 matrix_expr_evaluate_m_d (const char *name,
2682                           gsl_matrix *subs[], size_t n_subs,
2683                           matrix_proto_m_d *f)
2684 {
2685   assert (n_subs == 1);
2686
2687   double d;
2688   return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2689 }
2690
2691 static gsl_matrix *
2692 matrix_expr_evaluate_m_dd (const char *name,
2693                            gsl_matrix *subs[], size_t n_subs,
2694                            matrix_proto_m_dd *f)
2695 {
2696   assert (n_subs == 2);
2697
2698   double d[2];
2699   return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2700 }
2701
2702 static gsl_matrix *
2703 matrix_expr_evaluate_m_ddd (const char *name,
2704                            gsl_matrix *subs[], size_t n_subs,
2705                            matrix_proto_m_ddd *f)
2706 {
2707   assert (n_subs == 3);
2708
2709   double d[3];
2710   return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2711 }
2712
2713 static gsl_matrix *
2714 matrix_expr_evaluate_m_m (const char *name UNUSED,
2715                           gsl_matrix *subs[], size_t n_subs,
2716                           matrix_proto_m_m *f)
2717 {
2718   assert (n_subs == 1);
2719   return f (subs[0]);
2720 }
2721
2722 static gsl_matrix *
2723 matrix_expr_evaluate_m_md (const char *name UNUSED,
2724                            gsl_matrix *subs[], size_t n_subs,
2725                            matrix_proto_m_md *f)
2726 {
2727   assert (n_subs == 2);
2728   if (!check_scalar_arg (name, subs, 1))
2729     return NULL;
2730   return f (subs[0], to_scalar (subs[1]));
2731 }
2732
2733 static gsl_matrix *
2734 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2735                             gsl_matrix *subs[], size_t n_subs,
2736                             matrix_proto_m_mdd *f)
2737 {
2738   assert (n_subs == 3);
2739   if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2740     return NULL;
2741   return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2742 }
2743
2744 static gsl_matrix *
2745 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2746                            gsl_matrix *subs[], size_t n_subs,
2747                            matrix_proto_m_mm *f)
2748 {
2749   assert (n_subs == 2);
2750   return f (subs[0], subs[1]);
2751 }
2752
2753 static gsl_matrix *
2754 matrix_expr_evaluate_m_v (const char *name,
2755                           gsl_matrix *subs[], size_t n_subs,
2756                           matrix_proto_m_v *f)
2757 {
2758   assert (n_subs == 1);
2759   if (!check_vector_arg (name, subs, 0))
2760     return NULL;
2761   gsl_vector v = to_vector (subs[0]);
2762   return f (&v);
2763 }
2764
2765 static gsl_matrix *
2766 matrix_expr_evaluate_d_m (const char *name UNUSED,
2767                           gsl_matrix *subs[], size_t n_subs,
2768                           matrix_proto_d_m *f)
2769 {
2770   assert (n_subs == 1);
2771
2772   gsl_matrix *m = gsl_matrix_alloc (1, 1);
2773   gsl_matrix_set (m, 0, 0, f (subs[0]));
2774   return m;
2775 }
2776
2777 static gsl_matrix *
2778 matrix_expr_evaluate_m_any (const char *name UNUSED,
2779                           gsl_matrix *subs[], size_t n_subs,
2780                           matrix_proto_m_any *f)
2781 {
2782   return f (subs, n_subs);
2783 }
2784
2785 static gsl_matrix *
2786 matrix_expr_evaluate_IDENT (const char *name,
2787                             gsl_matrix *subs[], size_t n_subs,
2788                             matrix_proto_IDENT *f)
2789 {
2790   assert (n_subs <= 2);
2791
2792   double d[2];
2793   if (!to_scalar_args (name, subs, n_subs, d))
2794     return NULL;
2795
2796   return f (d[0], d[n_subs - 1]);
2797 }
2798
2799 static gsl_matrix *
2800 matrix_expr_evaluate (const struct matrix_expr *e)
2801 {
2802   if (e->op == MOP_NUMBER)
2803     {
2804       gsl_matrix *m = gsl_matrix_alloc (1, 1);
2805       gsl_matrix_set (m, 0, 0, e->number);
2806       return m;
2807     }
2808   else if (e->op == MOP_VARIABLE)
2809     {
2810       const gsl_matrix *src = e->variable->value;
2811       if (!src)
2812         {
2813           msg (SE, _("Uninitialized variable %s used in expression."),
2814                e->variable->name);
2815           return NULL;
2816         }
2817
2818       gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2819       gsl_matrix_memcpy (dst, src);
2820       return dst;
2821     }
2822   else if (e->op == MOP_EOF)
2823     {
2824       struct dfm_reader *reader = read_file_open (e->eof);
2825       gsl_matrix *m = gsl_matrix_alloc (1, 1);
2826       gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2827       return m;
2828     }
2829
2830   enum { N_LOCAL = 3 };
2831   gsl_matrix *local_subs[N_LOCAL];
2832   gsl_matrix **subs = (e->n_subs < N_LOCAL
2833                        ? local_subs
2834                        : xmalloc (e->n_subs * sizeof *subs));
2835
2836   for (size_t i = 0; i < e->n_subs; i++)
2837     {
2838       subs[i] = matrix_expr_evaluate (e->subs[i]);
2839       if (!subs[i])
2840         {
2841           for (size_t j = 0; j < i; j++)
2842             gsl_matrix_free (subs[j]);
2843           if (subs != local_subs)
2844             free (subs);
2845           return NULL;
2846         }
2847     }
2848
2849   gsl_matrix *result = NULL;
2850   switch (e->op)
2851     {
2852 #define F(NAME, PROTOTYPE)                                              \
2853       case MOP_F_##NAME:                                                \
2854         result = matrix_expr_evaluate_##PROTOTYPE (#NAME,               \
2855                                                    subs, e->n_subs,     \
2856                                                    matrix_eval_##NAME); \
2857       break;
2858       MATRIX_FUNCTIONS
2859 #undef F
2860
2861     case MOP_NEGATE:
2862       gsl_matrix_scale (subs[0], -1.0);
2863       result = subs[0];
2864       break;
2865
2866     case MOP_ADD_ELEMS:
2867     case MOP_SUB_ELEMS:
2868     case MOP_MUL_ELEMS:
2869     case MOP_DIV_ELEMS:
2870     case MOP_EXP_ELEMS:
2871     case MOP_GT:
2872     case MOP_GE:
2873     case MOP_LT:
2874     case MOP_LE:
2875     case MOP_EQ:
2876     case MOP_NE:
2877     case MOP_AND:
2878     case MOP_OR:
2879     case MOP_XOR:
2880       result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2881       break;
2882
2883     case MOP_NOT:
2884       result = matrix_expr_evaluate_not (subs[0]);
2885       break;
2886
2887     case MOP_SEQ:
2888       result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2889       break;
2890
2891     case MOP_SEQ_BY:
2892       result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2893       break;
2894
2895     case MOP_MUL_MAT:
2896       result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2897       break;
2898
2899     case MOP_EXP_MAT:
2900       result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2901       break;
2902
2903     case MOP_PASTE_HORZ:
2904       result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2905       break;
2906
2907     case MOP_PASTE_VERT:
2908       result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2909       break;
2910
2911     case MOP_EMPTY:
2912       result = gsl_matrix_alloc (0, 0);
2913       break;
2914
2915     case MOP_VEC_INDEX:
2916       result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2917       break;
2918
2919     case MOP_VEC_ALL:
2920       result = matrix_expr_evaluate_vec_all (subs[0]);
2921       break;
2922
2923     case MOP_MAT_INDEX:
2924       result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2925       break;
2926
2927     case MOP_ROW_INDEX:
2928       result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2929       break;
2930
2931     case MOP_COL_INDEX:
2932       result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
2933       break;
2934
2935     case MOP_NUMBER:
2936     case MOP_VARIABLE:
2937     case MOP_EOF:
2938       NOT_REACHED ();
2939     }
2940
2941   for (size_t i = 0; i < e->n_subs; i++)
2942     if (subs[i] != result)
2943       gsl_matrix_free (subs[i]);
2944   if (subs != local_subs)
2945     free (subs);
2946   return result;
2947 }
2948
2949 static bool
2950 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
2951                              double *d)
2952 {
2953   gsl_matrix *m = matrix_expr_evaluate (e);
2954   if (!m)
2955     return false;
2956
2957   if (!is_scalar (m))
2958     {
2959       msg (SE, _("Expression for %s must evaluate to scalar, "
2960                  "not a %zu×%zu matrix."),
2961            context, m->size1, m->size2);
2962       gsl_matrix_free (m);
2963       return false;
2964     }
2965
2966   *d = to_scalar (m);
2967   gsl_matrix_free (m);
2968   return true;
2969 }
2970
2971 static bool
2972 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
2973                               long int *integer)
2974 {
2975   double d;
2976   if (!matrix_expr_evaluate_scalar (e, context, &d))
2977     return false;
2978
2979   d = trunc (d);
2980   if (d < LONG_MIN || d > LONG_MAX)
2981     {
2982       msg (SE, _("Expression for %s is outside the integer range."), context);
2983       return false;
2984     }
2985   *integer = d;
2986   return true;
2987 }
2988 \f
2989 struct matrix_lvalue
2990   {
2991     struct matrix_var *var;
2992     struct matrix_expr *indexes[2];
2993     size_t n_indexes;
2994   };
2995
2996 static void
2997 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
2998 {
2999   if (lvalue)
3000     {
3001       for (size_t i = 0; i < lvalue->n_indexes; i++)
3002         matrix_expr_destroy (lvalue->indexes[i]);
3003       free (lvalue);
3004     }
3005 }
3006
3007 static struct matrix_lvalue *
3008 matrix_lvalue_parse (struct matrix_state *s)
3009 {
3010   struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3011   if (!lex_force_id (s->lexer))
3012     goto error;
3013   lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3014   if (lex_next_token (s->lexer, 1) == T_LPAREN)
3015     {
3016       if (!lvalue->var)
3017         {
3018           msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3019           free (lvalue);
3020           return NULL;
3021         }
3022
3023       lex_get_n (s->lexer, 2);
3024
3025       if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3026         goto error;
3027       lvalue->n_indexes++;
3028
3029       if (lex_match (s->lexer, T_COMMA))
3030         {
3031           if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3032             goto error;
3033           lvalue->n_indexes++;
3034         }
3035       if (!lex_force_match (s->lexer, T_RPAREN))
3036         goto error;
3037     }
3038   else
3039     {
3040       if (!lvalue->var)
3041         lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3042       lex_get (s->lexer);
3043     }
3044   return lvalue;
3045
3046 error:
3047   matrix_lvalue_destroy (lvalue);
3048   return NULL;
3049 }
3050
3051 static bool
3052 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3053                                enum index_type index_type, size_t other_size,
3054                                struct index_vector *iv)
3055 {
3056   gsl_matrix *m;
3057   if (e)
3058     {
3059       m = matrix_expr_evaluate (e);
3060       if (!m)
3061         return false;
3062     }
3063   else
3064     m = NULL;
3065
3066   bool ok = matrix_normalize_index_vector (m, size, index_type,
3067                                            other_size, iv);
3068   gsl_matrix_free (m);
3069   return ok;
3070 }
3071
3072 static bool
3073 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3074                              struct index_vector *iv, gsl_matrix *sm)
3075 {
3076   gsl_vector dv = to_vector (lvalue->var->value);
3077
3078   /* Convert source matrix 'sm' to source vector 'sv'. */
3079   if (!is_vector (sm))
3080     {
3081       msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3082            sm->size1, sm->size2);
3083       return false;
3084     }
3085   gsl_vector sv = to_vector (sm);
3086   if (iv->n != sv.size)
3087     {
3088       msg (SE, _("Can't assign %zu-element vector "
3089                  "to %zu-element subvector."), sv.size, iv->n);
3090       return false;
3091     }
3092
3093   /* Assign elements. */
3094   for (size_t x = 0; x < iv->n; x++)
3095     gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3096   return true;
3097 }
3098
3099 static bool
3100 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3101                              struct index_vector *iv0,
3102                              struct index_vector *iv1, gsl_matrix *sm)
3103 {
3104   gsl_matrix *dm = lvalue->var->value;
3105
3106   /* Convert source matrix 'sm' to source vector 'sv'. */
3107   if (iv0->n != sm->size1)
3108     {
3109       msg (SE, _("Row index vector for assignment to %s has %zu elements "
3110                  "but source matrix has %zu rows."),
3111            lvalue->var->name, iv0->n, sm->size1);
3112       return false;
3113     }
3114   if (iv1->n != sm->size2)
3115     {
3116       msg (SE, _("Column index vector for assignment to %s has %zu elements "
3117                  "but source matrix has %zu columns."),
3118            lvalue->var->name, iv1->n, sm->size2);
3119       return false;
3120     }
3121
3122   /* Assign elements. */
3123   for (size_t y = 0; y < iv0->n; y++)
3124     {
3125       size_t f0 = iv0->indexes[y];
3126       for (size_t x = 0; x < iv1->n; x++)
3127         {
3128           size_t f1 = iv1->indexes[x];
3129           gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3130         }
3131     }
3132   return true;
3133 }
3134
3135 /* Takes ownership of M. */
3136 static bool
3137 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3138                       struct index_vector *iv0, struct index_vector *iv1,
3139                       gsl_matrix *sm)
3140 {
3141   if (!lvalue->n_indexes)
3142     {
3143       gsl_matrix_free (lvalue->var->value);
3144       lvalue->var->value = sm;
3145       return true;
3146     }
3147   else
3148     {
3149       bool ok = (lvalue->n_indexes == 1
3150                  ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3151                  : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3152       free (iv0->indexes);
3153       free (iv1->indexes);
3154       gsl_matrix_free (sm);
3155       return ok;
3156     }
3157 }
3158
3159 static bool
3160 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3161                         struct index_vector *iv0,
3162                         struct index_vector *iv1)
3163 {
3164   *iv0 = INDEX_VECTOR_INIT;
3165   *iv1 = INDEX_VECTOR_INIT;
3166   if (!lvalue->n_indexes)
3167     return true;
3168
3169   /* Validate destination matrix exists and has the right shape. */
3170   gsl_matrix *dm = lvalue->var->value;
3171   if (!dm)
3172     {
3173       msg (SE, _("Undefined variable %s."), lvalue->var->name);
3174       return false;
3175     }
3176   else if (dm->size1 == 0 || dm->size2 == 0)
3177     {
3178       msg (SE, _("Cannot index %zu×%zu matrix."),
3179            dm->size1, dm->size2);
3180       return false;
3181     }
3182   else if (lvalue->n_indexes == 1)
3183     {
3184       if (!is_vector (dm))
3185         {
3186           msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3187                dm->size1, dm->size2, lvalue->var->name);
3188           return false;
3189         }
3190       return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3191                                             MAX (dm->size1, dm->size2),
3192                                             IV_VECTOR, 0, iv0);
3193     }
3194   else
3195     {
3196       assert (lvalue->n_indexes == 2);
3197       if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3198                                           IV_ROW, dm->size2, iv0))
3199         return false;
3200
3201       if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3202                                           IV_COLUMN, dm->size1, iv1))
3203         {
3204           free (iv0->indexes);
3205           return false;
3206         }
3207       return true;
3208     }
3209 }
3210
3211 /* Takes ownership of M. */
3212 static bool
3213 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3214 {
3215   struct index_vector iv0, iv1;
3216   if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3217     {
3218       gsl_matrix_free (sm);
3219       return false;
3220     }
3221
3222   return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3223 }
3224
3225 \f
3226 /* Matrix command. */
3227
3228 struct matrix_cmds
3229   {
3230     struct matrix_cmd **commands;
3231     size_t n;
3232   };
3233
3234 static void matrix_cmds_uninit (struct matrix_cmds *);
3235
3236 struct matrix_cmd
3237   {
3238     enum matrix_cmd_type
3239       {
3240         MCMD_COMPUTE,
3241         MCMD_PRINT,
3242         MCMD_DO_IF,
3243         MCMD_LOOP,
3244         MCMD_BREAK,
3245         MCMD_DISPLAY,
3246         MCMD_RELEASE,
3247         MCMD_SAVE,
3248         MCMD_READ,
3249         MCMD_WRITE,
3250         MCMD_GET,
3251         MCMD_MSAVE,
3252         MCMD_MGET,
3253         MCMD_EIGEN,
3254         MCMD_SETDIAG,
3255         MCMD_SVD,
3256       }
3257     type;
3258
3259     union
3260       {
3261         struct compute_command
3262           {
3263             struct matrix_lvalue *lvalue;
3264             struct matrix_expr *rvalue;
3265           }
3266         compute;
3267
3268         struct print_command
3269           {
3270             struct matrix_expr *expression;
3271             bool use_default_format;
3272             struct fmt_spec format;
3273             char *title;
3274             int space;          /* -1 means NEWPAGE. */
3275
3276             struct print_labels
3277               {
3278                 struct string_array literals; /* CLABELS/RLABELS. */
3279                 struct matrix_expr *expr;     /* CNAMES/RNAMES. */
3280               }
3281             *rlabels, *clabels;
3282           }
3283         print;
3284
3285         struct do_if_command
3286           {
3287             struct do_if_clause
3288               {
3289                 struct matrix_expr *condition;
3290                 struct matrix_cmds commands;
3291               }
3292             *clauses;
3293             size_t n_clauses;
3294           }
3295         do_if;
3296
3297         struct loop_command
3298           {
3299             /* Index. */
3300             struct matrix_var *var;
3301             struct matrix_expr *start;
3302             struct matrix_expr *end;
3303             struct matrix_expr *increment;
3304
3305             /* Loop conditions. */
3306             struct matrix_expr *top_condition;
3307             struct matrix_expr *bottom_condition;
3308
3309             /* Commands. */
3310             struct matrix_cmds commands;
3311           }
3312         loop;
3313
3314         struct display_command
3315           {
3316             struct matrix_state *state;
3317             bool dictionary;
3318             bool status;
3319           }
3320         display;
3321
3322         struct release_command
3323           {
3324             struct matrix_var **vars;
3325             size_t n_vars;
3326           }
3327         release;
3328
3329         struct save_command
3330           {
3331             struct matrix_expr *expression;
3332             struct file_handle *outfile;
3333             struct string_array *variables;
3334             struct matrix_expr *names;
3335             struct stringi_set strings;
3336           }
3337         save;
3338
3339         struct read_command
3340           {
3341             struct read_file *rf;
3342             struct matrix_lvalue *dst;
3343             struct matrix_expr *size;
3344             int c1, c2;
3345             enum fmt_type format;
3346             int w;
3347             int d;
3348             bool symmetric;
3349             bool reread;
3350           }
3351         read;
3352
3353         struct write_command
3354           {
3355             struct matrix_expr *expression;
3356             struct file_handle *outfile;
3357             char *encoding;
3358             int c1, c2;
3359             enum fmt_type format;
3360             int w;
3361             int d;
3362             bool triangular;
3363             bool hold;          /* XXX */
3364           }
3365         write;
3366
3367         struct get_command
3368           {
3369             struct matrix_lvalue *dst;
3370             struct file_handle *file;
3371             char *encoding;
3372             struct string_array variables;
3373             struct matrix_var *names;
3374
3375             /* Treatment of missing values. */
3376             struct
3377               {
3378                 enum
3379                   {
3380                     MGET_ERROR,  /* Flag error on command. */
3381                     MGET_ACCEPT, /* Accept without change, user-missing only. */
3382                     MGET_OMIT,   /* Drop this case. */
3383                     MGET_RECODE  /* Recode to 'substitute'. */
3384                   }
3385                 treatment;
3386                 double substitute; /* MGET_RECODE only. */
3387               }
3388             user, system;
3389           }
3390         get;
3391
3392         struct msave_command
3393           {
3394             struct msave_common *common;
3395             char *varname_;
3396             struct matrix_expr *expr;
3397             const char *rowtype;
3398             struct matrix_expr *factors;
3399             struct matrix_expr *splits;
3400           }
3401          msave;
3402
3403         struct mget_command
3404           {
3405             struct matrix_state *state;
3406             struct file_handle *file;
3407             struct stringi_set rowtypes;
3408           }
3409         mget;
3410
3411         struct eigen_command
3412           {
3413             struct matrix_expr *expr;
3414             struct matrix_var *evec;
3415             struct matrix_var *eval;
3416           }
3417         eigen;
3418
3419         struct setdiag_command
3420           {
3421             struct matrix_var *dst;
3422             struct matrix_expr *expr;
3423           }
3424         setdiag;
3425
3426         struct svd_command
3427           {
3428             struct matrix_expr *expr;
3429             struct matrix_var *u;
3430             struct matrix_var *s;
3431             struct matrix_var *v;
3432           }
3433         svd;
3434       };
3435   };
3436
3437 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3438 static bool matrix_cmd_execute (struct matrix_cmd *);
3439 static void matrix_cmd_destroy (struct matrix_cmd *);
3440
3441 \f
3442 static struct matrix_cmd *
3443 matrix_parse_compute (struct matrix_state *s)
3444 {
3445   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3446   *cmd = (struct matrix_cmd) {
3447     .type = MCMD_COMPUTE,
3448     .compute = { .lvalue = NULL },
3449   };
3450
3451   struct compute_command *compute = &cmd->compute;
3452   compute->lvalue = matrix_lvalue_parse (s);
3453   if (!compute->lvalue)
3454     goto error;
3455
3456   if (!lex_force_match (s->lexer, T_EQUALS))
3457     goto error;
3458
3459   compute->rvalue = matrix_parse_expr (s);
3460   if (!compute->rvalue)
3461     goto error;
3462
3463   return cmd;
3464
3465 error:
3466   matrix_cmd_destroy (cmd);
3467   return NULL;
3468 }
3469
3470 static void
3471 matrix_cmd_execute_compute (struct compute_command *compute)
3472 {
3473   gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3474   if (!value)
3475     return;
3476
3477   matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3478 }
3479 \f
3480 static void
3481 print_labels_destroy (struct print_labels *labels)
3482 {
3483   if (labels)
3484     {
3485       string_array_destroy (&labels->literals);
3486       matrix_expr_destroy (labels->expr);
3487       free (labels);
3488     }
3489 }
3490
3491 static void
3492 parse_literal_print_labels (struct matrix_state *s,
3493                             struct print_labels **labelsp)
3494 {
3495   lex_match (s->lexer, T_EQUALS);
3496   print_labels_destroy (*labelsp);
3497   *labelsp = xzalloc (sizeof **labelsp);
3498   while (lex_token (s->lexer) != T_SLASH
3499          && lex_token (s->lexer) != T_ENDCMD
3500          && lex_token (s->lexer) != T_STOP)
3501     {
3502       struct string label = DS_EMPTY_INITIALIZER;
3503       while (lex_token (s->lexer) != T_COMMA
3504              && lex_token (s->lexer) != T_SLASH
3505              && lex_token (s->lexer) != T_ENDCMD
3506              && lex_token (s->lexer) != T_STOP)
3507         {
3508           if (!ds_is_empty (&label))
3509             ds_put_byte (&label, ' ');
3510
3511           if (lex_token (s->lexer) == T_STRING)
3512             ds_put_cstr (&label, lex_tokcstr (s->lexer));
3513           else
3514             {
3515               char *rep = lex_next_representation (s->lexer, 0, 0);
3516               ds_put_cstr (&label, rep);
3517               free (rep);
3518             }
3519           lex_get (s->lexer);
3520         }
3521       string_array_append_nocopy (&(*labelsp)->literals,
3522                                   ds_steal_cstr (&label));
3523
3524       if (!lex_match (s->lexer, T_COMMA))
3525         break;
3526     }
3527 }
3528
3529 static bool
3530 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3531 {
3532   lex_match (s->lexer, T_EQUALS);
3533   struct matrix_expr *e = matrix_parse_exp (s);
3534   if (!e)
3535     return false;
3536
3537   print_labels_destroy (*labelsp);
3538   *labelsp = xzalloc (sizeof **labelsp);
3539   (*labelsp)->expr = e;
3540   return true;
3541 }
3542
3543 static struct matrix_cmd *
3544 matrix_parse_print (struct matrix_state *s)
3545 {
3546   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3547   *cmd = (struct matrix_cmd) {
3548     .type = MCMD_PRINT,
3549     .print = {
3550       .use_default_format = true,
3551     }
3552   };
3553
3554   if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3555     {
3556       size_t depth = 0;
3557       for (size_t i = 0; ; i++)
3558         {
3559           enum token_type t = lex_next_token (s->lexer, i);
3560           if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3561             depth++;
3562           else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3563             depth--;
3564           else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3565             {
3566               if (i > 0)
3567                 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3568               break;
3569             }
3570         }
3571
3572       cmd->print.expression = matrix_parse_exp (s);
3573       if (!cmd->print.expression)
3574         goto error;
3575     }
3576
3577   while (lex_match (s->lexer, T_SLASH))
3578     {
3579       if (lex_match_id (s->lexer, "FORMAT"))
3580         {
3581           lex_match (s->lexer, T_EQUALS);
3582           if (!parse_format_specifier (s->lexer, &cmd->print.format))
3583             goto error;
3584           cmd->print.use_default_format = false;
3585         }
3586       else if (lex_match_id (s->lexer, "TITLE"))
3587         {
3588           lex_match (s->lexer, T_EQUALS);
3589           if (!lex_force_string (s->lexer))
3590             goto error;
3591           free (cmd->print.title);
3592           cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3593           lex_get (s->lexer);
3594         }
3595       else if (lex_match_id (s->lexer, "SPACE"))
3596         {
3597           lex_match (s->lexer, T_EQUALS);
3598           if (lex_match_id (s->lexer, "NEWPAGE"))
3599             cmd->print.space = -1;
3600           else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3601             {
3602               cmd->print.space = lex_integer (s->lexer);
3603               lex_get (s->lexer);
3604             }
3605           else
3606             goto error;
3607         }
3608       else if (lex_match_id (s->lexer, "RLABELS"))
3609         parse_literal_print_labels (s, &cmd->print.rlabels);
3610       else if (lex_match_id (s->lexer, "CLABELS"))
3611         parse_literal_print_labels (s, &cmd->print.clabels);
3612       else if (lex_match_id (s->lexer, "RNAMES"))
3613         {
3614           if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3615             goto error;
3616         }
3617       else if (lex_match_id (s->lexer, "CNAMES"))
3618         {
3619           if (!parse_expr_print_labels (s, &cmd->print.clabels))
3620             goto error;
3621         }
3622       else
3623         {
3624           lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3625                                "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3626           goto error;
3627         }
3628
3629     }
3630   return cmd;
3631
3632 error:
3633   matrix_cmd_destroy (cmd);
3634   return NULL;
3635 }
3636
3637 static bool
3638 matrix_is_integer (const gsl_matrix *m)
3639 {
3640   for (size_t y = 0; y < m->size1; y++)
3641     for (size_t x = 0; x < m->size2; x++)
3642       {
3643         double d = gsl_matrix_get (m, y, x);
3644         if (d != floor (d))
3645           return false;
3646       }
3647   return true;
3648 }
3649
3650 static double
3651 matrix_max_magnitude (const gsl_matrix *m)
3652 {
3653   double max = 0.0;
3654   for (size_t y = 0; y < m->size1; y++)
3655     for (size_t x = 0; x < m->size2; x++)
3656       {
3657         double d = fabs (gsl_matrix_get (m, y, x));
3658         if (d > max)
3659           max = d;
3660       }
3661   return max;
3662 }
3663
3664 static bool
3665 format_fits (struct fmt_spec format, double x)
3666 {
3667   char *s = data_out (&(union value) { .f = x }, NULL,
3668                       &format, settings_get_fmt_settings ());
3669   bool fits = *s != '*' && !strchr (s, 'E');
3670   free (s);
3671   return fits;
3672 }
3673
3674 static struct fmt_spec
3675 default_format (const gsl_matrix *m, int *log_scale)
3676 {
3677   *log_scale = 0;
3678
3679   double max = matrix_max_magnitude (m);
3680
3681   if (matrix_is_integer (m))
3682     for (int w = 1; w <= 12; w++)
3683       {
3684         struct fmt_spec format = { .type = FMT_F, .w = w };
3685         if (format_fits (format, -max))
3686           return format;
3687       }
3688
3689   if (max >= 1e9 || max <= 1e-4)
3690     {
3691       char s[64];
3692       snprintf (s, sizeof s, "%.1e", max);
3693
3694       const char *e = strchr (s, 'e');
3695       if (e)
3696         *log_scale = atoi (e + 1);
3697     }
3698
3699   return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3700 }
3701
3702 static char *
3703 trimmed_string (double d)
3704 {
3705   struct substring s = ss_buffer ((char *) &d, sizeof d);
3706   ss_rtrim (&s, ss_cstr (" "));
3707   return ss_xstrdup (s);
3708 }
3709
3710 static struct string_array *
3711 print_labels_get (const struct print_labels *labels, size_t n,
3712                   const char *prefix, bool truncate)
3713 {
3714   if (!labels)
3715     return NULL;
3716
3717   struct string_array *out = xzalloc (sizeof *out);
3718   if (labels->literals.n)
3719     string_array_clone (out, &labels->literals);
3720   else if (labels->expr)
3721     {
3722       gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3723       if (m && is_vector (m))
3724         {
3725           gsl_vector v = to_vector (m);
3726           for (size_t i = 0; i < v.size; i++)
3727             string_array_append_nocopy (out, trimmed_string (
3728                                           gsl_vector_get (&v, i)));
3729         }
3730       gsl_matrix_free (m);
3731     }
3732
3733   while (out->n < n)
3734     {
3735       if (labels->expr)
3736         string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3737                                                     out->n + 1));
3738       else
3739         string_array_append (out, "");
3740     }
3741   while (out->n > n)
3742     string_array_delete (out, out->n - 1);
3743
3744   if (truncate)
3745     for (size_t i = 0; i < out->n; i++)
3746       {
3747         char *s = out->strings[i];
3748         s[strnlen (s, 8)] = '\0';
3749       }
3750
3751   return out;
3752 }
3753
3754 static void
3755 matrix_cmd_print_space (int space)
3756 {
3757   if (space < 0)
3758     output_item_submit (page_break_item_create ());
3759   for (int i = 0; i < space; i++)
3760     output_log ("%s", "");
3761 }
3762
3763 static void
3764 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3765                        struct fmt_spec format, int log_scale)
3766 {
3767   matrix_cmd_print_space (print->space);
3768   if (print->title)
3769     output_log ("%s", print->title);
3770   if (log_scale != 0)
3771     output_log ("  10 ** %d   X", log_scale);
3772
3773   struct string_array *clabels = print_labels_get (print->clabels,
3774                                                    m->size2, "col", true);
3775   if (clabels && format.w < 8)
3776     format.w = 8;
3777   struct string_array *rlabels = print_labels_get (print->rlabels,
3778                                                    m->size1, "row", true);
3779
3780   if (clabels)
3781     {
3782       struct string line = DS_EMPTY_INITIALIZER;
3783       if (rlabels)
3784         ds_put_byte_multiple (&line, ' ', 8);
3785       for (size_t x = 0; x < m->size2; x++)
3786         ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3787       output_log_nocopy (ds_steal_cstr (&line));
3788     }
3789
3790   double scale = pow (10.0, log_scale);
3791   bool numeric = fmt_is_numeric (format.type);
3792   for (size_t y = 0; y < m->size1; y++)
3793     {
3794       struct string line = DS_EMPTY_INITIALIZER;
3795       if (rlabels)
3796         ds_put_format (&line, "%-8s", rlabels->strings[y]);
3797
3798       for (size_t x = 0; x < m->size2; x++)
3799         {
3800           double f = gsl_matrix_get (m, y, x);
3801           char *s = (numeric
3802                      ? data_out (&(union value) { .f = f / scale}, NULL,
3803                                  &format, settings_get_fmt_settings ())
3804                      : trimmed_string (f));
3805           ds_put_format (&line, " %s", s);
3806           free (s);
3807         }
3808       output_log_nocopy (ds_steal_cstr (&line));
3809     }
3810
3811   string_array_destroy (rlabels);
3812   free (rlabels);
3813   string_array_destroy (clabels);
3814   free (clabels);
3815 }
3816
3817 static void
3818 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3819                         const struct print_labels *print_labels, size_t n,
3820                         const char *name, const char *prefix)
3821 {
3822   struct string_array *labels = print_labels_get (print_labels, n, prefix,
3823                                                   false);
3824   struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3825   for (size_t i = 0; i < n; i++)
3826     pivot_category_create_leaf (
3827       d->root, (labels
3828                 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3829                 : pivot_value_new_integer (i + 1)));
3830   if (!labels)
3831     d->hide_all_labels = true;
3832   string_array_destroy (labels);
3833   free (labels);
3834 }
3835
3836 static void
3837 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3838                          struct fmt_spec format, int log_scale)
3839 {
3840   struct pivot_table *table = pivot_table_create__ (
3841     pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3842     "Matrix Print");
3843
3844   create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3845                           N_("Rows"), "row");
3846   create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3847                           N_("Columns"), "col");
3848
3849   struct pivot_footnote *footnote = NULL;
3850   if (log_scale != 0)
3851     {
3852       char *s = xasprintf ("× 10**%d", log_scale);
3853       footnote = pivot_table_create_footnote (
3854         table, pivot_value_new_user_text_nocopy (s));
3855     }
3856
3857   double scale = pow (10.0, log_scale);
3858   bool numeric = fmt_is_numeric (format.type);
3859   for (size_t y = 0; y < m->size1; y++)
3860     for (size_t x = 0; x < m->size2; x++)
3861       {
3862         double f = gsl_matrix_get (m, y, x);
3863         struct pivot_value *v;
3864         if (numeric)
3865           {
3866             v = pivot_value_new_number (f / scale);
3867             v->numeric.format = format;
3868           }
3869         else
3870           v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3871         if (footnote)
3872           pivot_value_add_footnote (v, footnote);
3873         pivot_table_put2 (table, y, x, v);
3874       }
3875
3876   pivot_table_submit (table);
3877 }
3878
3879 static void
3880 matrix_cmd_execute_print (const struct print_command *print)
3881 {
3882   if (print->expression)
3883     {
3884       gsl_matrix *m = matrix_expr_evaluate (print->expression);
3885       if (!m)
3886         return;
3887
3888       int log_scale = 0;
3889       struct fmt_spec format = (print->use_default_format
3890                                 ? default_format (m, &log_scale)
3891                                 : print->format);
3892
3893       if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3894         matrix_cmd_print_text (print, m, format, log_scale);
3895       else
3896         matrix_cmd_print_tables (print, m, format, log_scale);
3897
3898       gsl_matrix_free (m);
3899     }
3900   else
3901     {
3902       matrix_cmd_print_space (print->space);
3903       if (print->title)
3904         output_log ("%s", print->title);
3905     }
3906 }
3907 \f
3908 /* DO IF. */
3909
3910 static bool
3911 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3912                        const char *command_name,
3913                        const char *stop1, const char *stop2)
3914 {
3915   lex_end_of_command (s->lexer);
3916   lex_discard_rest_of_command (s->lexer);
3917
3918   size_t allocated = 0;
3919   for (;;)
3920     {
3921       while (lex_token (s->lexer) == T_ENDCMD)
3922         lex_get (s->lexer);
3923
3924       if (lex_at_phrase (s->lexer, stop1)
3925           || (stop2 && lex_at_phrase (s->lexer, stop2)))
3926         return true;
3927
3928       if (lex_at_phrase (s->lexer, "END MATRIX"))
3929         {
3930           msg (SE, _("Premature END MATRIX within %s."), command_name);
3931           return false;
3932         }
3933
3934       struct matrix_cmd *cmd = matrix_parse_command (s);
3935       if (!cmd)
3936         return false;
3937
3938       if (c->n >= allocated)
3939         c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
3940       c->commands[c->n++] = cmd;
3941     }
3942 }
3943
3944 static bool
3945 matrix_parse_do_if_clause (struct matrix_state *s,
3946                            struct do_if_command *ifc,
3947                            bool parse_condition,
3948                            size_t *allocated_clauses)
3949 {
3950   if (ifc->n_clauses >= *allocated_clauses)
3951     ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
3952                                sizeof *ifc->clauses);
3953   struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
3954   *c = (struct do_if_clause) { .condition = NULL };
3955
3956   if (parse_condition)
3957     {
3958       c->condition = matrix_parse_expr (s);
3959       if (!c->condition)
3960         return false;
3961     }
3962
3963   return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
3964 }
3965
3966 static struct matrix_cmd *
3967 matrix_parse_do_if (struct matrix_state *s)
3968 {
3969   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3970   *cmd = (struct matrix_cmd) {
3971     .type = MCMD_DO_IF,
3972     .do_if = { .n_clauses = 0 }
3973   };
3974
3975   size_t allocated_clauses = 0;
3976   do
3977     {
3978       if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
3979         goto error;
3980     }
3981   while (lex_match_phrase (s->lexer, "ELSE IF"));
3982
3983   if (lex_match_id (s->lexer, "ELSE")
3984       && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
3985     goto error;
3986
3987   if (!lex_match_phrase (s->lexer, "END IF"))
3988     NOT_REACHED ();
3989   return cmd;
3990
3991 error:
3992   matrix_cmd_destroy (cmd);
3993   return NULL;
3994 }
3995
3996 static bool
3997 matrix_cmd_execute_do_if (struct do_if_command *cmd)
3998 {
3999   for (size_t i = 0; i < cmd->n_clauses; i++)
4000     {
4001       struct do_if_clause *c = &cmd->clauses[i];
4002       if (c->condition)
4003         {
4004           double d;
4005           if (!matrix_expr_evaluate_scalar (c->condition,
4006                                             i ? "ELSE IF" : "DO IF",
4007                                             &d) || d <= 0)
4008             continue;
4009         }
4010
4011       for (size_t j = 0; j < c->commands.n; j++)
4012         if (!matrix_cmd_execute (c->commands.commands[j]))
4013           return false;
4014       return true;
4015     }
4016   return true;
4017 }
4018 \f
4019 static struct matrix_cmd *
4020 matrix_parse_loop (struct matrix_state *s)
4021 {
4022   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4023   *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4024
4025   struct loop_command *loop = &cmd->loop;
4026   if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4027     {
4028       struct substring name = lex_tokss (s->lexer);
4029       loop->var = matrix_var_lookup (s, name);
4030       if (!loop->var)
4031         loop->var = matrix_var_create (s, name);
4032
4033       lex_get (s->lexer);
4034       lex_get (s->lexer);
4035
4036       loop->start = matrix_parse_expr (s);
4037       if (!loop->start || !lex_force_match (s->lexer, T_TO))
4038         goto error;
4039
4040       loop->end = matrix_parse_expr (s);
4041       if (!loop->end)
4042         goto error;
4043
4044       if (lex_match (s->lexer, T_BY))
4045         {
4046           loop->increment = matrix_parse_expr (s);
4047           if (!loop->increment)
4048             goto error;
4049         }
4050     }
4051
4052   if (lex_match_id (s->lexer, "IF"))
4053     {
4054       loop->top_condition = matrix_parse_expr (s);
4055       if (!loop->top_condition)
4056         goto error;
4057     }
4058
4059   bool was_in_loop = s->in_loop;
4060   s->in_loop = true;
4061   bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4062                                    "END LOOP", NULL);
4063   s->in_loop = was_in_loop;
4064   if (!ok)
4065     goto error;
4066
4067   if (!lex_match_phrase (s->lexer, "END LOOP"))
4068     NOT_REACHED ();
4069
4070   if (lex_match_id (s->lexer, "IF"))
4071     {
4072       loop->bottom_condition = matrix_parse_expr (s);
4073       if (!loop->bottom_condition)
4074         goto error;
4075     }
4076
4077   return cmd;
4078
4079 error:
4080   matrix_cmd_destroy (cmd);
4081   return NULL;
4082 }
4083
4084 static void
4085 matrix_cmd_execute_loop (struct loop_command *cmd)
4086 {
4087   long int value = 0;
4088   long int end = 0;
4089   long int increment = 1;
4090   if (cmd->var)
4091     {
4092       if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4093           || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4094           || (cmd->increment
4095               && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4096                                                 &increment)))
4097         return;
4098
4099       if (increment > 0 ? value > end
4100           : increment < 0 ? value < end
4101           : true)
4102         return;
4103     }
4104
4105   int mxloops = settings_get_mxloops ();
4106   for (int i = 0; i < mxloops; i++)
4107     {
4108       if (cmd->var)
4109         {
4110           if (cmd->var->value
4111               && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4112             {
4113               gsl_matrix_free (cmd->var->value);
4114               cmd->var->value = NULL;
4115             }
4116           if (!cmd->var->value)
4117             cmd->var->value = gsl_matrix_alloc (1, 1);
4118
4119           gsl_matrix_set (cmd->var->value, 0, 0, value);
4120         }
4121
4122       if (cmd->top_condition)
4123         {
4124           double d;
4125           if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4126                                             &d) || d <= 0)
4127             return;
4128         }
4129
4130       for (size_t j = 0; j < cmd->commands.n; j++)
4131         if (!matrix_cmd_execute (cmd->commands.commands[j]))
4132           return;
4133
4134       if (cmd->bottom_condition)
4135         {
4136           double d;
4137           if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4138                                             "END LOOP IF", &d) || d > 0)
4139             return;
4140         }
4141
4142       if (cmd->var)
4143         {
4144           value += increment;
4145           if (increment > 0 ? value > end : value < end)
4146             return;
4147         }
4148     }
4149 }
4150 \f
4151 static struct matrix_cmd *
4152 matrix_parse_break (struct matrix_state *s)
4153 {
4154   if (!s->in_loop)
4155     {
4156       msg (SE, _("BREAK not inside LOOP."));
4157       return NULL;
4158     }
4159
4160   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4161   *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4162   return cmd;
4163 }
4164 \f
4165 static struct matrix_cmd *
4166 matrix_parse_display (struct matrix_state *s)
4167 {
4168   bool dictionary = false;
4169   bool status = false;
4170   while (lex_token (s->lexer) == T_ID)
4171     {
4172       if (lex_match_id (s->lexer, "DICTIONARY"))
4173         dictionary = true;
4174       else if (lex_match_id (s->lexer, "STATUS"))
4175         status = true;
4176       else
4177         {
4178           lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4179           return NULL;
4180         }
4181     }
4182   if (!dictionary && !status)
4183     dictionary = status = true;
4184
4185   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4186   *cmd = (struct matrix_cmd) {
4187     .type = MCMD_DISPLAY,
4188     .display = {
4189       .state = s,
4190       .dictionary = dictionary,
4191       .status = status,
4192     }
4193   };
4194   return cmd;
4195 }
4196
4197 static int
4198 compare_matrix_var_pointers (const void *a_, const void *b_)
4199 {
4200   const struct matrix_var *const *ap = a_;
4201   const struct matrix_var *const *bp = b_;
4202   const struct matrix_var *a = *ap;
4203   const struct matrix_var *b = *bp;
4204   return strcmp (a->name, b->name);
4205 }
4206
4207 static void
4208 matrix_cmd_execute_display (struct display_command *cmd)
4209 {
4210   const struct matrix_state *s = cmd->state;
4211
4212   struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4213   pivot_dimension_create (
4214     table, PIVOT_AXIS_COLUMN, N_("Property"),
4215     N_("Rows"), N_("Columns"), N_("Size (kB)"));
4216
4217   const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4218   size_t n_vars = 0;
4219
4220   const struct matrix_var *var;
4221   HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4222     if (var->value)
4223       vars[n_vars++] = var;
4224   qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4225
4226   struct pivot_dimension *rows = pivot_dimension_create (
4227     table, PIVOT_AXIS_ROW, N_("Variable"));
4228   for (size_t i = 0; i < n_vars; i++)
4229     {
4230       const struct matrix_var *var = vars[i];
4231       pivot_category_create_leaf (
4232         rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4233
4234       size_t r = var->value->size1;
4235       size_t c = var->value->size2;
4236       double values[] = { r, c, r * c * sizeof (double) / 1024 };
4237       for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4238         pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4239     }
4240   pivot_table_submit (table);
4241 }
4242 \f
4243 static struct matrix_cmd *
4244 matrix_parse_release (struct matrix_state *s)
4245 {
4246   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4247   *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4248
4249   size_t allocated_vars = 0;
4250   while (lex_token (s->lexer) == T_ID)
4251     {
4252       struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4253       if (var)
4254         {
4255           if (cmd->release.n_vars >= allocated_vars)
4256             cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4257                                             sizeof *cmd->release.vars);
4258           cmd->release.vars[cmd->release.n_vars++] = var;
4259         }
4260       else
4261         lex_error (s->lexer, _("Variable name expected."));
4262       lex_get (s->lexer);
4263
4264       if (!lex_match (s->lexer, T_COMMA))
4265         break;
4266     }
4267
4268   return cmd;
4269 }
4270
4271 static void
4272 matrix_cmd_execute_release (struct release_command *cmd)
4273 {
4274   for (size_t i = 0; i < cmd->n_vars; i++)
4275     {
4276       struct matrix_var *var = cmd->vars[i];
4277       gsl_matrix_free (var->value);
4278       var->value = NULL;
4279     }
4280 }
4281 \f
4282 static struct matrix_cmd *
4283 matrix_parse_save (struct matrix_state *s)
4284 {
4285   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4286   *cmd = (struct matrix_cmd) {
4287     .type = MCMD_SAVE,
4288     .save = {
4289       .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4290     }
4291   };
4292
4293   struct save_command *save = &cmd->save;
4294   save->expression = matrix_parse_exp (s);
4295   if (!save->expression)
4296     goto error;
4297
4298   while (lex_match (s->lexer, T_SLASH))
4299     {
4300       if (lex_match_id (s->lexer, "OUTFILE"))
4301         {
4302           lex_match (s->lexer, T_EQUALS);
4303
4304           fh_unref (save->outfile);
4305           save->outfile = (lex_match (s->lexer, T_ASTERISK)
4306                            ? fh_inline_file ()
4307                            : fh_parse (s->lexer, FH_REF_FILE, s->session));
4308           if (!save->outfile)
4309             goto error;
4310         }
4311       else if (lex_match_id (s->lexer, "VARIABLES"))
4312         {
4313           lex_match (s->lexer, T_EQUALS);
4314
4315           char **names;
4316           size_t n;
4317           struct dictionary *d = dict_create (get_default_encoding ());
4318           bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4319                                           PV_NO_SCRATCH | PV_NO_DUPLICATE);
4320           dict_unref (d);
4321           if (!ok)
4322             goto error;
4323
4324           string_array_destroy (save->variables);
4325           if (!save->variables)
4326             save->variables = xmalloc (sizeof *save->variables);
4327           *save->variables = (struct string_array) {
4328             .strings = names,
4329             .n = n,
4330             .allocated = n,
4331           };
4332         }
4333       else if (lex_match_id (s->lexer, "NAMES"))
4334         {
4335           lex_match (s->lexer, T_EQUALS);
4336           matrix_expr_destroy (save->names);
4337           save->names = matrix_parse_exp (s);
4338           if (!save->names)
4339             goto error;
4340         }
4341       else if (lex_match_id (s->lexer, "STRINGS"))
4342         {
4343           lex_match (s->lexer, T_EQUALS);
4344           while (lex_token (s->lexer) == T_ID)
4345             {
4346               stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4347               lex_get (s->lexer);
4348               if (!lex_match (s->lexer, T_COMMA))
4349                 break;
4350             }
4351         }
4352       else
4353         {
4354           lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4355                                "STRINGS");
4356           goto error;
4357         }
4358     }
4359
4360   if (!save->outfile)
4361     {
4362       if (s->prev_save_outfile)
4363         save->outfile = fh_ref (s->prev_save_outfile);
4364       else
4365         {
4366           lex_sbc_missing ("OUTFILE");
4367           goto error;
4368         }
4369     }
4370   fh_unref (s->prev_save_outfile);
4371   s->prev_save_outfile = fh_ref (save->outfile);
4372
4373   if (save->variables && save->names)
4374     {
4375       msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4376       matrix_expr_destroy (save->names);
4377       save->names = NULL;
4378     }
4379
4380   return cmd;
4381
4382 error:
4383   matrix_cmd_destroy (cmd);
4384   return NULL;
4385 }
4386
4387 static void
4388 matrix_cmd_execute_save (const struct save_command *save)
4389 {
4390   assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4391
4392   gsl_matrix *m = matrix_expr_evaluate (save->expression);
4393   if (!m)
4394     return;
4395
4396   bool ok = true;
4397   struct dictionary *dict = dict_create (get_default_encoding ());
4398   struct string_array names = { .n = 0 };
4399   if (save->variables)
4400     string_array_clone (&names, save->variables);
4401   else if (save->names)
4402     {
4403       gsl_matrix *nm = matrix_expr_evaluate (save->names);
4404       if (nm && is_vector (nm))
4405         {
4406           gsl_vector nv = to_vector (nm);
4407           for (size_t i = 0; i < nv.size; i++)
4408             {
4409               char *name = trimmed_string (gsl_vector_get (&nv, i));
4410               if (dict_id_is_valid (dict, name, true))
4411                 string_array_append_nocopy (&names, name);
4412               else
4413                 ok = false;
4414             }
4415         }
4416     }
4417
4418   struct stringi_set strings;
4419   stringi_set_clone (&strings, &save->strings);
4420
4421   for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4422     {
4423       char tmp_name[64];
4424       const char *name;
4425       if (i < names.n)
4426         name = names.strings[i];
4427       else
4428         {
4429           snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4430           name = tmp_name;
4431         }
4432
4433       int width = stringi_set_delete (&strings, name) ? 8 : 0;
4434       struct variable *var = dict_create_var (dict, name, width);
4435       if (!var)
4436         {
4437           msg (SE, _("Duplicate variable name %s in SAVE statement."),
4438                name);
4439           ok = false;
4440         }
4441     }
4442
4443   if (!stringi_set_is_empty (&strings))
4444     {
4445       const char *example = stringi_set_node_get_string (
4446         stringi_set_first (&strings));
4447       msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4448                          "with that name was found on SAVE.",
4449                          "STRINGS specified %2$zu variables, including %1$s, "
4450                          "whose names were not found on SAVE.",
4451                          stringi_set_count (&strings)),
4452            example, stringi_set_count (&strings));
4453       ok = false;
4454     }
4455   stringi_set_destroy (&strings);
4456   string_array_destroy (&names);
4457
4458   if (!ok)
4459     {
4460       gsl_matrix_free (m);
4461       dict_unref (dict);
4462       return;
4463     }
4464
4465   struct casewriter *writer = any_writer_open (save->outfile, dict);
4466   if (!writer)
4467     {
4468       gsl_matrix_free (m);
4469       dict_unref (dict);
4470       return;
4471     }
4472
4473   const struct caseproto *proto = dict_get_proto (dict);
4474   for (size_t y = 0; y < m->size1; y++)
4475     {
4476       struct ccase *c = case_create (proto);
4477       for (size_t x = 0; x < m->size2; x++)
4478         {
4479           double d = gsl_matrix_get (m, y, x);
4480           int width = caseproto_get_width (proto, x);
4481           union value *value = case_data_rw_idx (c, x);
4482           if (width == 0)
4483             value->f = d;
4484           else
4485             memcpy (value->s, &d, width);
4486         }
4487       casewriter_write (writer, c);
4488     }
4489   casewriter_destroy (writer);
4490
4491   gsl_matrix_free (m);
4492   dict_unref (dict);
4493 }
4494 \f
4495 static struct matrix_cmd *
4496 matrix_parse_read (struct matrix_state *s)
4497 {
4498   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4499   *cmd = (struct matrix_cmd) {
4500     .type = MCMD_READ,
4501     .read = { .format = FMT_F },
4502   };
4503
4504   struct file_handle *fh = NULL;
4505   char *encoding = NULL;
4506   struct read_command *read = &cmd->read;
4507   read->dst = matrix_lvalue_parse (s);
4508   if (!read->dst)
4509     goto error;
4510
4511   int by = 0;
4512   int repetitions = 0;
4513   int record_width = 0;
4514   bool seen_format = false;
4515   while (lex_match (s->lexer, T_SLASH))
4516     {
4517       if (lex_match_id (s->lexer, "FILE"))
4518         {
4519           lex_match (s->lexer, T_EQUALS);
4520
4521           fh_unref (fh);
4522           fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4523           if (!fh)
4524             goto error;
4525         }
4526       else if (lex_match_id (s->lexer, "ENCODING"))
4527         {
4528           lex_match (s->lexer, T_EQUALS);
4529           if (!lex_force_string (s->lexer))
4530             goto error;
4531
4532           free (encoding);
4533           encoding = ss_xstrdup (lex_tokss (s->lexer));
4534
4535           lex_get (s->lexer);
4536         }
4537       else if (lex_match_id (s->lexer, "FIELD"))
4538         {
4539           lex_match (s->lexer, T_EQUALS);
4540
4541           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4542             goto error;
4543           read->c1 = lex_integer (s->lexer);
4544           lex_get (s->lexer);
4545           if (!lex_force_match (s->lexer, T_TO)
4546               || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4547             goto error;
4548           read->c2 = lex_integer (s->lexer) + 1;
4549           lex_get (s->lexer);
4550
4551           record_width = read->c2 - read->c1;
4552           if (lex_match (s->lexer, T_BY))
4553             {
4554               if (!lex_force_int_range (s->lexer, "BY", 1,
4555                                         read->c2 - read->c1))
4556                 goto error;
4557               by = lex_integer (s->lexer);
4558               lex_get (s->lexer);
4559
4560               if (record_width % by)
4561                 {
4562                   msg (SE, _("BY %d does not evenly divide record width %d."),
4563                        by, record_width);
4564                   goto error;
4565                 }
4566             }
4567           else
4568             by = 0;
4569         }
4570       else if (lex_match_id (s->lexer, "SIZE"))
4571         {
4572           lex_match (s->lexer, T_EQUALS);
4573           read->size = matrix_parse_exp (s);
4574           if (!read->size)
4575             goto error;
4576         }
4577       else if (lex_match_id (s->lexer, "MODE"))
4578         {
4579           lex_match (s->lexer, T_EQUALS);
4580           if (lex_match_id (s->lexer, "RECTANGULAR"))
4581             read->symmetric = false;
4582           else if (lex_match_id (s->lexer, "SYMMETRIC"))
4583             read->symmetric = true;
4584           else
4585             {
4586               lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4587               goto error;
4588             }
4589         }
4590       else if (lex_match_id (s->lexer, "REREAD"))
4591         read->reread = true;
4592       else if (lex_match_id (s->lexer, "FORMAT"))
4593         {
4594           if (seen_format)
4595             {
4596               lex_sbc_only_once ("FORMAT");
4597               goto error;
4598             }
4599           seen_format = true;
4600
4601           lex_match (s->lexer, T_EQUALS);
4602
4603           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4604             goto error;
4605
4606           const char *p = lex_tokcstr (s->lexer);
4607           if (c_isdigit (p[0]))
4608             {
4609               repetitions = atoi (p);
4610               p += strspn (p, "0123456789");
4611               if (!fmt_from_name (p, &read->format))
4612                 {
4613                   lex_error (s->lexer, _("Unknown format %s."), p);
4614                   goto error;
4615                 }
4616               lex_get (s->lexer);
4617             }
4618           else if (!fmt_from_name (p, &read->format))
4619             {
4620               struct fmt_spec format;
4621               if (!parse_format_specifier (s->lexer, &format))
4622                 goto error;
4623               read->format = format.type;
4624               read->w = format.w;
4625               read->d = format.d;
4626             }
4627         }
4628       else
4629         {
4630           lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4631                                "REREAD", "FORMAT");
4632           goto error;
4633         }
4634     }
4635
4636   if (!read->c1)
4637     {
4638       lex_sbc_missing ("FIELD");
4639       goto error;
4640     }
4641
4642   if (!read->dst->n_indexes && !read->size)
4643     {
4644       msg (SE, _("SIZE is required for reading data into a full matrix "
4645                  "(as opposed to a submatrix)."));
4646       goto error;
4647     }
4648
4649   if (!fh)
4650     {
4651       if (s->prev_read_file)
4652         fh = fh_ref (s->prev_read_file);
4653       else
4654         {
4655           lex_sbc_missing ("FILE");
4656           goto error;
4657         }
4658     }
4659   fh_unref (s->prev_read_file);
4660   s->prev_read_file = fh_ref (fh);
4661
4662   read->rf = read_file_create (s, fh);
4663   if (encoding)
4664     {
4665       free (read->rf->encoding);
4666       read->rf->encoding = encoding;
4667       encoding = NULL;
4668     }
4669
4670   /* Field width may be specified in multiple ways:
4671
4672      1. BY on FIELD.
4673      2. The format on FORMAT.
4674      3. The repetition factor on FORMAT.
4675
4676      (2) and (3) are mutually exclusive.
4677
4678      If more than one of these is present, they must agree.  If none of them is
4679      present, then free-field format is used.
4680    */
4681   if (repetitions > record_width)
4682     {
4683       msg (SE, _("%d repetitions cannot fit in record width %d."),
4684            repetitions, record_width);
4685       goto error;
4686     }
4687   int w = (repetitions ? record_width / repetitions
4688            : read->w ? read->w
4689            : by);
4690   if (by && w != by)
4691     {
4692       if (repetitions)
4693         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4694                    "which implies field width %d, "
4695                    "but BY specifies field width %d."),
4696              repetitions, record_width, w, by);
4697       else
4698         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4699              w, by);
4700       goto error;
4701     }
4702   read->w = w;
4703   return cmd;
4704
4705 error:
4706   matrix_cmd_destroy (cmd);
4707   free (encoding);
4708   return NULL;
4709 }
4710
4711 static void
4712 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4713                        gsl_matrix *m, struct substring p, size_t y, size_t x)
4714 {
4715   const char *input_encoding = dfm_reader_get_encoding (reader);
4716   union value v;
4717   char *error = data_in (p, input_encoding, read->format,
4718                          settings_get_fmt_settings (), &v, 0, NULL);
4719   /* XXX report error if value is missing */
4720   if (error)
4721     msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4722   else
4723     {
4724       gsl_matrix_set (m, y, x, v.f);
4725       if (read->symmetric && x != y)
4726         gsl_matrix_set (m, x, y, v.f);
4727     }
4728 }
4729
4730 static bool
4731 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4732                   struct substring *line)
4733 {
4734   if (dfm_eof (reader))
4735     {
4736       msg (SE, _("Unexpected end of file reading matrix data."));
4737       return false;
4738     }
4739   dfm_expand_tabs (reader);
4740   *line = ss_substr (dfm_get_record (reader),
4741                      read->c1 - 1, read->c2 - read->c1);
4742   return true;
4743 }
4744
4745 static void
4746 matrix_read (struct read_command *read, struct dfm_reader *reader,
4747              gsl_matrix *m)
4748 {
4749   for (size_t y = 0; y < m->size1; y++)
4750     {
4751       size_t nx = read->symmetric ? y + 1 : m->size2;
4752
4753       struct substring line = ss_empty ();
4754       for (size_t x = 0; x < nx; x++)
4755         {
4756           struct substring p;
4757           if (!read->w)
4758             {
4759               for (;;)
4760                 {
4761                   ss_ltrim (&line, ss_cstr (" ,"));
4762                   if (!ss_is_empty (line))
4763                     break;
4764                   if (!matrix_read_line (read, reader, &line))
4765                     return;
4766                   dfm_forward_record (reader);
4767                 }
4768
4769               ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4770             }
4771           else
4772             {
4773               if (!matrix_read_line (read, reader, &line))
4774                 return;
4775               size_t fields_per_line = (read->c2 - read->c1) / read->w;
4776               int f = x % fields_per_line;
4777               if (f == fields_per_line - 1)
4778                 dfm_forward_record (reader);
4779
4780               p = ss_substr (line, read->w * f, read->w);
4781             }
4782
4783           matrix_read_set_field (read, reader, m, p, y, x);
4784         }
4785
4786       if (read->w)
4787         dfm_forward_record (reader);
4788       else
4789         {
4790           ss_ltrim (&line, ss_cstr (" ,"));
4791           if (!ss_is_empty (line))
4792             msg (SW, _("Trailing garbage on line \"%.*s\""),
4793                  (int) line.length, line.string);
4794         }
4795     }
4796 }
4797
4798 static void
4799 matrix_cmd_execute_read (struct read_command *read)
4800 {
4801   struct index_vector iv0, iv1;
4802   if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4803     return;
4804
4805   size_t size[2] = { SIZE_MAX, SIZE_MAX };
4806   if (read->size)
4807     {
4808       gsl_matrix *m = matrix_expr_evaluate (read->size);
4809       if (!m)
4810         return;
4811
4812       if (!is_vector (m))
4813         {
4814           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4815                      "not a %zu×%zu matrix."), m->size1, m->size2);
4816           gsl_matrix_free (m);
4817           free (iv0.indexes);
4818           free (iv1.indexes);
4819           return;
4820         }
4821
4822       gsl_vector v = to_vector (m);
4823       double d[2];
4824       if (v.size == 1)
4825         {
4826           d[0] = gsl_vector_get (&v, 0);
4827           d[1] = 1;
4828         }
4829       else if (v.size == 2)
4830         {
4831           d[0] = gsl_vector_get (&v, 0);
4832           d[1] = gsl_vector_get (&v, 1);
4833         }
4834       else
4835         {
4836           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4837                      "not a %zu×%zu matrix."),
4838                m->size1, m->size2),
4839           gsl_matrix_free (m);
4840           free (iv0.indexes);
4841           free (iv1.indexes);
4842           return;
4843         }
4844       gsl_matrix_free (m);
4845
4846       if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4847         {
4848           msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4849           free (iv0.indexes);
4850           free (iv1.indexes);
4851           return;
4852         }
4853
4854       size[0] = d[0];
4855       size[1] = d[1];
4856     }
4857
4858   if (read->dst->n_indexes)
4859     {
4860       size_t submatrix_size[2];
4861       if (read->dst->n_indexes == 2)
4862         {
4863           submatrix_size[0] = iv0.n;
4864           submatrix_size[1] = iv1.n;
4865         }
4866       else if (read->dst->var->value->size1 == 1)
4867         {
4868           submatrix_size[0] = 1;
4869           submatrix_size[1] = iv0.n;
4870         }
4871       else
4872         {
4873           submatrix_size[0] = iv0.n;
4874           submatrix_size[1] = 1;
4875         }
4876
4877       if (read->size)
4878         {
4879           if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4880             {
4881               msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4882                          "%zu×%zu."),
4883                    size[0], size[1],
4884                    submatrix_size[0], submatrix_size[1]);
4885               free (iv0.indexes);
4886               free (iv1.indexes);
4887               return;
4888             }
4889         }
4890       else
4891         {
4892           size[0] = submatrix_size[0];
4893           size[1] = submatrix_size[1];
4894         }
4895     }
4896
4897   struct dfm_reader *reader = read_file_open (read->rf);
4898   if (read->reread)
4899     dfm_reread_record (reader, 1);
4900
4901   if (read->symmetric && size[0] != size[1])
4902     {
4903       msg (SE, _("Cannot read non-square %zu×%zu matrix "
4904                  "using READ with MODE=SYMMETRIC."),
4905            size[0], size[1]);
4906       free (iv0.indexes);
4907       free (iv1.indexes);
4908       return;
4909     }
4910   gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4911   matrix_read (read, reader, tmp);
4912   matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4913 }
4914 \f
4915 static struct matrix_cmd *
4916 matrix_parse_write (struct matrix_state *s)
4917 {
4918   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4919   *cmd = (struct matrix_cmd) {
4920     .type = MCMD_WRITE,
4921     .write = { .format = FMT_F },
4922   };
4923
4924   struct write_command *write = &cmd->write;
4925   write->expression = matrix_parse_exp (s);
4926   if (!write->expression)
4927     goto error;
4928
4929   int by = 0;
4930   int repetitions = 0;
4931   int record_width = 0;
4932   bool seen_format = false;
4933   while (lex_match (s->lexer, T_SLASH))
4934     {
4935       if (lex_match_id (s->lexer, "OUTFILE"))
4936         {
4937           lex_match (s->lexer, T_EQUALS);
4938
4939           fh_unref (write->outfile);
4940           write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
4941           if (!write->outfile)
4942             goto error;
4943         }
4944       else if (lex_match_id (s->lexer, "ENCODING"))
4945         {
4946           lex_match (s->lexer, T_EQUALS);
4947           if (!lex_force_string (s->lexer))
4948             goto error;
4949
4950           free (write->encoding);
4951           write->encoding = ss_xstrdup (lex_tokss (s->lexer));
4952
4953           lex_get (s->lexer);
4954         }
4955       else if (lex_match_id (s->lexer, "FIELD"))
4956         {
4957           lex_match (s->lexer, T_EQUALS);
4958
4959           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4960             goto error;
4961           write->c1 = lex_integer (s->lexer);
4962           lex_get (s->lexer);
4963           if (!lex_force_match (s->lexer, T_TO)
4964               || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
4965             goto error;
4966           write->c2 = lex_integer (s->lexer) + 1;
4967           lex_get (s->lexer);
4968
4969           record_width = write->c2 - write->c1;
4970           if (lex_match (s->lexer, T_BY))
4971             {
4972               if (!lex_force_int_range (s->lexer, "BY", 1,
4973                                         write->c2 - write->c1))
4974                 goto error;
4975               by = lex_integer (s->lexer);
4976               lex_get (s->lexer);
4977
4978               if (record_width % by)
4979                 {
4980                   msg (SE, _("BY %d does not evenly divide record width %d."),
4981                        by, record_width);
4982                   goto error;
4983                 }
4984             }
4985           else
4986             by = 0;
4987         }
4988       else if (lex_match_id (s->lexer, "MODE"))
4989         {
4990           lex_match (s->lexer, T_EQUALS);
4991           if (lex_match_id (s->lexer, "RECTANGULAR"))
4992             write->triangular = false;
4993           else if (lex_match_id (s->lexer, "TRIANGULAR"))
4994             write->triangular = true;
4995           else
4996             {
4997               lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
4998               goto error;
4999             }
5000         }
5001       else if (lex_match_id (s->lexer, "HOLD"))
5002         write->hold = true;
5003       else if (lex_match_id (s->lexer, "FORMAT"))
5004         {
5005           if (seen_format)
5006             {
5007               lex_sbc_only_once ("FORMAT");
5008               goto error;
5009             }
5010           seen_format = true;
5011
5012           lex_match (s->lexer, T_EQUALS);
5013
5014           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5015             goto error;
5016
5017           const char *p = lex_tokcstr (s->lexer);
5018           if (c_isdigit (p[0]))
5019             {
5020               repetitions = atoi (p);
5021               p += strspn (p, "0123456789");
5022               if (!fmt_from_name (p, &write->format))
5023                 {
5024                   lex_error (s->lexer, _("Unknown format %s."), p);
5025                   goto error;
5026                 }
5027               lex_get (s->lexer);
5028             }
5029           else if (!fmt_from_name (p, &write->format))
5030             {
5031               struct fmt_spec format;
5032               if (!parse_format_specifier (s->lexer, &format))
5033                 goto error;
5034               write->format = format.type;
5035               write->w = format.w;
5036               write->d = format.d;
5037             }
5038         }
5039       else
5040         {
5041           lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5042                                "HOLD", "FORMAT");
5043           goto error;
5044         }
5045     }
5046
5047   if (!write->c1)
5048     {
5049       lex_sbc_missing ("FIELD");
5050       goto error;
5051     }
5052
5053   if (!write->outfile)
5054     {
5055       if (s->prev_write_outfile)
5056         write->outfile = fh_ref (s->prev_write_outfile);
5057       else
5058         {
5059           lex_sbc_missing ("OUTFILE");
5060           goto error;
5061         }
5062     }
5063   fh_unref (s->prev_write_outfile);
5064   s->prev_write_outfile = fh_ref (write->outfile);
5065
5066   /* Field width may be specified in multiple ways:
5067
5068      1. BY on FIELD.
5069      2. The format on FORMAT.
5070      3. The repetition factor on FORMAT.
5071
5072      (2) and (3) are mutually exclusive.
5073
5074      If more than one of these is present, they must agree.  If none of them is
5075      present, then free-field format is used.
5076    */
5077   if (repetitions > record_width)
5078     {
5079       msg (SE, _("%d repetitions cannot fit in record width %d."),
5080            repetitions, record_width);
5081       goto error;
5082     }
5083   int w = (repetitions ? record_width / repetitions
5084            : write->w ? write->w
5085            : by);
5086   if (by && w != by)
5087     {
5088       if (repetitions)
5089         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5090                    "which implies field width %d, "
5091                    "but BY specifies field width %d."),
5092              repetitions, record_width, w, by);
5093       else
5094         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5095              w, by);
5096       goto error;
5097     }
5098   write->w = w;
5099   return cmd;
5100
5101 error:
5102   matrix_cmd_destroy (cmd);
5103   return NULL;
5104 }
5105
5106 static void
5107 matrix_cmd_execute_write (struct write_command *write)
5108 {
5109   gsl_matrix *m = matrix_expr_evaluate (write->expression);
5110   if (!m)
5111     return;
5112
5113   if (write->triangular && m->size1 != m->size2)
5114     {
5115       msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5116                  "the matrix to be written has dimensions %zu×%zu."),
5117            m->size1, m->size2);
5118       gsl_matrix_free (m);
5119       return;
5120     }
5121
5122   struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
5123   if (!writer)
5124     {
5125       gsl_matrix_free (m);
5126       return;
5127     }
5128
5129   const struct fmt_settings *settings = settings_get_fmt_settings ();
5130   struct fmt_spec format = {
5131     .type = write->format,
5132     .w = write->w ? write->w : 40,
5133     .d = write->d
5134   };
5135   struct u8_line line = U8_LINE_INITIALIZER;
5136   for (size_t y = 0; y < m->size1; y++)
5137     {
5138       size_t nx = write->triangular ? y + 1 : m->size2;
5139       int x0 = write->c1;
5140       for (size_t x = 0; x < nx; x++)
5141         {
5142           /* XXX string values */
5143           union value v = { .f = gsl_matrix_get (m, y, x) };
5144           char *s = (write->w
5145                      ? data_out (&v, NULL, &format, settings)
5146                      : data_out_stretchy (&v, NULL, &format, settings, NULL));
5147           size_t len = strlen (s);
5148           int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5149           if (width + x0 > write->c2)
5150             {
5151               dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5152               u8_line_clear (&line);
5153               x0 = write->c1;
5154             }
5155           u8_line_put (&line, x0, x0 + width, s, len);
5156           free (s);
5157
5158           x0 += write->w ? write->w : width + 1;
5159         }
5160
5161       dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5162       u8_line_clear (&line);
5163     }
5164   u8_line_destroy (&line);
5165   dfm_close_writer (writer);
5166
5167   gsl_matrix_free (m);
5168 }
5169 \f
5170 static struct matrix_cmd *
5171 matrix_parse_get (struct matrix_state *s)
5172 {
5173   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5174   *cmd = (struct matrix_cmd) {
5175     .type = MCMD_GET,
5176     .get = {
5177       .user = { .treatment = MGET_ERROR },
5178       .system = { .treatment = MGET_ERROR },
5179     }
5180   };
5181
5182   struct get_command *get = &cmd->get;
5183   get->dst = matrix_lvalue_parse (s);
5184   if (!get->dst)
5185     goto error;
5186
5187   while (lex_match (s->lexer, T_SLASH))
5188     {
5189       if (lex_match_id (s->lexer, "FILE"))
5190         {
5191           if (get->variables.n)
5192             {
5193               lex_error (s->lexer, _("FILE must precede VARIABLES"));
5194               goto error;
5195             }
5196           lex_match (s->lexer, T_EQUALS);
5197
5198           fh_unref (get->file);
5199           get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5200           if (!get->file)
5201             goto error;
5202         }
5203       else if (lex_match_id (s->lexer, "ENCODING"))
5204         {
5205           if (get->variables.n)
5206             {
5207               lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5208               goto error;
5209             }
5210           lex_match (s->lexer, T_EQUALS);
5211           if (!lex_force_string (s->lexer))
5212             goto error;
5213
5214           free (get->encoding);
5215           get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5216
5217           lex_get (s->lexer);
5218         }
5219       else if (lex_match_id (s->lexer, "VARIABLES"))
5220         {
5221           lex_match (s->lexer, T_EQUALS);
5222
5223           struct dictionary *dict = NULL;
5224           if (!get->file)
5225             {
5226               dict = dataset_dict (s->dataset);
5227               if (dict_get_var_cnt (dict) == 0)
5228                 {
5229                   lex_error (s->lexer, _("GET cannot read empty active file."));
5230                   goto error;
5231                 }
5232             }
5233           else
5234             {
5235               struct casereader *reader = any_reader_open_and_decode (
5236                 get->file, get->encoding, &dict, NULL);
5237               if (!reader)
5238                 goto error;
5239               casereader_destroy (reader);
5240             }
5241
5242           struct variable **vars;
5243           size_t n_vars;
5244           bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5245                                      PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5246           if (!ok)
5247             {
5248               dict_unref (dict);
5249               goto error;
5250             }
5251
5252           string_array_clear (&get->variables);
5253           for (size_t i = 0; i < n_vars; i++)
5254             string_array_append (&get->variables, var_get_name (vars[i]));
5255           free (vars);
5256           dict_unref (dict);
5257         }
5258       else if (lex_match_id (s->lexer, "NAMES"))
5259         {
5260           lex_match (s->lexer, T_EQUALS);
5261           if (!lex_force_id (s->lexer))
5262             goto error;
5263
5264           struct substring name = lex_tokss (s->lexer);
5265           get->names = matrix_var_lookup (s, name);
5266           if (!get->names)
5267             get->names = matrix_var_create (s, name);
5268           lex_get (s->lexer);
5269         }
5270       else if (lex_match_id (s->lexer, "MISSING"))
5271         {
5272           lex_match (s->lexer, T_EQUALS);
5273           if (lex_match_id (s->lexer, "ACCEPT"))
5274             get->user.treatment = MGET_ACCEPT;
5275           else if (lex_match_id (s->lexer, "OMIT"))
5276             get->user.treatment = MGET_OMIT;
5277           else if (lex_is_number (s->lexer))
5278             {
5279               get->user.treatment = MGET_RECODE;
5280               get->user.substitute = lex_number (s->lexer);
5281               lex_get (s->lexer);
5282             }
5283           else
5284             {
5285               lex_error (s->lexer, NULL);
5286               goto error;
5287             }
5288         }
5289       else if (lex_match_id (s->lexer, "SYSMIS"))
5290         {
5291           lex_match (s->lexer, T_EQUALS);
5292           if (lex_match_id (s->lexer, "OMIT"))
5293             get->user.treatment = MGET_OMIT;
5294           else if (lex_is_number (s->lexer))
5295             {
5296               get->user.treatment = MGET_RECODE;
5297               get->user.substitute = lex_number (s->lexer);
5298               lex_get (s->lexer);
5299             }
5300           else
5301             {
5302               lex_error (s->lexer, NULL);
5303               goto error;
5304             }
5305         }
5306       else
5307         {
5308           lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5309                                "MISSING", "SYSMIS");
5310           goto error;
5311         }
5312     }
5313   return cmd;
5314
5315 error:
5316   matrix_cmd_destroy (cmd);
5317   return NULL;
5318 }
5319
5320 static void
5321 matrix_cmd_execute_get (struct get_command *get)
5322 {
5323   assert (get->file);           /* XXX */
5324
5325   struct dictionary *dict;
5326   struct casereader *reader = any_reader_open_and_decode (
5327     get->file, get->encoding, &dict, NULL);
5328   if (!reader)
5329     return;
5330
5331   const struct variable **vars = xnmalloc (
5332     get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5333     sizeof *vars);
5334   size_t n_vars = 0;
5335
5336   if (get->variables.n)
5337     {
5338       for (size_t i = 0; i < get->variables.n; i++)
5339         {
5340           const char *name = get->variables.strings[i];
5341           const struct variable *var = dict_lookup_var (dict, name);
5342           if (!var)
5343             {
5344               msg (SE, _("GET: Data file does not contain variable %s."),
5345                    name);
5346               dict_unref (dict);
5347               free (vars);
5348               return;
5349             }
5350           if (!var_is_numeric (var))
5351             {
5352               msg (SE, _("GET: Variable %s is not numeric."), name);
5353               dict_unref (dict);
5354               free (vars);
5355               return;
5356             }
5357           vars[n_vars++] = var;
5358         }
5359     }
5360   else
5361     {
5362       for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5363         {
5364           const struct variable *var = dict_get_var (dict, i);
5365           if (!var_is_numeric (var))
5366             {
5367               msg (SE, _("GET: Variable %s is not numeric."),
5368                    var_get_name (var));
5369               dict_unref (dict);
5370               free (vars);
5371               return;
5372             }
5373           vars[n_vars++] = var;
5374         }
5375     }
5376
5377   size_t n_rows = 0;
5378   gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5379   long long int casenum = 1;
5380   bool error = false;
5381   for (struct ccase *c = casereader_read (reader); c;
5382        c = casereader_read (reader), casenum++)
5383     {
5384       if (n_rows >= m->size1)
5385         {
5386           gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5387           for (size_t y = 0; y < n_rows; y++)
5388             for (size_t x = 0; x < n_vars; x++)
5389               gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5390           gsl_matrix_free (m);
5391           m = p;
5392         }
5393
5394       bool keep = true;
5395       for (size_t x = 0; x < n_vars; x++)
5396         {
5397           const struct variable *var = vars[x];
5398           double d = case_num (c, var);
5399           if (d == SYSMIS)
5400             {
5401               if (get->system.treatment == MGET_RECODE)
5402                 d = get->system.substitute;
5403               else if (get->system.treatment == MGET_OMIT)
5404                 keep = false;
5405               else
5406                 {
5407                   msg (SE, _("GET: Variable %s in case %lld "
5408                              "is system-missing."),
5409                        var_get_name (var), casenum);
5410                   error = true;
5411                 }
5412             }
5413           else if (var_is_num_missing (var, d, MV_USER))
5414             {
5415               if (get->user.treatment == MGET_RECODE)
5416                 d = get->user.substitute;
5417               else if (get->user.treatment == MGET_OMIT)
5418                 keep = false;
5419               else if (get->user.treatment != MGET_ACCEPT)
5420                 {
5421                   msg (SE, _("GET: Variable %s in case %lld has user-missing "
5422                              "value %g."),
5423                        var_get_name (var), casenum, d);
5424                   error = true;
5425                 }
5426             }
5427           gsl_matrix_set (m, n_rows, x, d);
5428         }
5429       case_unref (c);
5430       if (error)
5431         break;
5432       if (keep)
5433         n_rows++;
5434     }
5435   casereader_destroy (reader);
5436   if (!error)
5437     {
5438       m->size1 = n_rows;
5439       matrix_lvalue_evaluate_and_assign (get->dst, m);
5440     }
5441   else
5442     gsl_matrix_free (m);
5443   dict_unref (dict);
5444   free (vars);
5445 }
5446 \f
5447 static const char *
5448 match_rowtype (struct lexer *lexer)
5449 {
5450   static const char *rowtypes[] = {
5451     "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5452   };
5453   size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5454
5455   for (size_t i = 0; i < n_rowtypes; i++)
5456     if (lex_match_id (lexer, rowtypes[i]))
5457       return rowtypes[i];
5458
5459   lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5460   return NULL;
5461 }
5462
5463 static bool
5464 parse_var_names (struct lexer *lexer, struct string_array *sa)
5465 {
5466   lex_match (lexer, T_EQUALS);
5467
5468   string_array_clear (sa);
5469
5470   struct dictionary *dict = dict_create (get_default_encoding ());
5471   dict_create_var_assert (dict, "ROWTYPE_", 8);
5472   dict_create_var_assert (dict, "VARNAME_", 8);
5473   char **names;
5474   size_t n_names;
5475   bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5476                                   PV_NO_DUPLICATE | PV_NO_SCRATCH);
5477   dict_unref (dict);
5478
5479   if (ok)
5480     {
5481       string_array_clear (sa);
5482       sa->strings = names;
5483       sa->n = sa->allocated = n_names;
5484     }
5485   return ok;
5486 }
5487
5488 static void
5489 msave_common_uninit (struct msave_common *common)
5490 {
5491   if (common)
5492     {
5493       fh_unref (common->outfile);
5494       string_array_destroy (&common->variables);
5495       string_array_destroy (&common->fnames);
5496       string_array_destroy (&common->snames);
5497     }
5498 }
5499
5500 static bool
5501 compare_variables (const char *keyword,
5502                    const struct string_array *new,
5503                    const struct string_array *old)
5504 {
5505   if (new->n)
5506     {
5507       if (!old->n)
5508         {
5509           msg (SE, _("%s may only be specified on MSAVE if it was specified "
5510                      "on the first MSAVE within MATRIX."), keyword);
5511           return false;
5512         }
5513       else if (!string_array_equal_case (old, new))
5514         {
5515           msg (SE, _("%s must specify the same variables each time within "
5516                      "a given MATRIX."), keyword);
5517           return false;
5518         }
5519     }
5520   return true;
5521 }
5522
5523 static struct matrix_cmd *
5524 matrix_parse_msave (struct matrix_state *s)
5525 {
5526   struct msave_common common = { .outfile = NULL };
5527   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5528   *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5529
5530   struct msave_command *msave = &cmd->msave;
5531   if (lex_token (s->lexer) == T_ID)
5532     msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5533   msave->expr = matrix_parse_exp (s);
5534   if (!msave->expr)
5535     return NULL;
5536
5537   while (lex_match (s->lexer, T_SLASH))
5538     {
5539       if (lex_match_id (s->lexer, "TYPE"))
5540         {
5541           lex_match (s->lexer, T_EQUALS);
5542
5543           msave->rowtype = match_rowtype (s->lexer);
5544           if (!msave->rowtype)
5545             goto error;
5546         }
5547       else if (lex_match_id (s->lexer, "OUTFILE"))
5548         {
5549           lex_match (s->lexer, T_EQUALS);
5550
5551           fh_unref (common.outfile);
5552           common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5553           if (!common.outfile)
5554             goto error;
5555         }
5556       else if (lex_match_id (s->lexer, "VARIABLES"))
5557         {
5558           if (!parse_var_names (s->lexer, &common.variables))
5559             goto error;
5560         }
5561       else if (lex_match_id (s->lexer, "FNAMES"))
5562         {
5563           if (!parse_var_names (s->lexer, &common.fnames))
5564             goto error;
5565         }
5566       else if (lex_match_id (s->lexer, "SNAMES"))
5567         {
5568           if (!parse_var_names (s->lexer, &common.snames))
5569             goto error;
5570         }
5571       else if (lex_match_id (s->lexer, "SPLIT"))
5572         {
5573           lex_match (s->lexer, T_EQUALS);
5574
5575           matrix_expr_destroy (msave->splits);
5576           msave->splits = matrix_parse_exp (s);
5577           if (!msave->splits)
5578             goto error;
5579         }
5580       else if (lex_match_id (s->lexer, "FACTOR"))
5581         {
5582           lex_match (s->lexer, T_EQUALS);
5583
5584           matrix_expr_destroy (msave->factors);
5585           msave->factors = matrix_parse_exp (s);
5586           if (!msave->factors)
5587             goto error;
5588         }
5589       else
5590         {
5591           lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5592                                "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5593           goto error;
5594         }
5595     }
5596   if (!msave->rowtype)
5597     {
5598       lex_sbc_missing ("TYPE");
5599       goto error;
5600     }
5601   common.has_splits = msave->splits || common.snames.n;
5602   common.has_factors = msave->factors || common.fnames.n;
5603
5604   struct msave_common *c = s->common ? s->common : &common;
5605   if (c->fnames.n && !msave->factors)
5606     {
5607       msg (SE, _("FNAMES requires FACTOR."));
5608       goto error;
5609     }
5610   if (c->snames.n && !msave->splits)
5611     {
5612       msg (SE, _("SNAMES requires SPLIT."));
5613       goto error;
5614     }
5615   if (c->has_factors && !common.has_factors)
5616     {
5617       msg (SE, _("%s is required because it was present on the first "
5618                  "MSAVE in this MATRIX command."), "FACTOR");
5619       goto error;
5620     }
5621   if (c->has_splits && !common.has_splits)
5622     {
5623       msg (SE, _("%s is required because it was present on the first "
5624                  "MSAVE in this MATRIX command."), "SPLIT");
5625       goto error;
5626     }
5627
5628   if (!s->common)
5629     {
5630       if (!common.outfile)
5631         {
5632           lex_sbc_missing ("OUTFILE");
5633           goto error;
5634         }
5635       s->common = xmemdup (&common, sizeof common);
5636     }
5637   else
5638     {
5639       if (common.outfile)
5640         {
5641           bool same = common.outfile == s->common->outfile;
5642           fh_unref (common.outfile);
5643           if (!same)
5644             {
5645               msg (SE, _("OUTFILE must name the same file on each MSAVE "
5646                          "within a single MATRIX command."));
5647               goto error;
5648             }
5649         }
5650       if (!compare_variables ("VARIABLES",
5651                               &common.variables, &s->common->variables)
5652           || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5653           || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5654         goto error;
5655       msave_common_uninit (&common);
5656     }
5657   msave->common = s->common;
5658   if (!msave->varname_)
5659     msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5660   return cmd;
5661
5662 error:
5663   msave_common_uninit (&common);
5664   matrix_cmd_destroy (cmd);
5665   return NULL;
5666 }
5667
5668 static gsl_vector *
5669 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5670 {
5671   gsl_matrix *m = matrix_expr_evaluate (e);
5672   if (!m)
5673     return NULL;
5674
5675   if (!is_vector (m))
5676     {
5677       msg (SE, _("%s expression must evaluate to vector, "
5678                  "not a %zu×%zu matrix."),
5679            name, m->size1, m->size2);
5680       gsl_matrix_free (m);
5681       return NULL;
5682     }
5683
5684   return matrix_to_vector (m);
5685 }
5686
5687 static const char *
5688 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5689 {
5690   for (size_t i = 0; i < vars->n; i++)
5691     if (!dict_create_var (d, vars->strings[i], 0))
5692       return vars->strings[i];
5693   return NULL;
5694 }
5695
5696 static struct dictionary *
5697 msave_create_dict (const struct msave_common *common)
5698 {
5699   struct dictionary *dict = dict_create (get_default_encoding ());
5700
5701   const char *dup_split = msave_add_vars (dict, &common->snames);
5702   if (dup_split)
5703     {
5704       msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5705       goto error;
5706     }
5707
5708   dict_create_var_assert (dict, "ROWTYPE_", 8);
5709
5710   const char *dup_factor = msave_add_vars (dict, &common->fnames);
5711   if (dup_factor)
5712     {
5713       msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5714       goto error;
5715     }
5716
5717   dict_create_var_assert (dict, "VARNAME_", 8);
5718
5719   const char *dup_var = msave_add_vars (dict, &common->variables);
5720   if (dup_var)
5721     {
5722       msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5723       goto error;
5724     }
5725
5726   return dict;
5727
5728 error:
5729   dict_unref (dict);
5730   return NULL;
5731 }
5732
5733 static void
5734 matrix_cmd_execute_msave (struct msave_command *msave)
5735 {
5736   struct msave_common *common = msave->common;
5737   gsl_matrix *m = NULL;
5738   gsl_vector *factors = NULL;
5739   gsl_vector *splits = NULL;
5740
5741   m = matrix_expr_evaluate (msave->expr);
5742   if (!m)
5743     goto error;
5744
5745   if (!common->variables.n)
5746     for (size_t i = 0; i < m->size2; i++)
5747       string_array_append_nocopy (&common->variables,
5748                                   xasprintf ("COL%zu", i + 1));
5749
5750   if (m->size2 != common->variables.n)
5751     {
5752       msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5753            m->size2, common->variables.n);
5754       goto error;
5755     }
5756
5757   if (msave->factors)
5758     {
5759       factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5760       if (!factors)
5761         goto error;
5762
5763       if (!common->fnames.n)
5764         for (size_t i = 0; i < factors->size; i++)
5765           string_array_append_nocopy (&common->fnames,
5766                                       xasprintf ("FAC%zu", i + 1));
5767
5768       if (factors->size != common->fnames.n)
5769         {
5770           msg (SE, _("There are %zu factor variables, "
5771                      "but %zu split values were supplied."),
5772                common->fnames.n, factors->size);
5773           goto error;
5774         }
5775     }
5776
5777   if (msave->splits)
5778     {
5779       splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5780       if (!splits)
5781         goto error;
5782
5783       if (!common->fnames.n)
5784         for (size_t i = 0; i < splits->size; i++)
5785           string_array_append_nocopy (&common->fnames,
5786                                       xasprintf ("SPL%zu", i + 1));
5787
5788       if (splits->size != common->fnames.n)
5789         {
5790           msg (SE, _("There are %zu split variables, "
5791                      "but %zu split values were supplied."),
5792                common->fnames.n, splits->size);
5793           goto error;
5794         }
5795     }
5796
5797   if (!common->writer)
5798     {
5799       struct dictionary *dict = msave_create_dict (common);
5800       if (!dict)
5801         goto error;
5802
5803       common->writer = any_writer_open (common->outfile, dict);
5804       if (!common->writer)
5805         {
5806           dict_unref (dict);
5807           goto error;
5808         }
5809
5810       common->dict = dict;
5811     }
5812
5813   for (size_t y = 0; y < m->size1; y++)
5814     {
5815       struct ccase *c = case_create (dict_get_proto (common->dict));
5816       size_t idx = 0;
5817
5818       /* Split variables */
5819       if (splits)
5820         for (size_t i = 0; i < splits->size; i++)
5821           case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5822
5823       /* ROWTYPE_. */
5824       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5825                          msave->rowtype, ' ');
5826
5827       /* Factors. */
5828       if (factors)
5829         for (size_t i = 0; i < factors->size; i++)
5830           *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5831
5832       /* VARNAME_. */
5833       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5834                          msave->varname_, ' ');
5835
5836       /* Continuous variables. */
5837       for (size_t x = 0; x < m->size2; x++)
5838         case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5839       casewriter_write (common->writer, c);
5840     }
5841
5842 error:
5843   gsl_matrix_free (m);
5844   gsl_vector_free (factors);
5845   gsl_vector_free (splits);
5846 }
5847 \f
5848 static struct matrix_cmd *
5849 matrix_parse_mget (struct matrix_state *s)
5850 {
5851   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5852   *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5853
5854   struct mget_command *mget = &cmd->mget;
5855
5856   for (;;)
5857     {
5858       if (lex_match_id (s->lexer, "FILE"))
5859         {
5860           lex_match (s->lexer, T_EQUALS);
5861
5862           fh_unref (mget->file);
5863           mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5864           if (!mget->file)
5865             goto error;
5866         }
5867       else if (lex_match_id (s->lexer, "TYPE"))
5868         {
5869           lex_match (s->lexer, T_EQUALS);
5870           while (lex_token (s->lexer) != T_SLASH
5871                  && lex_token (s->lexer) != T_ENDCMD)
5872             {
5873               const char *rowtype = match_rowtype (s->lexer);
5874               if (!rowtype)
5875                 goto error;
5876
5877               stringi_set_insert (&mget->rowtypes, rowtype);
5878             }
5879         }
5880       else
5881         {
5882           lex_error_expecting (s->lexer, "FILE", "TYPE");
5883           goto error;
5884         }
5885       if (lex_token (s->lexer) == T_ENDCMD)
5886         break;
5887
5888       if (!lex_force_match (s->lexer, T_SLASH))
5889         goto error;
5890     }
5891   return cmd;
5892
5893 error:
5894   matrix_cmd_destroy (cmd);
5895   return NULL;
5896 }
5897
5898 static const struct variable *
5899 get_a8_var (const struct dictionary *d, const char *name)
5900 {
5901   const struct variable *v = dict_lookup_var (d, name);
5902   if (!v)
5903     {
5904       msg (SE, _("Matrix data file lacks %s variable."), name);
5905       return NULL;
5906     }
5907   if (var_get_width (v) != 8)
5908     {
5909       msg (SE, _("%s variable in matrix data file must be "
5910                  "8-byte string, but it has width %d."),
5911            name, var_get_width (v));
5912       return NULL;
5913     }
5914   return v;
5915 }
5916
5917 static bool
5918 is_rowtype (const union value *v, const char *rowtype)
5919 {
5920   struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
5921   ss_rtrim (&vs, ss_cstr (" "));
5922   return ss_equals_case (vs, ss_cstr (rowtype));
5923 }
5924
5925 static void
5926 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
5927                         const struct dictionary *d,
5928                         const struct variable *rowtype_var,
5929                         struct matrix_state *s, size_t si, size_t fi,
5930                         size_t cs, size_t cn)
5931 {
5932   if (!n_rows)
5933     return;
5934
5935   const union value *rowtype_ = case_data (rows[0], rowtype_var);
5936   const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
5937                              : is_rowtype (rowtype_, "CORR") ? "CR"
5938                              : is_rowtype (rowtype_, "MEAN") ? "MN"
5939                              : is_rowtype (rowtype_, "STDDEV") ? "SD"
5940                              : is_rowtype (rowtype_, "N") ? "NC"
5941                              : "CN");
5942
5943   struct string name = DS_EMPTY_INITIALIZER;
5944   ds_put_cstr (&name, name_prefix);
5945   if (fi > 0)
5946     ds_put_format (&name, "F%zu", fi);
5947   if (si > 0)
5948     ds_put_format (&name, "S%zu", si);
5949
5950   struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
5951   if (!mv)
5952     mv = matrix_var_create (s, ds_ss (&name));
5953   else if (mv->value)
5954     {
5955       msg (SW, _("Matrix data file contains variable with existing name %s."),
5956            ds_cstr (&name));
5957       goto exit;
5958     }
5959
5960   gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
5961   size_t n_missing = 0;
5962   for (size_t y = 0; y < n_rows; y++)
5963     {
5964       for (size_t x = 0; x < cn; x++)
5965         {
5966           struct variable *var = dict_get_var (d, cs + x);
5967           double value = case_num (rows[y], var);
5968           if (var_is_num_missing (var, value, MV_ANY))
5969             {
5970               n_missing++;
5971               value = 0.0;
5972             }
5973           gsl_matrix_set (m, y, x, value);
5974         }
5975     }
5976
5977   if (n_missing)
5978     msg (SE, ngettext ("Matrix data file variable %s contains a missing "
5979                        "value, which was treated as zero.",
5980                        "Matrix data file variable %s contains %zu missing "
5981                        "values, which were treated as zero.", n_missing),
5982          ds_cstr (&name), n_missing);
5983   mv->value = m;
5984
5985 exit:
5986   ds_destroy (&name);
5987   for (size_t y = 0; y < n_rows; y++)
5988     case_unref (rows[y]);
5989 }
5990
5991 static bool
5992 var_changed (const struct ccase *ca, const struct ccase *cb,
5993              const struct variable *var)
5994 {
5995   return !value_equal (case_data (ca, var), case_data (cb, var),
5996                        var_get_width (var));
5997 }
5998
5999 static bool
6000 vars_changed (const struct ccase *ca, const struct ccase *cb,
6001               const struct dictionary *d,
6002               size_t first_var, size_t n_vars)
6003 {
6004   for (size_t i = 0; i < n_vars; i++)
6005     {
6006       const struct variable *v = dict_get_var (d, first_var + i);
6007       if (var_changed (ca, cb, v))
6008         return true;
6009     }
6010   return false;
6011 }
6012
6013 static void
6014 matrix_cmd_execute_mget (struct mget_command *mget)
6015 {
6016   struct dictionary *d;
6017   struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6018                                                      &d, NULL);
6019   if (!r)
6020     return;
6021
6022   const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6023   const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6024   if (!rowtype_ || !varname_)
6025     goto exit;
6026
6027   if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6028     {
6029       msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6030       goto exit;
6031     }
6032   if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6033     {
6034       msg (SE, _("Matrix data file contains no continuous variables."));
6035       goto exit;
6036     }
6037
6038   for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6039     {
6040       const struct variable *v = dict_get_var (d, i);
6041       if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6042         {
6043           msg (SE,
6044                _("Matrix data file contains unexpected string variable %s."),
6045                var_get_name (v));
6046           goto exit;
6047         }
6048     }
6049
6050   /* SPLIT variables. */
6051   size_t ss = 0;
6052   size_t sn = var_get_dict_index (rowtype_);
6053   struct ccase *sc = NULL;
6054   size_t si = 0;
6055
6056   /* FACTOR variables. */
6057   size_t fs = var_get_dict_index (rowtype_) + 1;
6058   size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6059   struct ccase *fc = NULL;
6060   size_t fi = 0;
6061
6062   /* Continuous variables. */
6063   size_t cs = var_get_dict_index (varname_) + 1;
6064   size_t cn = dict_get_var_cnt (d) - cs;
6065   struct ccase *cc = NULL;
6066
6067   /* Matrix. */
6068   struct ccase **rows = NULL;
6069   size_t allocated_rows = 0;
6070   size_t n_rows = 0;
6071
6072   struct ccase *c;
6073   while ((c = casereader_read (r)) != NULL)
6074     {
6075       bool sd = vars_changed (sc, c, d, ss, sn);
6076       bool fd = sd || vars_changed (fc, c, d, fs, fn);
6077       bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6078       if (sd)
6079         {
6080           si++;
6081           case_unref (sc);
6082           sc = case_ref (c);
6083         }
6084       if (fd)
6085         {
6086           fi++;
6087           case_unref (fc);
6088           fc = case_ref (c);
6089         }
6090       if (md)
6091         {
6092           matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6093                                   mget->state, si, fi, cs, cn);
6094           n_rows = 0;
6095           case_unref (cc);
6096           cc = case_ref (c);
6097         }
6098
6099       if (n_rows >= allocated_rows)
6100         rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6101       rows[n_rows++] = c;
6102     }
6103   matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6104                           mget->state, si, fi, cs, cn);
6105   free (rows);
6106
6107 exit:
6108   casereader_destroy (r);
6109 }
6110 \f
6111 static bool
6112 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6113 {
6114   if (!lex_force_id (s->lexer))
6115     return false;
6116
6117   *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6118   if (!*varp)
6119     *varp = matrix_var_create (s, lex_tokss (s->lexer));
6120   lex_get (s->lexer);
6121   return true;
6122 }
6123
6124 static struct matrix_cmd *
6125 matrix_parse_eigen (struct matrix_state *s)
6126 {
6127   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6128   *cmd = (struct matrix_cmd) {
6129     .type = MCMD_EIGEN,
6130     .eigen = { .expr = NULL }
6131   };
6132
6133   struct eigen_command *eigen = &cmd->eigen;
6134   if (!lex_force_match (s->lexer, T_LPAREN))
6135     goto error;
6136   eigen->expr = matrix_parse_expr (s);
6137   if (!eigen->expr
6138       || !lex_force_match (s->lexer, T_COMMA)
6139       || !matrix_parse_dst_var (s, &eigen->evec)
6140       || !lex_force_match (s->lexer, T_COMMA)
6141       || !matrix_parse_dst_var (s, &eigen->eval)
6142       || !lex_force_match (s->lexer, T_RPAREN))
6143     goto error;
6144
6145   return cmd;
6146
6147 error:
6148   matrix_cmd_destroy (cmd);
6149   return NULL;
6150 }
6151
6152 static void
6153 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6154 {
6155   gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6156   if (!A)
6157     return;
6158   if (!is_symmetric (A))
6159     {
6160       msg (SE, _("Argument of EIGEN must be symmetric."));
6161       gsl_matrix_free (A);
6162       return;
6163     }
6164
6165   gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6166   gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6167   gsl_vector v_eval = to_vector (eval);
6168   gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6169   gsl_eigen_symmv (A, &v_eval, evec, w);
6170   gsl_eigen_symmv_free (w);
6171
6172   gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6173
6174   gsl_matrix_free (A);
6175
6176   gsl_matrix_free (eigen->eval->value);
6177   eigen->eval->value = eval;
6178
6179   gsl_matrix_free (eigen->evec->value);
6180   eigen->evec->value = evec;
6181 }
6182 \f
6183 static struct matrix_cmd *
6184 matrix_parse_setdiag (struct matrix_state *s)
6185 {
6186   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6187   *cmd = (struct matrix_cmd) {
6188     .type = MCMD_SETDIAG,
6189     .setdiag = { .dst = NULL }
6190   };
6191
6192   struct setdiag_command *setdiag = &cmd->setdiag;
6193   if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6194     goto error;
6195
6196   setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6197   if (!setdiag->dst)
6198     {
6199       lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6200       goto error;
6201     }
6202   lex_get (s->lexer);
6203
6204   if (!lex_force_match (s->lexer, T_COMMA))
6205     goto error;
6206
6207   setdiag->expr = matrix_parse_expr (s);
6208   if (!setdiag->expr)
6209     goto error;
6210
6211   if (!lex_force_match (s->lexer, T_RPAREN))
6212     goto error;
6213
6214   return cmd;
6215
6216 error:
6217   matrix_cmd_destroy (cmd);
6218   return NULL;
6219 }
6220
6221 static void
6222 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6223 {
6224   gsl_matrix *dst = setdiag->dst->value;
6225   if (!dst)
6226     {
6227       msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6228            setdiag->dst->name);
6229       return;
6230     }
6231
6232   gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6233   if (!src)
6234     return;
6235
6236   size_t n = MIN (dst->size1, dst->size2);
6237   if (is_scalar (src))
6238     {
6239       double d = to_scalar (src);
6240       for (size_t i = 0; i < n; i++)
6241         gsl_matrix_set (dst, i, i, d);
6242     }
6243   else if (is_vector (src))
6244     {
6245       gsl_vector v = to_vector (src);
6246       for (size_t i = 0; i < n && i < v.size; i++)
6247         gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6248     }
6249   else
6250     msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6251                "not a %zu×%zu matrix."),
6252          src->size1, src->size2);
6253   gsl_matrix_free (src);
6254 }
6255 \f
6256 static struct matrix_cmd *
6257 matrix_parse_svd (struct matrix_state *s)
6258 {
6259   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6260   *cmd = (struct matrix_cmd) {
6261     .type = MCMD_SVD,
6262     .svd = { .expr = NULL }
6263   };
6264
6265   struct svd_command *svd = &cmd->svd;
6266   if (!lex_force_match (s->lexer, T_LPAREN))
6267     goto error;
6268   svd->expr = matrix_parse_expr (s);
6269   if (!svd->expr
6270       || !lex_force_match (s->lexer, T_COMMA)
6271       || !matrix_parse_dst_var (s, &svd->u)
6272       || !lex_force_match (s->lexer, T_COMMA)
6273       || !matrix_parse_dst_var (s, &svd->s)
6274       || !lex_force_match (s->lexer, T_COMMA)
6275       || !matrix_parse_dst_var (s, &svd->v)
6276       || !lex_force_match (s->lexer, T_RPAREN))
6277     goto error;
6278
6279   return cmd;
6280
6281 error:
6282   matrix_cmd_destroy (cmd);
6283   return NULL;
6284 }
6285
6286 static void
6287 matrix_cmd_execute_svd (struct svd_command *svd)
6288 {
6289   gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6290   if (!m)
6291     return;
6292
6293   if (m->size1 >= m->size2)
6294     {
6295       gsl_matrix *A = m;
6296       gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6297       gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6298       gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6299       gsl_vector *work = gsl_vector_alloc (A->size2);
6300       gsl_linalg_SV_decomp (A, V, &Sv, work);
6301       gsl_vector_free (work);
6302
6303       matrix_var_set (svd->u, A);
6304       matrix_var_set (svd->s, S);
6305       matrix_var_set (svd->v, V);
6306     }
6307   else
6308     {
6309       gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6310       gsl_matrix_transpose_memcpy (At, m);
6311       gsl_matrix_free (m);
6312
6313       gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6314       gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6315       gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6316       gsl_vector *work = gsl_vector_alloc (At->size2);
6317       gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6318       gsl_vector_free (work);
6319
6320       matrix_var_set (svd->v, At);
6321       matrix_var_set (svd->s, St);
6322       matrix_var_set (svd->u, Vt);
6323     }
6324 }
6325 \f
6326 static bool
6327 matrix_cmd_execute (struct matrix_cmd *cmd)
6328 {
6329   switch (cmd->type)
6330     {
6331     case MCMD_COMPUTE:
6332       matrix_cmd_execute_compute (&cmd->compute);
6333       break;
6334
6335     case MCMD_PRINT:
6336       matrix_cmd_execute_print (&cmd->print);
6337       break;
6338
6339     case MCMD_DO_IF:
6340       return matrix_cmd_execute_do_if (&cmd->do_if);
6341
6342     case MCMD_LOOP:
6343       matrix_cmd_execute_loop (&cmd->loop);
6344       break;
6345
6346     case MCMD_BREAK:
6347       return false;
6348
6349     case MCMD_DISPLAY:
6350       matrix_cmd_execute_display (&cmd->display);
6351       break;
6352
6353     case MCMD_RELEASE:
6354       matrix_cmd_execute_release (&cmd->release);
6355       break;
6356
6357     case MCMD_SAVE:
6358       matrix_cmd_execute_save (&cmd->save);
6359       break;
6360
6361     case MCMD_READ:
6362       matrix_cmd_execute_read (&cmd->read);
6363       break;
6364
6365     case MCMD_WRITE:
6366       matrix_cmd_execute_write (&cmd->write);
6367       break;
6368
6369     case MCMD_GET:
6370       matrix_cmd_execute_get (&cmd->get);
6371       break;
6372
6373     case MCMD_MSAVE:
6374       matrix_cmd_execute_msave (&cmd->msave);
6375       break;
6376
6377     case MCMD_MGET:
6378       matrix_cmd_execute_mget (&cmd->mget);
6379       break;
6380
6381     case MCMD_EIGEN:
6382       matrix_cmd_execute_eigen (&cmd->eigen);
6383       break;
6384
6385     case MCMD_SETDIAG:
6386       matrix_cmd_execute_setdiag (&cmd->setdiag);
6387       break;
6388
6389     case MCMD_SVD:
6390       matrix_cmd_execute_svd (&cmd->svd);
6391       break;
6392     }
6393
6394   return true;
6395 }
6396
6397 static void
6398 matrix_cmds_uninit (struct matrix_cmds *cmds)
6399 {
6400   for (size_t i = 0; i < cmds->n; i++)
6401     matrix_cmd_destroy (cmds->commands[i]);
6402   free (cmds->commands);
6403 }
6404
6405 static void
6406 matrix_cmd_destroy (struct matrix_cmd *cmd)
6407 {
6408   if (!cmd)
6409     return;
6410
6411   switch (cmd->type)
6412     {
6413     case MCMD_COMPUTE:
6414       matrix_lvalue_destroy (cmd->compute.lvalue);
6415       matrix_expr_destroy (cmd->compute.rvalue);
6416       break;
6417
6418     case MCMD_PRINT:
6419       matrix_expr_destroy (cmd->print.expression);
6420       free (cmd->print.title);
6421       print_labels_destroy (cmd->print.rlabels);
6422       print_labels_destroy (cmd->print.clabels);
6423       break;
6424
6425     case MCMD_DO_IF:
6426       for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6427         {
6428           matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6429           matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6430         }
6431       free (cmd->do_if.clauses);
6432       break;
6433
6434     case MCMD_LOOP:
6435       matrix_expr_destroy (cmd->loop.start);
6436       matrix_expr_destroy (cmd->loop.end);
6437       matrix_expr_destroy (cmd->loop.increment);
6438       matrix_expr_destroy (cmd->loop.top_condition);
6439       matrix_expr_destroy (cmd->loop.bottom_condition);
6440       matrix_cmds_uninit (&cmd->loop.commands);
6441       break;
6442
6443     case MCMD_BREAK:
6444       break;
6445
6446     case MCMD_DISPLAY:
6447       break;
6448
6449     case MCMD_RELEASE:
6450       free (cmd->release.vars);
6451       break;
6452
6453     case MCMD_SAVE:
6454       matrix_expr_destroy (cmd->save.expression);
6455       fh_unref (cmd->save.outfile);
6456       string_array_destroy (cmd->save.variables);
6457       matrix_expr_destroy (cmd->save.names);
6458       stringi_set_destroy (&cmd->save.strings);
6459       break;
6460
6461     case MCMD_READ:
6462       matrix_lvalue_destroy (cmd->read.dst);
6463       matrix_expr_destroy (cmd->read.size);
6464       break;
6465
6466     case MCMD_WRITE:
6467       matrix_expr_destroy (cmd->write.expression);
6468       free (cmd->write.encoding);
6469       break;
6470
6471     case MCMD_GET:
6472       matrix_lvalue_destroy (cmd->get.dst);
6473       fh_unref (cmd->get.file);
6474       free (cmd->get.encoding);
6475       string_array_destroy (&cmd->get.variables);
6476       break;
6477
6478     case MCMD_MSAVE:
6479       free (cmd->msave.varname_);
6480       matrix_expr_destroy (cmd->msave.expr);
6481       matrix_expr_destroy (cmd->msave.factors);
6482       matrix_expr_destroy (cmd->msave.splits);
6483       break;
6484
6485     case MCMD_MGET:
6486       fh_unref (cmd->mget.file);
6487       stringi_set_destroy (&cmd->mget.rowtypes);
6488       break;
6489
6490     case MCMD_EIGEN:
6491       matrix_expr_destroy (cmd->eigen.expr);
6492       break;
6493
6494     case MCMD_SETDIAG:
6495       matrix_expr_destroy (cmd->setdiag.expr);
6496       break;
6497
6498     case MCMD_SVD:
6499       matrix_expr_destroy (cmd->svd.expr);
6500       break;
6501     }
6502   free (cmd);
6503 }
6504
6505 struct matrix_command_name
6506   {
6507     const char *name;
6508     struct matrix_cmd *(*parse) (struct matrix_state *);
6509   };
6510
6511 static const struct matrix_command_name *
6512 matrix_parse_command_name (struct lexer *lexer)
6513 {
6514   static const struct matrix_command_name commands[] = {
6515     { "COMPUTE", matrix_parse_compute },
6516     { "CALL EIGEN", matrix_parse_eigen },
6517     { "CALL SETDIAG", matrix_parse_setdiag },
6518     { "CALL SVD", matrix_parse_svd },
6519     { "PRINT", matrix_parse_print },
6520     { "DO IF", matrix_parse_do_if },
6521     { "LOOP", matrix_parse_loop },
6522     { "BREAK", matrix_parse_break },
6523     { "READ", matrix_parse_read },
6524     { "WRITE", matrix_parse_write },
6525     { "GET", matrix_parse_get },
6526     { "SAVE", matrix_parse_save },
6527     { "MGET", matrix_parse_mget },
6528     { "MSAVE", matrix_parse_msave },
6529     { "DISPLAY", matrix_parse_display },
6530     { "RELEASE", matrix_parse_release },
6531   };
6532   static size_t n = sizeof commands / sizeof *commands;
6533
6534   for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6535     if (lex_match_phrase (lexer, c->name))
6536       return c;
6537   return NULL;
6538 }
6539
6540 static struct matrix_cmd *
6541 matrix_parse_command (struct matrix_state *s)
6542 {
6543   size_t nesting_level = SIZE_MAX;
6544
6545   struct matrix_cmd *c = NULL;
6546   const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6547   if (!cmd)
6548     lex_error (s->lexer, _("Unknown matrix command."));
6549   else if (!cmd->parse)
6550     lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6551                cmd->name);
6552   else
6553     {
6554       nesting_level = output_open_group (
6555         group_item_create_nocopy (utf8_to_title (cmd->name),
6556                                   utf8_to_title (cmd->name)));
6557       c = cmd->parse (s);
6558     }
6559
6560   if (c)
6561     lex_end_of_command (s->lexer);
6562   lex_discard_rest_of_command (s->lexer);
6563   if (nesting_level != SIZE_MAX)
6564     output_close_groups (nesting_level);
6565
6566   return c;
6567 }
6568
6569 int
6570 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6571 {
6572   if (!lex_force_match (lexer, T_ENDCMD))
6573     return CMD_FAILURE;
6574
6575   struct matrix_state state = {
6576     .session = dataset_session (ds),
6577     .lexer = lexer,
6578     .vars = HMAP_INITIALIZER (state.vars),
6579   };
6580
6581   for (;;)
6582     {
6583       while (lex_match (lexer, T_ENDCMD))
6584         continue;
6585       if (lex_token (lexer) == T_STOP)
6586         {
6587           msg (SE, _("Unexpected end of input expecting matrix command."));
6588           break;
6589         }
6590
6591       if (lex_match_phrase (lexer, "END MATRIX"))
6592         break;
6593
6594       struct matrix_cmd *c = matrix_parse_command (&state);
6595       if (c)
6596         {
6597           matrix_cmd_execute (c);
6598           matrix_cmd_destroy (c);
6599         }
6600     }
6601
6602   struct matrix_var *var, *next;
6603   HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6604     {
6605       free (var->name);
6606       gsl_matrix_free (var->value);
6607       hmap_delete (&state.vars, &var->hmap_node);
6608       free (var);
6609     }
6610   hmap_destroy (&state.vars);
6611   fh_unref (state.prev_save_outfile);
6612   fh_unref (state.prev_write_outfile);
6613   if (state.common)
6614     {
6615       dict_unref (state.common->dict);
6616       casewriter_destroy (state.common->writer);
6617       free (state.common);
6618     }
6619   fh_unref (state.prev_read_file);
6620   for (size_t i = 0; i < state.n_read_files; i++)
6621     read_file_destroy (state.read_files[i]);
6622   free (state.read_files);
6623
6624   return CMD_SUCCESS;
6625 }