pintos-os.org Git - pspp/blob - 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, _("Argument to GSCH with dimensions (%zu,%zu) 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 matrices with dimensions (%zu,%zu) "
2257                  "and (%zu,%zu)."),
2258            matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2259       return NULL;
2260     }
2261 }
2262
2263 static gsl_matrix *
2264 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2265 {
2266   if (is_scalar (a) || is_scalar (b))
2267     return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2268
2269   if (a->size2 != b->size1)
2270     {
2271       msg (SE, _("Matrices with dimensions (%zu,%zu) and (%zu,%zu) are "
2272                  "not conformable for multiplication."),
2273            a->size1, a->size2, b->size1, b->size2);
2274       return NULL;
2275     }
2276
2277   gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2278   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2279   return c;
2280 }
2281
2282 static void
2283 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2284 {
2285   gsl_matrix *tmp = *a;
2286   *a = *b;
2287   *b = tmp;
2288 }
2289
2290 static void
2291 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2292             gsl_matrix **tmp)
2293 {
2294   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2295   swap_matrix (z, tmp);
2296 }
2297
2298 static void
2299 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2300 {
2301   mul_matrix (x, *x, *x, tmp);
2302 }
2303
2304 static gsl_matrix *
2305 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2306 {
2307   gsl_matrix *x = x_;
2308   if (x->size1 != x->size2)
2309     {
2310       msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2311                  "the left-hand size, not one with dimensions (%zu,%zu)."),
2312            x->size1, x->size2);
2313       return NULL;
2314     }
2315   if (!is_scalar (b))
2316     {
2317       msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2318                  "right-hand side, not a matrix with dimensions (%zu,%zu)."),
2319            b->size1, b->size2);
2320       return NULL;
2321     }
2322   double bf = to_scalar (b);
2323   if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2324     {
2325       msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2326                  "or outside the valid range."), bf);
2327       return NULL;
2328     }
2329   long int bl = bf;
2330
2331   gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2332   gsl_matrix *y = y_;
2333   gsl_matrix_set_identity (y);
2334   if (bl == 0)
2335     return y;
2336
2337   gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2338   gsl_matrix *t = t_;
2339   for (unsigned long int n = labs (bl); n > 1; n /= 2)
2340     if (n & 1)
2341       {
2342         mul_matrix (&y, x, y, &t);
2343         square_matrix (&x, &t);
2344       }
2345     else
2346       square_matrix (&x, &t);
2347
2348   mul_matrix (&y, x, y, &t);
2349   if (bf < 0)
2350     invert_matrix (y);
2351
2352   /* Garbage collection.
2353
2354      There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2355      a permutation of them.  We are returning one of them; that one must not be
2356      destroyed.  We must not destroy 'x_' because the caller owns it. */
2357   if (y != y_)
2358     gsl_matrix_free (y_);
2359   if (y != t_)
2360     gsl_matrix_free (t_);
2361
2362   return y;
2363 }
2364
2365 static gsl_matrix *
2366 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2367                           gsl_matrix *by_)
2368 {
2369   if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2370     {
2371       msg (SE, _("All operands of : operator must be scalars."));
2372       return NULL;
2373     }
2374
2375   long int start = to_scalar (start_);
2376   long int end = to_scalar (end_);
2377   long int by = by_ ? to_scalar (by_) : 1;
2378
2379   if (!by)
2380     {
2381       msg (SE, _("The increment operand to : must be nonzero."));
2382       return NULL;
2383     }
2384
2385   long int n = (end > start && by > 0 ? (end - start + by) / by
2386                 : end < start && by < 0 ? (start - end - by) / -by
2387                 : 0);
2388   gsl_matrix *m = gsl_matrix_alloc (1, n);
2389   for (long int i = 0; i < n; i++)
2390     gsl_matrix_set (m, 0, i, start + i * by);
2391   return m;
2392 }
2393
2394 static gsl_matrix *
2395 matrix_expr_evaluate_not (gsl_matrix *a)
2396 {
2397   for (size_t r = 0; r < a->size1; r++)
2398     for (size_t c = 0; c < a->size2; c++)
2399       {
2400         double *ae = gsl_matrix_ptr (a, r, c);
2401         *ae = !(*ae > 0);
2402       }
2403   return a;
2404 }
2405
2406 static gsl_matrix *
2407 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2408 {
2409   if (a->size1 != b->size1)
2410     {
2411       if (!a->size1 || !a->size2)
2412         return b;
2413       else if (!b->size1 || !b->size2)
2414         return a;
2415
2416       msg (SE, _("All columns in a matrix must have the same number of rows, "
2417                  "but this tries to paste matrices with %zu and %zu rows."),
2418            a->size1, b->size1);
2419       return NULL;
2420     }
2421
2422   gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2423   for (size_t y = 0; y < a->size1; y++)
2424     {
2425       for (size_t x = 0; x < a->size2; x++)
2426         gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2427       for (size_t x = 0; x < b->size2; x++)
2428         gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2429     }
2430   return c;
2431 }
2432
2433 static gsl_matrix *
2434 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2435 {
2436   if (a->size2 != b->size2)
2437     {
2438       if (!a->size1 || !a->size2)
2439         return b;
2440       else if (!b->size1 || !b->size2)
2441         return a;
2442
2443       msg (SE, _("All rows in a matrix must have the same number of columns, "
2444                  "but this tries to stack matrices with %zu and %zu columns."),
2445            a->size2, b->size2);
2446       return NULL;
2447     }
2448
2449   gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2450   for (size_t x = 0; x < a->size2; x++)
2451     {
2452       for (size_t y = 0; y < a->size1; y++)
2453         gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2454       for (size_t y = 0; y < b->size1; y++)
2455         gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2456     }
2457   return c;
2458 }
2459
2460 static gsl_vector *
2461 matrix_to_vector (gsl_matrix *m)
2462 {
2463   assert (m->owner);
2464   gsl_vector v = to_vector (m);
2465   assert (v.block == m->block || !v.block);
2466   assert (!v.owner);
2467   v.owner = 1;
2468   m->owner = 0;
2469   return xmemdup (&v, sizeof v);
2470 }
2471
2472 enum index_type {
2473   IV_ROW,
2474   IV_COLUMN,
2475   IV_VECTOR
2476 };
2477
2478 struct index_vector
2479   {
2480     size_t *indexes;
2481     size_t n;
2482   };
2483 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2484
2485 static bool
2486 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2487                                enum index_type index_type, size_t other_size,
2488                                struct index_vector *iv)
2489 {
2490   if (m)
2491     {
2492       if (!is_vector (m))
2493         {
2494           switch (index_type)
2495             {
2496             case IV_VECTOR:
2497               msg (SE, _("Vector index must be scalar or vector, not a "
2498                          "matrix with dimensions (%zu,%zu)."),
2499                    m->size1, m->size2);
2500               break;
2501
2502             case IV_ROW:
2503               msg (SE, _("Matrix row index must be scalar or vector, not a "
2504                          "matrix with dimensions (%zu,%zu)."),
2505                    m->size1, m->size2);
2506               break;
2507
2508             case IV_COLUMN:
2509               msg (SE, _("Matrix column index must be scalar or vector, not a "
2510                          "matrix with dimensions (%zu,%zu)."),
2511                    m->size1, m->size2);
2512               break;
2513             }
2514           return false;
2515         }
2516
2517       gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2518       *iv = (struct index_vector) {
2519         .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2520         .n = v.size,
2521       };
2522       for (size_t i = 0; i < v.size; i++)
2523         {
2524           double index = gsl_vector_get (&v, i);
2525           if (index < 1 || index >= size + 1)
2526             {
2527               switch (index_type)
2528                 {
2529                 case IV_VECTOR:
2530                   msg (SE, _("Index %g is out of range for vector "
2531                              "with %zu elements."), index, size);
2532                   break;
2533
2534                 case IV_ROW:
2535                   msg (SE, _("%g is not a valid row index for a matrix "
2536                              "with dimensions (%zu,%zu)."),
2537                        index, size, other_size);
2538                   break;
2539
2540                 case IV_COLUMN:
2541                   msg (SE, _("%g is not a valid column index for a matrix "
2542                              "with dimensions (%zu,%zu)."),
2543                        index, other_size, size);
2544                   break;
2545                 }
2546
2547               free (iv->indexes);
2548               return false;
2549             }
2550           iv->indexes[i] = index - 1;
2551         }
2552       return true;
2553     }
2554   else
2555     {
2556       *iv = (struct index_vector) {
2557         .indexes = xnmalloc (size, sizeof *iv->indexes),
2558         .n = size,
2559       };
2560       for (size_t i = 0; i < size; i++)
2561         iv->indexes[i] = i;
2562       return true;
2563     }
2564 }
2565
2566 static gsl_matrix *
2567 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2568 {
2569   if (!is_vector (sm))
2570     {
2571       msg (SE, _("Vector index operator must be applied to vector, "
2572                  "not a matrix with dimensions (%zu,%zu)."),
2573            sm->size1, sm->size2);
2574       return NULL;
2575     }
2576
2577   return sm;
2578 }
2579
2580 static gsl_matrix *
2581 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2582 {
2583   if (!matrix_expr_evaluate_vec_all (sm))
2584     return NULL;
2585
2586   gsl_vector sv = to_vector (sm);
2587   struct index_vector iv;
2588   if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2589     return NULL;
2590
2591   gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2592                                      sm->size1 == 1 ? iv.n : 1);
2593   gsl_vector dv = to_vector (dm);
2594   for (size_t dx = 0; dx < iv.n; dx++)
2595     {
2596       size_t sx = iv.indexes[dx];
2597       gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2598     }
2599   free (iv.indexes);
2600
2601   return dm;
2602 }
2603
2604 static gsl_matrix *
2605 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2606                                 gsl_matrix *im1)
2607 {
2608   struct index_vector iv0;
2609   if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2610     return NULL;
2611
2612   struct index_vector iv1;
2613   if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2614                                       &iv1))
2615     {
2616       free (iv0.indexes);
2617       return NULL;
2618     }
2619
2620   gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2621   for (size_t dy = 0; dy < iv0.n; dy++)
2622     {
2623       size_t sy = iv0.indexes[dy];
2624
2625       for (size_t dx = 0; dx < iv1.n; dx++)
2626         {
2627           size_t sx = iv1.indexes[dx];
2628           gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2629         }
2630     }
2631   free (iv0.indexes);
2632   free (iv1.indexes);
2633   return dm;
2634 }
2635
2636 #define F(NAME, PROTOTYPE)                              \
2637   static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2638     const char *name, gsl_matrix *[], size_t,           \
2639     matrix_proto_##PROTOTYPE *);
2640 MATRIX_FUNCTIONS
2641 #undef F
2642
2643 static bool
2644 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2645 {
2646   if (!is_scalar (subs[index]))
2647     {
2648       msg (SE, _("Function %s argument %zu must be a scalar, "
2649                  "but it has dimensions (%zu,%zu)."),
2650            name, index + 1, subs[index]->size1, subs[index]->size2);
2651       return false;
2652     }
2653   return true;
2654 }
2655
2656 static bool
2657 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2658 {
2659   if (!is_vector (subs[index]))
2660     {
2661       msg (SE, _("Function %s argument %zu must be a vector, "
2662                  "but it has dimensions (%zu,%zu)."),
2663            name, index + 1, subs[index]->size1, subs[index]->size2);
2664       return false;
2665     }
2666   return true;
2667 }
2668
2669 static bool
2670 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2671 {
2672   for (size_t i = 0; i < n_subs; i++)
2673     {
2674       if (!check_scalar_arg (name, subs, i))
2675         return false;
2676       d[i] = to_scalar (subs[i]);
2677     }
2678   return true;
2679 }
2680
2681 static gsl_matrix *
2682 matrix_expr_evaluate_m_d (const char *name,
2683                           gsl_matrix *subs[], size_t n_subs,
2684                           matrix_proto_m_d *f)
2685 {
2686   assert (n_subs == 1);
2687
2688   double d;
2689   return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2690 }
2691
2692 static gsl_matrix *
2693 matrix_expr_evaluate_m_dd (const char *name,
2694                            gsl_matrix *subs[], size_t n_subs,
2695                            matrix_proto_m_dd *f)
2696 {
2697   assert (n_subs == 2);
2698
2699   double d[2];
2700   return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2701 }
2702
2703 static gsl_matrix *
2704 matrix_expr_evaluate_m_ddd (const char *name,
2705                            gsl_matrix *subs[], size_t n_subs,
2706                            matrix_proto_m_ddd *f)
2707 {
2708   assert (n_subs == 3);
2709
2710   double d[3];
2711   return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2712 }
2713
2714 static gsl_matrix *
2715 matrix_expr_evaluate_m_m (const char *name UNUSED,
2716                           gsl_matrix *subs[], size_t n_subs,
2717                           matrix_proto_m_m *f)
2718 {
2719   assert (n_subs == 1);
2720   return f (subs[0]);
2721 }
2722
2723 static gsl_matrix *
2724 matrix_expr_evaluate_m_md (const char *name UNUSED,
2725                            gsl_matrix *subs[], size_t n_subs,
2726                            matrix_proto_m_md *f)
2727 {
2728   assert (n_subs == 2);
2729   if (!check_scalar_arg (name, subs, 1))
2730     return NULL;
2731   return f (subs[0], to_scalar (subs[1]));
2732 }
2733
2734 static gsl_matrix *
2735 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2736                             gsl_matrix *subs[], size_t n_subs,
2737                             matrix_proto_m_mdd *f)
2738 {
2739   assert (n_subs == 3);
2740   if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2741     return NULL;
2742   return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2743 }
2744
2745 static gsl_matrix *
2746 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2747                            gsl_matrix *subs[], size_t n_subs,
2748                            matrix_proto_m_mm *f)
2749 {
2750   assert (n_subs == 2);
2751   return f (subs[0], subs[1]);
2752 }
2753
2754 static gsl_matrix *
2755 matrix_expr_evaluate_m_v (const char *name,
2756                           gsl_matrix *subs[], size_t n_subs,
2757                           matrix_proto_m_v *f)
2758 {
2759   assert (n_subs == 1);
2760   if (!check_vector_arg (name, subs, 0))
2761     return NULL;
2762   gsl_vector v = to_vector (subs[0]);
2763   return f (&v);
2764 }
2765
2766 static gsl_matrix *
2767 matrix_expr_evaluate_d_m (const char *name UNUSED,
2768                           gsl_matrix *subs[], size_t n_subs,
2769                           matrix_proto_d_m *f)
2770 {
2771   assert (n_subs == 1);
2772
2773   gsl_matrix *m = gsl_matrix_alloc (1, 1);
2774   gsl_matrix_set (m, 0, 0, f (subs[0]));
2775   return m;
2776 }
2777
2778 static gsl_matrix *
2779 matrix_expr_evaluate_m_any (const char *name UNUSED,
2780                           gsl_matrix *subs[], size_t n_subs,
2781                           matrix_proto_m_any *f)
2782 {
2783   return f (subs, n_subs);
2784 }
2785
2786 static gsl_matrix *
2787 matrix_expr_evaluate_IDENT (const char *name,
2788                             gsl_matrix *subs[], size_t n_subs,
2789                             matrix_proto_IDENT *f)
2790 {
2791   assert (n_subs <= 2);
2792
2793   double d[2];
2794   if (!to_scalar_args (name, subs, n_subs, d))
2795     return NULL;
2796
2797   return f (d[0], d[n_subs - 1]);
2798 }
2799
2800 static gsl_matrix *
2801 matrix_expr_evaluate (const struct matrix_expr *e)
2802 {
2803   if (e->op == MOP_NUMBER)
2804     {
2805       gsl_matrix *m = gsl_matrix_alloc (1, 1);
2806       gsl_matrix_set (m, 0, 0, e->number);
2807       return m;
2808     }
2809   else if (e->op == MOP_VARIABLE)
2810     {
2811       const gsl_matrix *src = e->variable->value;
2812       if (!src)
2813         {
2814           msg (SE, _("Uninitialized variable %s used in expression."),
2815                e->variable->name);
2816           return NULL;
2817         }
2818
2819       gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2820       gsl_matrix_memcpy (dst, src);
2821       return dst;
2822     }
2823   else if (e->op == MOP_EOF)
2824     {
2825       struct dfm_reader *reader = read_file_open (e->eof);
2826       gsl_matrix *m = gsl_matrix_alloc (1, 1);
2827       gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2828       return m;
2829     }
2830
2831   enum { N_LOCAL = 3 };
2832   gsl_matrix *local_subs[N_LOCAL];
2833   gsl_matrix **subs = (e->n_subs < N_LOCAL
2834                        ? local_subs
2835                        : xmalloc (e->n_subs * sizeof *subs));
2836
2837   for (size_t i = 0; i < e->n_subs; i++)
2838     {
2839       subs[i] = matrix_expr_evaluate (e->subs[i]);
2840       if (!subs[i])
2841         {
2842           for (size_t j = 0; j < i; j++)
2843             gsl_matrix_free (subs[j]);
2844           if (subs != local_subs)
2845             free (subs);
2846           return NULL;
2847         }
2848     }
2849
2850   gsl_matrix *result = NULL;
2851   switch (e->op)
2852     {
2853 #define F(NAME, PROTOTYPE)                                              \
2854       case MOP_F_##NAME:                                                \
2855         result = matrix_expr_evaluate_##PROTOTYPE (#NAME,               \
2856                                                    subs, e->n_subs,     \
2857                                                    matrix_eval_##NAME); \
2858       break;
2859       MATRIX_FUNCTIONS
2860 #undef F
2861
2862     case MOP_NEGATE:
2863       gsl_matrix_scale (subs[0], -1.0);
2864       result = subs[0];
2865       break;
2866
2867     case MOP_ADD_ELEMS:
2868     case MOP_SUB_ELEMS:
2869     case MOP_MUL_ELEMS:
2870     case MOP_DIV_ELEMS:
2871     case MOP_EXP_ELEMS:
2872     case MOP_GT:
2873     case MOP_GE:
2874     case MOP_LT:
2875     case MOP_LE:
2876     case MOP_EQ:
2877     case MOP_NE:
2878     case MOP_AND:
2879     case MOP_OR:
2880     case MOP_XOR:
2881       result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2882       break;
2883
2884     case MOP_NOT:
2885       result = matrix_expr_evaluate_not (subs[0]);
2886       break;
2887
2888     case MOP_SEQ:
2889       result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2890       break;
2891
2892     case MOP_SEQ_BY:
2893       result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2894       break;
2895
2896     case MOP_MUL_MAT:
2897       result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2898       break;
2899
2900     case MOP_EXP_MAT:
2901       result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2902       break;
2903
2904     case MOP_PASTE_HORZ:
2905       result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2906       break;
2907
2908     case MOP_PASTE_VERT:
2909       result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2910       break;
2911
2912     case MOP_EMPTY:
2913       result = gsl_matrix_alloc (0, 0);
2914       break;
2915
2916     case MOP_VEC_INDEX:
2917       result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2918       break;
2919
2920     case MOP_VEC_ALL:
2921       result = matrix_expr_evaluate_vec_all (subs[0]);
2922       break;
2923
2924     case MOP_MAT_INDEX:
2925       result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2926       break;
2927
2928     case MOP_ROW_INDEX:
2929       result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2930       break;
2931
2932     case MOP_COL_INDEX:
2933       result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
2934       break;
2935
2936     case MOP_NUMBER:
2937     case MOP_VARIABLE:
2938     case MOP_EOF:
2939       NOT_REACHED ();
2940     }
2941
2942   for (size_t i = 0; i < e->n_subs; i++)
2943     if (subs[i] != result)
2944       gsl_matrix_free (subs[i]);
2945   if (subs != local_subs)
2946     free (subs);
2947   return result;
2948 }
2949
2950 static bool
2951 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
2952                              double *d)
2953 {
2954   gsl_matrix *m = matrix_expr_evaluate (e);
2955   if (!m)
2956     return false;
2957
2958   if (!is_scalar (m))
2959     {
2960       msg (SE, _("Expression for %s must evaluate to scalar, "
2961                  "not a matrix with dimensions (%zu,%zu)."),
2962            context, m->size1, m->size2);
2963       gsl_matrix_free (m);
2964       return false;
2965     }
2966
2967   *d = to_scalar (m);
2968   gsl_matrix_free (m);
2969   return true;
2970 }
2971
2972 static bool
2973 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
2974                               long int *integer)
2975 {
2976   double d;
2977   if (!matrix_expr_evaluate_scalar (e, context, &d))
2978     return false;
2979
2980   d = trunc (d);
2981   if (d < LONG_MIN || d > LONG_MAX)
2982     {
2983       msg (SE, _("Expression for %s is outside the integer range."), context);
2984       return false;
2985     }
2986   *integer = d;
2987   return true;
2988 }
2989 \f
2990 struct matrix_lvalue
2991   {
2992     struct matrix_var *var;
2993     struct matrix_expr *indexes[2];
2994     size_t n_indexes;
2995   };
2996
2997 static void
2998 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
2999 {
3000   if (lvalue)
3001     {
3002       for (size_t i = 0; i < lvalue->n_indexes; i++)
3003         matrix_expr_destroy (lvalue->indexes[i]);
3004       free (lvalue);
3005     }
3006 }
3007
3008 static struct matrix_lvalue *
3009 matrix_lvalue_parse (struct matrix_state *s)
3010 {
3011   struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3012   if (!lex_force_id (s->lexer))
3013     goto error;
3014   lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3015   if (lex_next_token (s->lexer, 1) == T_LPAREN)
3016     {
3017       if (!lvalue->var)
3018         {
3019           msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3020           free (lvalue);
3021           return NULL;
3022         }
3023
3024       lex_get_n (s->lexer, 2);
3025
3026       if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3027         goto error;
3028       lvalue->n_indexes++;
3029
3030       if (lex_match (s->lexer, T_COMMA))
3031         {
3032           if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3033             goto error;
3034           lvalue->n_indexes++;
3035         }
3036       if (!lex_force_match (s->lexer, T_RPAREN))
3037         goto error;
3038     }
3039   else
3040     {
3041       if (!lvalue->var)
3042         lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3043       lex_get (s->lexer);
3044     }
3045   return lvalue;
3046
3047 error:
3048   matrix_lvalue_destroy (lvalue);
3049   return NULL;
3050 }
3051
3052 static bool
3053 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3054                                enum index_type index_type, size_t other_size,
3055                                struct index_vector *iv)
3056 {
3057   gsl_matrix *m;
3058   if (e)
3059     {
3060       m = matrix_expr_evaluate (e);
3061       if (!m)
3062         return false;
3063     }
3064   else
3065     m = NULL;
3066
3067   bool ok = matrix_normalize_index_vector (m, size, index_type,
3068                                            other_size, iv);
3069   gsl_matrix_free (m);
3070   return ok;
3071 }
3072
3073 static bool
3074 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3075                              struct index_vector *iv, gsl_matrix *sm)
3076 {
3077   gsl_vector dv = to_vector (lvalue->var->value);
3078
3079   /* Convert source matrix 'sm' to source vector 'sv'. */
3080   if (!is_vector (sm))
3081     {
3082       msg (SE, _("Can't assign matrix with dimensions (%zu,%zu) to subvector."),
3083            sm->size1, sm->size2);
3084       return false;
3085     }
3086   gsl_vector sv = to_vector (sm);
3087   if (iv->n != sv.size)
3088     {
3089       msg (SE, _("Can't assign vector with %zu elements "
3090                  "to subvector with %zu."), sv.size, iv->n);
3091       return false;
3092     }
3093
3094   /* Assign elements. */
3095   for (size_t x = 0; x < iv->n; x++)
3096     gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3097   return true;
3098 }
3099
3100 static bool
3101 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3102                              struct index_vector *iv0,
3103                              struct index_vector *iv1, gsl_matrix *sm)
3104 {
3105   gsl_matrix *dm = lvalue->var->value;
3106
3107   /* Convert source matrix 'sm' to source vector 'sv'. */
3108   if (iv0->n != sm->size1)
3109     {
3110       msg (SE, _("Row index vector for assignment to %s has %zu elements "
3111                  "but source matrix has %zu rows."),
3112            lvalue->var->name, iv0->n, sm->size1);
3113       return false;
3114     }
3115   if (iv1->n != sm->size2)
3116     {
3117       msg (SE, _("Column index vector for assignment to %s has %zu elements "
3118                  "but source matrix has %zu columns."),
3119            lvalue->var->name, iv1->n, sm->size2);
3120       return false;
3121     }
3122
3123   /* Assign elements. */
3124   for (size_t y = 0; y < iv0->n; y++)
3125     {
3126       size_t f0 = iv0->indexes[y];
3127       for (size_t x = 0; x < iv1->n; x++)
3128         {
3129           size_t f1 = iv1->indexes[x];
3130           gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3131         }
3132     }
3133   return true;
3134 }
3135
3136 /* Takes ownership of M. */
3137 static bool
3138 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3139                       struct index_vector *iv0, struct index_vector *iv1,
3140                       gsl_matrix *sm)
3141 {
3142   if (!lvalue->n_indexes)
3143     {
3144       gsl_matrix_free (lvalue->var->value);
3145       lvalue->var->value = sm;
3146       return true;
3147     }
3148   else
3149     {
3150       bool ok = (lvalue->n_indexes == 1
3151                  ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3152                  : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3153       free (iv0->indexes);
3154       free (iv1->indexes);
3155       gsl_matrix_free (sm);
3156       return ok;
3157     }
3158 }
3159
3160 static bool
3161 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3162                         struct index_vector *iv0,
3163                         struct index_vector *iv1)
3164 {
3165   *iv0 = INDEX_VECTOR_INIT;
3166   *iv1 = INDEX_VECTOR_INIT;
3167   if (!lvalue->n_indexes)
3168     return true;
3169
3170   /* Validate destination matrix exists and has the right shape. */
3171   gsl_matrix *dm = lvalue->var->value;
3172   if (!dm)
3173     {
3174       msg (SE, _("Undefined variable %s."), lvalue->var->name);
3175       return false;
3176     }
3177   else if (dm->size1 == 0 || dm->size2 == 0)
3178     {
3179       msg (SE, _("Cannot index matrix with dimensions (%zu,%zu)."),
3180            dm->size1, dm->size2);
3181       return false;
3182     }
3183   else if (lvalue->n_indexes == 1)
3184     {
3185       if (!is_vector (dm))
3186         {
3187           msg (SE, _("Can't use vector indexing on matrix %s with "
3188                      "dimensions (%zu,%zu)."),
3189                lvalue->var->name, dm->size1, dm->size2);
3190           return false;
3191         }
3192       return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3193                                             MAX (dm->size1, dm->size2),
3194                                             IV_VECTOR, 0, iv0);
3195     }
3196   else
3197     {
3198       assert (lvalue->n_indexes == 2);
3199       if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3200                                           IV_ROW, dm->size2, iv0))
3201         return false;
3202
3203       if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3204                                           IV_COLUMN, dm->size1, iv1))
3205         {
3206           free (iv0->indexes);
3207           return false;
3208         }
3209       return true;
3210     }
3211 }
3212
3213 /* Takes ownership of M. */
3214 static bool
3215 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3216 {
3217   struct index_vector iv0, iv1;
3218   if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3219     {
3220       gsl_matrix_free (sm);
3221       return false;
3222     }
3223
3224   return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3225 }
3226
3227 \f
3228 /* Matrix command. */
3229
3230 struct matrix_cmds
3231   {
3232     struct matrix_cmd **commands;
3233     size_t n;
3234   };
3235
3236 static void matrix_cmds_uninit (struct matrix_cmds *);
3237
3238 struct matrix_cmd
3239   {
3240     enum matrix_cmd_type
3241       {
3242         MCMD_COMPUTE,
3243         MCMD_PRINT,
3244         MCMD_DO_IF,
3245         MCMD_LOOP,
3246         MCMD_BREAK,
3247         MCMD_DISPLAY,
3248         MCMD_RELEASE,
3249         MCMD_SAVE,
3250         MCMD_READ,
3251         MCMD_WRITE,
3252         MCMD_GET,
3253         MCMD_MSAVE,
3254         MCMD_MGET,
3255         MCMD_EIGEN,
3256         MCMD_SETDIAG,
3257         MCMD_SVD,
3258       }
3259     type;
3260
3261     union
3262       {
3263         struct compute_command
3264           {
3265             struct matrix_lvalue *lvalue;
3266             struct matrix_expr *rvalue;
3267           }
3268         compute;
3269
3270         struct print_command
3271           {
3272             struct matrix_expr *expression;
3273             bool use_default_format;
3274             struct fmt_spec format;
3275             char *title;
3276             int space;          /* -1 means NEWPAGE. */
3277
3278             struct print_labels
3279               {
3280                 struct string_array literals; /* CLABELS/RLABELS. */
3281                 struct matrix_expr *expr;     /* CNAMES/RNAMES. */
3282               }
3283             *rlabels, *clabels;
3284           }
3285         print;
3286
3287         struct do_if_command
3288           {
3289             struct do_if_clause
3290               {
3291                 struct matrix_expr *condition;
3292                 struct matrix_cmds commands;
3293               }
3294             *clauses;
3295             size_t n_clauses;
3296           }
3297         do_if;
3298
3299         struct loop_command
3300           {
3301             /* Index. */
3302             struct matrix_var *var;
3303             struct matrix_expr *start;
3304             struct matrix_expr *end;
3305             struct matrix_expr *increment;
3306
3307             /* Loop conditions. */
3308             struct matrix_expr *top_condition;
3309             struct matrix_expr *bottom_condition;
3310
3311             /* Commands. */
3312             struct matrix_cmds commands;
3313           }
3314         loop;
3315
3316         struct display_command
3317           {
3318             struct matrix_state *state;
3319             bool dictionary;
3320             bool status;
3321           }
3322         display;
3323
3324         struct release_command
3325           {
3326             struct matrix_var **vars;
3327             size_t n_vars;
3328           }
3329         release;
3330
3331         struct save_command
3332           {
3333             struct matrix_expr *expression;
3334             struct file_handle *outfile;
3335             struct string_array *variables;
3336             struct matrix_expr *names;
3337             struct stringi_set strings;
3338           }
3339         save;
3340
3341         struct read_command
3342           {
3343             struct read_file *rf;
3344             struct matrix_lvalue *dst;
3345             struct matrix_expr *size;
3346             int c1, c2;
3347             enum fmt_type format;
3348             int w;
3349             int d;
3350             bool symmetric;
3351             bool reread;
3352           }
3353         read;
3354
3355         struct write_command
3356           {
3357             struct matrix_expr *expression;
3358             struct file_handle *outfile;
3359             char *encoding;
3360             int c1, c2;
3361             enum fmt_type format;
3362             int w;
3363             int d;
3364             bool triangular;
3365             bool hold;
3366           }
3367         write;
3368
3369         struct get_command
3370           {
3371             struct matrix_lvalue *dst;
3372             struct file_handle *file;
3373             char *encoding;
3374             struct string_array variables;
3375             struct matrix_var *names;
3376
3377             /* Treatment of missing values. */
3378             struct
3379               {
3380                 enum
3381                   {
3382                     MGET_ERROR,  /* Flag error on command. */
3383                     MGET_ACCEPT, /* Accept without change, user-missing only. */
3384                     MGET_OMIT,   /* Drop this case. */
3385                     MGET_RECODE  /* Recode to 'substitute'. */
3386                   }
3387                 treatment;
3388                 double substitute; /* MGET_RECODE only. */
3389               }
3390             user, system;
3391           }
3392         get;
3393
3394         struct msave_command
3395           {
3396             struct msave_common *common;
3397             char *varname_;
3398             struct matrix_expr *expr;
3399             const char *rowtype;
3400             struct matrix_expr *factors;
3401             struct matrix_expr *splits;
3402           }
3403          msave;
3404
3405         struct mget_command
3406           {
3407             struct matrix_state *state;
3408             struct file_handle *file;
3409             struct stringi_set rowtypes;
3410           }
3411         mget;
3412
3413         struct eigen_command
3414           {
3415             struct matrix_expr *expr;
3416             struct matrix_var *evec;
3417             struct matrix_var *eval;
3418           }
3419         eigen;
3420
3421         struct setdiag_command
3422           {
3423             struct matrix_var *dst;
3424             struct matrix_expr *expr;
3425           }
3426         setdiag;
3427
3428         struct svd_command
3429           {
3430             struct matrix_expr *expr;
3431             struct matrix_var *u;
3432             struct matrix_var *s;
3433             struct matrix_var *v;
3434           }
3435         svd;
3436       };
3437   };
3438
3439 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3440 static bool matrix_cmd_execute (struct matrix_cmd *);
3441 static void matrix_cmd_destroy (struct matrix_cmd *);
3442
3443 \f
3444 static struct matrix_cmd *
3445 matrix_parse_compute (struct matrix_state *s)
3446 {
3447   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3448   *cmd = (struct matrix_cmd) {
3449     .type = MCMD_COMPUTE,
3450     .compute = { .lvalue = NULL },
3451   };
3452
3453   struct compute_command *compute = &cmd->compute;
3454   compute->lvalue = matrix_lvalue_parse (s);
3455   if (!compute->lvalue)
3456     goto error;
3457
3458   if (!lex_force_match (s->lexer, T_EQUALS))
3459     goto error;
3460
3461   compute->rvalue = matrix_parse_expr (s);
3462   if (!compute->rvalue)
3463     goto error;
3464
3465   return cmd;
3466
3467 error:
3468   matrix_cmd_destroy (cmd);
3469   return NULL;
3470 }
3471
3472 static void
3473 matrix_cmd_execute_compute (struct compute_command *compute)
3474 {
3475   gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3476   if (!value)
3477     return;
3478
3479   matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3480 }
3481 \f
3482 static void
3483 print_labels_destroy (struct print_labels *labels)
3484 {
3485   if (labels)
3486     {
3487       string_array_destroy (&labels->literals);
3488       matrix_expr_destroy (labels->expr);
3489       free (labels);
3490     }
3491 }
3492
3493 static void
3494 parse_literal_print_labels (struct matrix_state *s,
3495                             struct print_labels **labelsp)
3496 {
3497   lex_match (s->lexer, T_EQUALS);
3498   print_labels_destroy (*labelsp);
3499   *labelsp = xzalloc (sizeof **labelsp);
3500   while (lex_token (s->lexer) != T_SLASH
3501          && lex_token (s->lexer) != T_ENDCMD
3502          && lex_token (s->lexer) != T_STOP)
3503     {
3504       struct string label = DS_EMPTY_INITIALIZER;
3505       while (lex_token (s->lexer) != T_COMMA
3506              && lex_token (s->lexer) != T_SLASH
3507              && lex_token (s->lexer) != T_ENDCMD
3508              && lex_token (s->lexer) != T_STOP)
3509         {
3510           if (!ds_is_empty (&label))
3511             ds_put_byte (&label, ' ');
3512
3513           if (lex_token (s->lexer) == T_STRING)
3514             ds_put_cstr (&label, lex_tokcstr (s->lexer));
3515           else
3516             {
3517               char *rep = lex_next_representation (s->lexer, 0, 0);
3518               ds_put_cstr (&label, rep);
3519               free (rep);
3520             }
3521           lex_get (s->lexer);
3522         }
3523       string_array_append_nocopy (&(*labelsp)->literals,
3524                                   ds_steal_cstr (&label));
3525
3526       if (!lex_match (s->lexer, T_COMMA))
3527         break;
3528     }
3529 }
3530
3531 static bool
3532 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3533 {
3534   lex_match (s->lexer, T_EQUALS);
3535   struct matrix_expr *e = matrix_parse_exp (s);
3536   if (!e)
3537     return false;
3538
3539   print_labels_destroy (*labelsp);
3540   *labelsp = xzalloc (sizeof **labelsp);
3541   (*labelsp)->expr = e;
3542   return true;
3543 }
3544
3545 static struct matrix_cmd *
3546 matrix_parse_print (struct matrix_state *s)
3547 {
3548   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3549   *cmd = (struct matrix_cmd) {
3550     .type = MCMD_PRINT,
3551     .print = {
3552       .use_default_format = true,
3553     }
3554   };
3555
3556   if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3557     {
3558       size_t depth = 0;
3559       for (size_t i = 0; ; i++)
3560         {
3561           enum token_type t = lex_next_token (s->lexer, i);
3562           if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3563             depth++;
3564           else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3565             depth--;
3566           else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3567             {
3568               if (i > 0)
3569                 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3570               break;
3571             }
3572         }
3573
3574       cmd->print.expression = matrix_parse_exp (s);
3575       if (!cmd->print.expression)
3576         goto error;
3577     }
3578
3579   while (lex_match (s->lexer, T_SLASH))
3580     {
3581       if (lex_match_id (s->lexer, "FORMAT"))
3582         {
3583           lex_match (s->lexer, T_EQUALS);
3584           if (!parse_format_specifier (s->lexer, &cmd->print.format))
3585             goto error;
3586           cmd->print.use_default_format = false;
3587         }
3588       else if (lex_match_id (s->lexer, "TITLE"))
3589         {
3590           lex_match (s->lexer, T_EQUALS);
3591           if (!lex_force_string (s->lexer))
3592             goto error;
3593           free (cmd->print.title);
3594           cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3595           lex_get (s->lexer);
3596         }
3597       else if (lex_match_id (s->lexer, "SPACE"))
3598         {
3599           lex_match (s->lexer, T_EQUALS);
3600           if (lex_match_id (s->lexer, "NEWPAGE"))
3601             cmd->print.space = -1;
3602           else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3603             {
3604               cmd->print.space = lex_integer (s->lexer);
3605               lex_get (s->lexer);
3606             }
3607           else
3608             goto error;
3609         }
3610       else if (lex_match_id (s->lexer, "RLABELS"))
3611         parse_literal_print_labels (s, &cmd->print.rlabels);
3612       else if (lex_match_id (s->lexer, "CLABELS"))
3613         parse_literal_print_labels (s, &cmd->print.clabels);
3614       else if (lex_match_id (s->lexer, "RNAMES"))
3615         {
3616           if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3617             goto error;
3618         }
3619       else if (lex_match_id (s->lexer, "CNAMES"))
3620         {
3621           if (!parse_expr_print_labels (s, &cmd->print.clabels))
3622             goto error;
3623         }
3624       else
3625         {
3626           lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3627                                "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3628           goto error;
3629         }
3630
3631     }
3632   return cmd;
3633
3634 error:
3635   matrix_cmd_destroy (cmd);
3636   return NULL;
3637 }
3638
3639 static bool
3640 matrix_is_integer (const gsl_matrix *m)
3641 {
3642   for (size_t y = 0; y < m->size1; y++)
3643     for (size_t x = 0; x < m->size2; x++)
3644       {
3645         double d = gsl_matrix_get (m, y, x);
3646         if (d != floor (d))
3647           return false;
3648       }
3649   return true;
3650 }
3651
3652 static double
3653 matrix_max_magnitude (const gsl_matrix *m)
3654 {
3655   double max = 0.0;
3656   for (size_t y = 0; y < m->size1; y++)
3657     for (size_t x = 0; x < m->size2; x++)
3658       {
3659         double d = fabs (gsl_matrix_get (m, y, x));
3660         if (d > max)
3661           max = d;
3662       }
3663   return max;
3664 }
3665
3666 static bool
3667 format_fits (struct fmt_spec format, double x)
3668 {
3669   char *s = data_out (&(union value) { .f = x }, NULL,
3670                       &format, settings_get_fmt_settings ());
3671   bool fits = *s != '*' && !strchr (s, 'E');
3672   free (s);
3673   return fits;
3674 }
3675
3676 static struct fmt_spec
3677 default_format (const gsl_matrix *m, int *log_scale)
3678 {
3679   *log_scale = 0;
3680
3681   double max = matrix_max_magnitude (m);
3682
3683   if (matrix_is_integer (m))
3684     for (int w = 1; w <= 12; w++)
3685       {
3686         struct fmt_spec format = { .type = FMT_F, .w = w };
3687         if (format_fits (format, -max))
3688           return format;
3689       }
3690
3691   if (max >= 1e9 || max <= 1e-4)
3692     {
3693       char s[64];
3694       snprintf (s, sizeof s, "%.1e", max);
3695
3696       const char *e = strchr (s, 'e');
3697       if (e)
3698         *log_scale = atoi (e + 1);
3699     }
3700
3701   return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3702 }
3703
3704 static char *
3705 trimmed_string (double d)
3706 {
3707   struct substring s = ss_buffer ((char *) &d, sizeof d);
3708   ss_rtrim (&s, ss_cstr (" "));
3709   return ss_xstrdup (s);
3710 }
3711
3712 static struct string_array *
3713 print_labels_get (const struct print_labels *labels, size_t n,
3714                   const char *prefix, bool truncate)
3715 {
3716   if (!labels)
3717     return NULL;
3718
3719   struct string_array *out = xzalloc (sizeof *out);
3720   if (labels->literals.n)
3721     string_array_clone (out, &labels->literals);
3722   else if (labels->expr)
3723     {
3724       gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3725       if (m && is_vector (m))
3726         {
3727           gsl_vector v = to_vector (m);
3728           for (size_t i = 0; i < v.size; i++)
3729             string_array_append_nocopy (out, trimmed_string (
3730                                           gsl_vector_get (&v, i)));
3731         }
3732       gsl_matrix_free (m);
3733     }
3734
3735   while (out->n < n)
3736     {
3737       if (labels->expr)
3738         string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3739                                                     out->n + 1));
3740       else
3741         string_array_append (out, "");
3742     }
3743   while (out->n > n)
3744     string_array_delete (out, out->n - 1);
3745
3746   if (truncate)
3747     for (size_t i = 0; i < out->n; i++)
3748       {
3749         char *s = out->strings[i];
3750         s[strnlen (s, 8)] = '\0';
3751       }
3752
3753   return out;
3754 }
3755
3756 static void
3757 matrix_cmd_print_space (int space)
3758 {
3759   if (space < 0)
3760     output_item_submit (page_break_item_create ());
3761   for (int i = 0; i < space; i++)
3762     output_log ("%s", "");
3763 }
3764
3765 static void
3766 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3767                        struct fmt_spec format, int log_scale)
3768 {
3769   matrix_cmd_print_space (print->space);
3770   if (print->title)
3771     output_log ("%s", print->title);
3772   if (log_scale != 0)
3773     output_log ("  10 ** %d   X", log_scale);
3774
3775   struct string_array *clabels = print_labels_get (print->clabels,
3776                                                    m->size2, "col", true);
3777   if (clabels && format.w < 8)
3778     format.w = 8;
3779   struct string_array *rlabels = print_labels_get (print->rlabels,
3780                                                    m->size1, "row", true);
3781
3782   if (clabels)
3783     {
3784       struct string line = DS_EMPTY_INITIALIZER;
3785       if (rlabels)
3786         ds_put_byte_multiple (&line, ' ', 8);
3787       for (size_t x = 0; x < m->size2; x++)
3788         ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3789       output_log_nocopy (ds_steal_cstr (&line));
3790     }
3791
3792   double scale = pow (10.0, log_scale);
3793   bool numeric = fmt_is_numeric (format.type);
3794   for (size_t y = 0; y < m->size1; y++)
3795     {
3796       struct string line = DS_EMPTY_INITIALIZER;
3797       if (rlabels)
3798         ds_put_format (&line, "%-8s", rlabels->strings[y]);
3799
3800       for (size_t x = 0; x < m->size2; x++)
3801         {
3802           double f = gsl_matrix_get (m, y, x);
3803           char *s = (numeric
3804                      ? data_out (&(union value) { .f = f / scale}, NULL,
3805                                  &format, settings_get_fmt_settings ())
3806                      : trimmed_string (f));
3807           ds_put_format (&line, " %s", s);
3808           free (s);
3809         }
3810       output_log_nocopy (ds_steal_cstr (&line));
3811     }
3812
3813   string_array_destroy (rlabels);
3814   free (rlabels);
3815   string_array_destroy (clabels);
3816   free (clabels);
3817 }
3818
3819 static void
3820 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3821                         const struct print_labels *print_labels, size_t n,
3822                         const char *name, const char *prefix)
3823 {
3824   struct string_array *labels = print_labels_get (print_labels, n, prefix,
3825                                                   false);
3826   struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3827   for (size_t i = 0; i < n; i++)
3828     pivot_category_create_leaf (
3829       d->root, (labels
3830                 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3831                 : pivot_value_new_integer (i + 1)));
3832   if (!labels)
3833     d->hide_all_labels = true;
3834   string_array_destroy (labels);
3835   free (labels);
3836 }
3837
3838 static void
3839 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3840                          struct fmt_spec format, int log_scale)
3841 {
3842   struct pivot_table *table = pivot_table_create__ (
3843     pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3844     "Matrix Print");
3845
3846   create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3847                           N_("Rows"), "row");
3848   create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3849                           N_("Columns"), "col");
3850
3851   struct pivot_footnote *footnote = NULL;
3852   if (log_scale != 0)
3853     {
3854       char *s = xasprintf ("× 10**%d", log_scale);
3855       footnote = pivot_table_create_footnote (
3856         table, pivot_value_new_user_text_nocopy (s));
3857     }
3858
3859   double scale = pow (10.0, log_scale);
3860   bool numeric = fmt_is_numeric (format.type);
3861   for (size_t y = 0; y < m->size1; y++)
3862     for (size_t x = 0; x < m->size2; x++)
3863       {
3864         double f = gsl_matrix_get (m, y, x);
3865         struct pivot_value *v;
3866         if (numeric)
3867           {
3868             v = pivot_value_new_number (f / scale);
3869             v->numeric.format = format;
3870           }
3871         else
3872           v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3873         if (footnote)
3874           pivot_value_add_footnote (v, footnote);
3875         pivot_table_put2 (table, y, x, v);
3876       }
3877
3878   pivot_table_submit (table);
3879 }
3880
3881 static void
3882 matrix_cmd_execute_print (const struct print_command *print)
3883 {
3884   if (print->expression)
3885     {
3886       gsl_matrix *m = matrix_expr_evaluate (print->expression);
3887       if (!m)
3888         return;
3889
3890       int log_scale = 0;
3891       struct fmt_spec format = (print->use_default_format
3892                                 ? default_format (m, &log_scale)
3893                                 : print->format);
3894
3895       if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3896         matrix_cmd_print_text (print, m, format, log_scale);
3897       else
3898         matrix_cmd_print_tables (print, m, format, log_scale);
3899
3900       gsl_matrix_free (m);
3901     }
3902   else
3903     {
3904       matrix_cmd_print_space (print->space);
3905       if (print->title)
3906         output_log ("%s", print->title);
3907     }
3908 }
3909 \f
3910 /* DO IF. */
3911
3912 static bool
3913 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3914                        const char *command_name,
3915                        const char *stop1, const char *stop2)
3916 {
3917   lex_end_of_command (s->lexer);
3918   lex_discard_rest_of_command (s->lexer);
3919
3920   size_t allocated = 0;
3921   for (;;)
3922     {
3923       while (lex_token (s->lexer) == T_ENDCMD)
3924         lex_get (s->lexer);
3925
3926       if (lex_at_phrase (s->lexer, stop1)
3927           || (stop2 && lex_at_phrase (s->lexer, stop2)))
3928         return true;
3929
3930       if (lex_at_phrase (s->lexer, "END MATRIX"))
3931         {
3932           msg (SE, _("Premature END MATRIX within %s."), command_name);
3933           return false;
3934         }
3935
3936       struct matrix_cmd *cmd = matrix_parse_command (s);
3937       if (!cmd)
3938         return false;
3939
3940       if (c->n >= allocated)
3941         c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
3942       c->commands[c->n++] = cmd;
3943     }
3944 }
3945
3946 static bool
3947 matrix_parse_do_if_clause (struct matrix_state *s,
3948                            struct do_if_command *ifc,
3949                            bool parse_condition,
3950                            size_t *allocated_clauses)
3951 {
3952   if (ifc->n_clauses >= *allocated_clauses)
3953     ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
3954                                sizeof *ifc->clauses);
3955   struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
3956   *c = (struct do_if_clause) { .condition = NULL };
3957
3958   if (parse_condition)
3959     {
3960       c->condition = matrix_parse_expr (s);
3961       if (!c->condition)
3962         return false;
3963     }
3964
3965   return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
3966 }
3967
3968 static struct matrix_cmd *
3969 matrix_parse_do_if (struct matrix_state *s)
3970 {
3971   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3972   *cmd = (struct matrix_cmd) {
3973     .type = MCMD_DO_IF,
3974     .do_if = { .n_clauses = 0 }
3975   };
3976
3977   size_t allocated_clauses = 0;
3978   do
3979     {
3980       if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
3981         goto error;
3982     }
3983   while (lex_match_phrase (s->lexer, "ELSE IF"));
3984
3985   if (lex_match_id (s->lexer, "ELSE")
3986       && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
3987     goto error;
3988
3989   if (!lex_match_phrase (s->lexer, "END IF"))
3990     NOT_REACHED ();
3991   return cmd;
3992
3993 error:
3994   matrix_cmd_destroy (cmd);
3995   return NULL;
3996 }
3997
3998 static bool
3999 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4000 {
4001   for (size_t i = 0; i < cmd->n_clauses; i++)
4002     {
4003       struct do_if_clause *c = &cmd->clauses[i];
4004       if (c->condition)
4005         {
4006           double d;
4007           if (!matrix_expr_evaluate_scalar (c->condition,
4008                                             i ? "ELSE IF" : "DO IF",
4009                                             &d) || d <= 0)
4010             continue;
4011         }
4012
4013       for (size_t j = 0; j < c->commands.n; j++)
4014         if (!matrix_cmd_execute (c->commands.commands[j]))
4015           return false;
4016       return true;
4017     }
4018   return true;
4019 }
4020 \f
4021 static struct matrix_cmd *
4022 matrix_parse_loop (struct matrix_state *s)
4023 {
4024   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4025   *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4026
4027   struct loop_command *loop = &cmd->loop;
4028   if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4029     {
4030       struct substring name = lex_tokss (s->lexer);
4031       loop->var = matrix_var_lookup (s, name);
4032       if (!loop->var)
4033         loop->var = matrix_var_create (s, name);
4034
4035       lex_get (s->lexer);
4036       lex_get (s->lexer);
4037
4038       loop->start = matrix_parse_expr (s);
4039       if (!loop->start || !lex_force_match (s->lexer, T_TO))
4040         goto error;
4041
4042       loop->end = matrix_parse_expr (s);
4043       if (!loop->end)
4044         goto error;
4045
4046       if (lex_match (s->lexer, T_BY))
4047         {
4048           loop->increment = matrix_parse_expr (s);
4049           if (!loop->increment)
4050             goto error;
4051         }
4052     }
4053
4054   if (lex_match_id (s->lexer, "IF"))
4055     {
4056       loop->top_condition = matrix_parse_expr (s);
4057       if (!loop->top_condition)
4058         goto error;
4059     }
4060
4061   bool was_in_loop = s->in_loop;
4062   s->in_loop = true;
4063   bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4064                                    "END LOOP", NULL);
4065   s->in_loop = was_in_loop;
4066   if (!ok)
4067     goto error;
4068
4069   if (!lex_match_phrase (s->lexer, "END LOOP"))
4070     NOT_REACHED ();
4071
4072   if (lex_match_id (s->lexer, "IF"))
4073     {
4074       loop->bottom_condition = matrix_parse_expr (s);
4075       if (!loop->bottom_condition)
4076         goto error;
4077     }
4078
4079   return cmd;
4080
4081 error:
4082   matrix_cmd_destroy (cmd);
4083   return NULL;
4084 }
4085
4086 static void
4087 matrix_cmd_execute_loop (struct loop_command *cmd)
4088 {
4089   long int value = 0;
4090   long int end = 0;
4091   long int increment = 1;
4092   if (cmd->var)
4093     {
4094       if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4095           || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4096           || (cmd->increment
4097               && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4098                                                 &increment)))
4099         return;
4100
4101       if (increment > 0 ? value > end
4102           : increment < 0 ? value < end
4103           : true)
4104         return;
4105     }
4106
4107   int mxloops = settings_get_mxloops ();
4108   for (int i = 0; i < mxloops; i++)
4109     {
4110       if (cmd->var)
4111         {
4112           if (cmd->var->value
4113               && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4114             {
4115               gsl_matrix_free (cmd->var->value);
4116               cmd->var->value = NULL;
4117             }
4118           if (!cmd->var->value)
4119             cmd->var->value = gsl_matrix_alloc (1, 1);
4120
4121           gsl_matrix_set (cmd->var->value, 0, 0, value);
4122         }
4123
4124       if (cmd->top_condition)
4125         {
4126           double d;
4127           if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4128                                             &d) || d <= 0)
4129             return;
4130         }
4131
4132       for (size_t j = 0; j < cmd->commands.n; j++)
4133         if (!matrix_cmd_execute (cmd->commands.commands[j]))
4134           return;
4135
4136       if (cmd->bottom_condition)
4137         {
4138           double d;
4139           if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4140                                             "END LOOP IF", &d) || d > 0)
4141             return;
4142         }
4143
4144       if (cmd->var)
4145         {
4146           value += increment;
4147           if (increment > 0 ? value > end : value < end)
4148             return;
4149         }
4150     }
4151 }
4152 \f
4153 static struct matrix_cmd *
4154 matrix_parse_break (struct matrix_state *s)
4155 {
4156   if (!s->in_loop)
4157     {
4158       msg (SE, _("BREAK not inside LOOP."));
4159       return NULL;
4160     }
4161
4162   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4163   *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4164   return cmd;
4165 }
4166 \f
4167 static struct matrix_cmd *
4168 matrix_parse_display (struct matrix_state *s)
4169 {
4170   bool dictionary = false;
4171   bool status = false;
4172   while (lex_token (s->lexer) == T_ID)
4173     {
4174       if (lex_match_id (s->lexer, "DICTIONARY"))
4175         dictionary = true;
4176       else if (lex_match_id (s->lexer, "STATUS"))
4177         status = true;
4178       else
4179         {
4180           lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4181           return NULL;
4182         }
4183     }
4184   if (!dictionary && !status)
4185     dictionary = status = true;
4186
4187   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4188   *cmd = (struct matrix_cmd) {
4189     .type = MCMD_DISPLAY,
4190     .display = {
4191       .state = s,
4192       .dictionary = dictionary,
4193       .status = status,
4194     }
4195   };
4196   return cmd;
4197 }
4198
4199 static int
4200 compare_matrix_var_pointers (const void *a_, const void *b_)
4201 {
4202   const struct matrix_var *const *ap = a_;
4203   const struct matrix_var *const *bp = b_;
4204   const struct matrix_var *a = *ap;
4205   const struct matrix_var *b = *bp;
4206   return strcmp (a->name, b->name);
4207 }
4208
4209 static void
4210 matrix_cmd_execute_display (struct display_command *cmd)
4211 {
4212   const struct matrix_state *s = cmd->state;
4213
4214   struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4215   pivot_dimension_create (
4216     table, PIVOT_AXIS_COLUMN, N_("Property"),
4217     N_("Rows"), N_("Columns"), N_("Size (kB)"));
4218
4219   const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4220   size_t n_vars = 0;
4221
4222   const struct matrix_var *var;
4223   HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4224     if (var->value)
4225       vars[n_vars++] = var;
4226   qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4227
4228   struct pivot_dimension *rows = pivot_dimension_create (
4229     table, PIVOT_AXIS_ROW, N_("Variable"));
4230   for (size_t i = 0; i < n_vars; i++)
4231     {
4232       const struct matrix_var *var = vars[i];
4233       pivot_category_create_leaf (
4234         rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4235
4236       size_t r = var->value->size1;
4237       size_t c = var->value->size2;
4238       double values[] = { r, c, r * c * sizeof (double) / 1024 };
4239       for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4240         pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4241     }
4242   pivot_table_submit (table);
4243 }
4244 \f
4245 static struct matrix_cmd *
4246 matrix_parse_release (struct matrix_state *s)
4247 {
4248   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4249   *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4250
4251   size_t allocated_vars = 0;
4252   while (lex_token (s->lexer) == T_ID)
4253     {
4254       struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4255       if (var)
4256         {
4257           if (cmd->release.n_vars >= allocated_vars)
4258             cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4259                                             sizeof *cmd->release.vars);
4260           cmd->release.vars[cmd->release.n_vars++] = var;
4261         }
4262       else
4263         lex_error (s->lexer, _("Variable name expected."));
4264       lex_get (s->lexer);
4265
4266       if (!lex_match (s->lexer, T_COMMA))
4267         break;
4268     }
4269
4270   return cmd;
4271 }
4272
4273 static void
4274 matrix_cmd_execute_release (struct release_command *cmd)
4275 {
4276   for (size_t i = 0; i < cmd->n_vars; i++)
4277     {
4278       struct matrix_var *var = cmd->vars[i];
4279       gsl_matrix_free (var->value);
4280       var->value = NULL;
4281     }
4282 }
4283 \f
4284 static struct matrix_cmd *
4285 matrix_parse_save (struct matrix_state *s)
4286 {
4287   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4288   *cmd = (struct matrix_cmd) {
4289     .type = MCMD_SAVE,
4290     .save = {
4291       .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4292     }
4293   };
4294
4295   struct save_command *save = &cmd->save;
4296   save->expression = matrix_parse_exp (s);
4297   if (!save->expression)
4298     goto error;
4299
4300   while (lex_match (s->lexer, T_SLASH))
4301     {
4302       if (lex_match_id (s->lexer, "OUTFILE"))
4303         {
4304           lex_match (s->lexer, T_EQUALS);
4305
4306           fh_unref (save->outfile);
4307           save->outfile = (lex_match (s->lexer, T_ASTERISK)
4308                            ? fh_inline_file ()
4309                            : fh_parse (s->lexer, FH_REF_FILE, s->session));
4310           if (!save->outfile)
4311             goto error;
4312         }
4313       else if (lex_match_id (s->lexer, "VARIABLES"))
4314         {
4315           lex_match (s->lexer, T_EQUALS);
4316
4317           char **names;
4318           size_t n;
4319           struct dictionary *d = dict_create (get_default_encoding ());
4320           bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4321                                           PV_NO_SCRATCH | PV_NO_DUPLICATE);
4322           dict_unref (d);
4323           if (!ok)
4324             goto error;
4325
4326           string_array_destroy (save->variables);
4327           if (!save->variables)
4328             save->variables = xmalloc (sizeof *save->variables);
4329           *save->variables = (struct string_array) {
4330             .strings = names,
4331             .n = n,
4332             .allocated = n,
4333           };
4334         }
4335       else if (lex_match_id (s->lexer, "NAMES"))
4336         {
4337           lex_match (s->lexer, T_EQUALS);
4338           matrix_expr_destroy (save->names);
4339           save->names = matrix_parse_exp (s);
4340           if (!save->names)
4341             goto error;
4342         }
4343       else if (lex_match_id (s->lexer, "STRINGS"))
4344         {
4345           lex_match (s->lexer, T_EQUALS);
4346           while (lex_token (s->lexer) == T_ID)
4347             {
4348               stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4349               lex_get (s->lexer);
4350               if (!lex_match (s->lexer, T_COMMA))
4351                 break;
4352             }
4353         }
4354       else
4355         {
4356           lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4357                                "STRINGS");
4358           goto error;
4359         }
4360     }
4361
4362   if (!save->outfile)
4363     {
4364       if (s->prev_save_outfile)
4365         save->outfile = fh_ref (s->prev_save_outfile);
4366       else
4367         {
4368           lex_sbc_missing ("OUTFILE");
4369           goto error;
4370         }
4371     }
4372   fh_unref (s->prev_save_outfile);
4373   s->prev_save_outfile = fh_ref (save->outfile);
4374
4375   if (save->variables && save->names)
4376     {
4377       msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4378       matrix_expr_destroy (save->names);
4379       save->names = NULL;
4380     }
4381
4382   return cmd;
4383
4384 error:
4385   matrix_cmd_destroy (cmd);
4386   return NULL;
4387 }
4388
4389 static void
4390 matrix_cmd_execute_save (const struct save_command *save)
4391 {
4392   assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4393
4394   gsl_matrix *m = matrix_expr_evaluate (save->expression);
4395   if (!m)
4396     return;
4397
4398   bool ok = true;
4399   struct dictionary *dict = dict_create (get_default_encoding ());
4400   struct string_array names = { .n = 0 };
4401   if (save->variables)
4402     string_array_clone (&names, save->variables);
4403   else if (save->names)
4404     {
4405       gsl_matrix *nm = matrix_expr_evaluate (save->names);
4406       if (nm && is_vector (nm))
4407         {
4408           gsl_vector nv = to_vector (nm);
4409           for (size_t i = 0; i < nv.size; i++)
4410             {
4411               char *name = trimmed_string (gsl_vector_get (&nv, i));
4412               if (dict_id_is_valid (dict, name, true))
4413                 string_array_append_nocopy (&names, name);
4414               else
4415                 ok = false;
4416             }
4417         }
4418     }
4419
4420   struct stringi_set strings;
4421   stringi_set_clone (&strings, &save->strings);
4422
4423   for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4424     {
4425       char tmp_name[64];
4426       const char *name;
4427       if (i < names.n)
4428         name = names.strings[i];
4429       else
4430         {
4431           snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4432           name = tmp_name;
4433         }
4434
4435       int width = stringi_set_delete (&strings, name) ? 8 : 0;
4436       struct variable *var = dict_create_var (dict, name, width);
4437       if (!var)
4438         {
4439           msg (SE, _("Duplicate variable name %s in SAVE statement."),
4440                name);
4441           ok = false;
4442         }
4443     }
4444
4445   if (!stringi_set_is_empty (&strings))
4446     {
4447       const char *example = stringi_set_node_get_string (
4448         stringi_set_first (&strings));
4449       msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4450                          "with that name was found on SAVE.",
4451                          "STRINGS specified %2$zu variables, including %1$s, "
4452                          "whose names were not found on SAVE.",
4453                          stringi_set_count (&strings)),
4454            example, stringi_set_count (&strings));
4455       ok = false;
4456     }
4457   stringi_set_destroy (&strings);
4458   string_array_destroy (&names);
4459
4460   if (!ok)
4461     {
4462       gsl_matrix_free (m);
4463       dict_unref (dict);
4464       return;
4465     }
4466
4467   struct casewriter *writer = any_writer_open (save->outfile, dict);
4468   if (!writer)
4469     {
4470       gsl_matrix_free (m);
4471       dict_unref (dict);
4472       return;
4473     }
4474
4475   const struct caseproto *proto = dict_get_proto (dict);
4476   for (size_t y = 0; y < m->size1; y++)
4477     {
4478       struct ccase *c = case_create (proto);
4479       for (size_t x = 0; x < m->size2; x++)
4480         {
4481           double d = gsl_matrix_get (m, y, x);
4482           int width = caseproto_get_width (proto, x);
4483           union value *value = case_data_rw_idx (c, x);
4484           if (width == 0)
4485             value->f = d;
4486           else
4487             memcpy (value->s, &d, width);
4488         }
4489       casewriter_write (writer, c);
4490     }
4491   casewriter_destroy (writer);
4492
4493   gsl_matrix_free (m);
4494   dict_unref (dict);
4495 }
4496 \f
4497 static struct matrix_cmd *
4498 matrix_parse_read (struct matrix_state *s)
4499 {
4500   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4501   *cmd = (struct matrix_cmd) {
4502     .type = MCMD_READ,
4503     .read = { .format = FMT_F },
4504   };
4505
4506   struct file_handle *fh = NULL;
4507   char *encoding = NULL;
4508   struct read_command *read = &cmd->read;
4509   read->dst = matrix_lvalue_parse (s);
4510   if (!read->dst)
4511     goto error;
4512
4513   int by = 0;
4514   int repetitions = 0;
4515   int record_width = 0;
4516   bool seen_format = false;
4517   while (lex_match (s->lexer, T_SLASH))
4518     {
4519       if (lex_match_id (s->lexer, "FILE"))
4520         {
4521           lex_match (s->lexer, T_EQUALS);
4522
4523           fh_unref (fh);
4524           fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4525           if (!fh)
4526             goto error;
4527         }
4528       else if (lex_match_id (s->lexer, "ENCODING"))
4529         {
4530           lex_match (s->lexer, T_EQUALS);
4531           if (!lex_force_string (s->lexer))
4532             goto error;
4533
4534           free (encoding);
4535           encoding = ss_xstrdup (lex_tokss (s->lexer));
4536
4537           lex_get (s->lexer);
4538         }
4539       else if (lex_match_id (s->lexer, "FIELD"))
4540         {
4541           lex_match (s->lexer, T_EQUALS);
4542
4543           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4544             goto error;
4545           read->c1 = lex_integer (s->lexer);
4546           lex_get (s->lexer);
4547           if (!lex_force_match (s->lexer, T_TO)
4548               || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4549             goto error;
4550           read->c2 = lex_integer (s->lexer) + 1;
4551           lex_get (s->lexer);
4552
4553           record_width = read->c2 - read->c1;
4554           if (lex_match (s->lexer, T_BY))
4555             {
4556               if (!lex_force_int_range (s->lexer, "BY", 1,
4557                                         read->c2 - read->c1))
4558                 goto error;
4559               by = lex_integer (s->lexer);
4560               lex_get (s->lexer);
4561
4562               if (record_width % by)
4563                 {
4564                   msg (SE, _("BY %d does not evenly divide record width %d."),
4565                        by, record_width);
4566                   goto error;
4567                 }
4568             }
4569           else
4570             by = 0;
4571         }
4572       else if (lex_match_id (s->lexer, "SIZE"))
4573         {
4574           lex_match (s->lexer, T_EQUALS);
4575           read->size = matrix_parse_exp (s);
4576           if (!read->size)
4577             goto error;
4578         }
4579       else if (lex_match_id (s->lexer, "MODE"))
4580         {
4581           lex_match (s->lexer, T_EQUALS);
4582           if (lex_match_id (s->lexer, "RECTANGULAR"))
4583             read->symmetric = false;
4584           else if (lex_match_id (s->lexer, "SYMMETRIC"))
4585             read->symmetric = true;
4586           else
4587             {
4588               lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4589               goto error;
4590             }
4591         }
4592       else if (lex_match_id (s->lexer, "REREAD"))
4593         read->reread = true;
4594       else if (lex_match_id (s->lexer, "FORMAT"))
4595         {
4596           if (seen_format)
4597             {
4598               lex_sbc_only_once ("FORMAT");
4599               goto error;
4600             }
4601           seen_format = true;
4602
4603           lex_match (s->lexer, T_EQUALS);
4604
4605           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4606             goto error;
4607
4608           const char *p = lex_tokcstr (s->lexer);
4609           if (c_isdigit (p[0]))
4610             {
4611               repetitions = atoi (p);
4612               p += strspn (p, "0123456789");
4613               if (!fmt_from_name (p, &read->format))
4614                 {
4615                   lex_error (s->lexer, _("Unknown format %s."), p);
4616                   goto error;
4617                 }
4618               lex_get (s->lexer);
4619             }
4620           else if (!fmt_from_name (p, &read->format))
4621             {
4622               struct fmt_spec format;
4623               if (!parse_format_specifier (s->lexer, &format))
4624                 goto error;
4625               read->format = format.type;
4626               read->w = format.w;
4627               read->d = format.d;
4628             }
4629         }
4630       else
4631         {
4632           lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4633                                "REREAD", "FORMAT");
4634           goto error;
4635         }
4636     }
4637
4638   if (!read->c1)
4639     {
4640       lex_sbc_missing ("FIELD");
4641       goto error;
4642     }
4643
4644   if (!read->dst->n_indexes && !read->size)
4645     {
4646       msg (SE, _("SIZE is required for reading data into a full matrix "
4647                  "(as opposed to a submatrix)."));
4648       goto error;
4649     }
4650
4651   if (!fh)
4652     {
4653       if (s->prev_read_file)
4654         fh = fh_ref (s->prev_read_file);
4655       else
4656         {
4657           lex_sbc_missing ("FILE");
4658           goto error;
4659         }
4660     }
4661   fh_unref (s->prev_read_file);
4662   s->prev_read_file = fh_ref (fh);
4663
4664   read->rf = read_file_create (s, fh);
4665   if (encoding)
4666     {
4667       free (read->rf->encoding);
4668       read->rf->encoding = encoding;
4669       encoding = NULL;
4670     }
4671
4672   /* Field width may be specified in multiple ways:
4673
4674      1. BY on FIELD.
4675      2. The format on FORMAT.
4676      3. The repetition factor on FORMAT.
4677
4678      (2) and (3) are mutually exclusive.
4679
4680      If more than one of these is present, they must agree.  If none of them is
4681      present, then free-field format is used.
4682    */
4683   if (repetitions > record_width)
4684     {
4685       msg (SE, _("%d repetitions cannot fit in record width %d."),
4686            repetitions, record_width);
4687       goto error;
4688     }
4689   int w = (repetitions ? record_width / repetitions
4690            : read->w ? read->w
4691            : by);
4692   if (by && w != by)
4693     {
4694       if (repetitions)
4695         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4696                    "which implies field width %d, "
4697                    "but BY specifies field width %d."),
4698              repetitions, record_width, w, by);
4699       else
4700         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4701              w, by);
4702       goto error;
4703     }
4704   read->w = w;
4705   return cmd;
4706
4707 error:
4708   matrix_cmd_destroy (cmd);
4709   free (encoding);
4710   return NULL;
4711 }
4712
4713 static void
4714 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4715                        gsl_matrix *m, struct substring p, size_t y, size_t x)
4716 {
4717   const char *input_encoding = dfm_reader_get_encoding (reader);
4718   union value v;
4719   char *error = data_in (p, input_encoding, read->format,
4720                          settings_get_fmt_settings (), &v, 0, NULL);
4721   /* XXX report error if value is missing */
4722   if (error)
4723     msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4724   else
4725     {
4726       gsl_matrix_set (m, y, x, v.f);
4727       if (read->symmetric && x != y)
4728         gsl_matrix_set (m, x, y, v.f);
4729     }
4730 }
4731
4732 static bool
4733 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4734                   struct substring *line)
4735 {
4736   if (dfm_eof (reader))
4737     {
4738       msg (SE, _("Unexpected end of file reading matrix data."));
4739       return false;
4740     }
4741   dfm_expand_tabs (reader);
4742   *line = ss_substr (dfm_get_record (reader),
4743                      read->c1 - 1, read->c2 - read->c1);
4744   return true;
4745 }
4746
4747 static void
4748 matrix_read (struct read_command *read, struct dfm_reader *reader,
4749              gsl_matrix *m)
4750 {
4751   for (size_t y = 0; y < m->size1; y++)
4752     {
4753       size_t nx = read->symmetric ? y + 1 : m->size2;
4754
4755       struct substring line = ss_empty ();
4756       for (size_t x = 0; x < nx; x++)
4757         {
4758           struct substring p;
4759           if (!read->w)
4760             {
4761               for (;;)
4762                 {
4763                   ss_ltrim (&line, ss_cstr (" ,"));
4764                   if (!ss_is_empty (line))
4765                     break;
4766                   if (!matrix_read_line (read, reader, &line))
4767                     return;
4768                   dfm_forward_record (reader);
4769                 }
4770
4771               ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4772             }
4773           else
4774             {
4775               if (!matrix_read_line (read, reader, &line))
4776                 return;
4777               size_t fields_per_line = (read->c2 - read->c1) / read->w;
4778               int f = x % fields_per_line;
4779               if (f == fields_per_line - 1)
4780                 dfm_forward_record (reader);
4781
4782               p = ss_substr (line, read->w * f, read->w);
4783             }
4784
4785           matrix_read_set_field (read, reader, m, p, y, x);
4786         }
4787
4788       if (read->w)
4789         dfm_forward_record (reader);
4790       else
4791         {
4792           ss_ltrim (&line, ss_cstr (" ,"));
4793           if (!ss_is_empty (line))
4794             msg (SW, _("Trailing garbage on line \"%.*s\""),
4795                  (int) line.length, line.string);
4796         }
4797     }
4798 }
4799
4800 static void
4801 matrix_cmd_execute_read (struct read_command *read)
4802 {
4803   struct index_vector iv0, iv1;
4804   if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4805     return;
4806
4807   size_t size[2] = { SIZE_MAX, SIZE_MAX };
4808   if (read->size)
4809     {
4810       gsl_matrix *m = matrix_expr_evaluate (read->size);
4811       if (!m)
4812         return;
4813
4814       if (!is_vector (m))
4815         {
4816           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4817           gsl_matrix_free (m);
4818           free (iv0.indexes);
4819           free (iv1.indexes);
4820           return;
4821         }
4822
4823       gsl_vector v = to_vector (m);
4824       double d[2];
4825       if (v.size == 1)
4826         {
4827           d[0] = gsl_vector_get (&v, 0);
4828           d[1] = 1;
4829         }
4830       else if (v.size == 2)
4831         {
4832           d[0] = gsl_vector_get (&v, 0);
4833           d[1] = gsl_vector_get (&v, 1);
4834         }
4835       else
4836         {
4837           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4838           gsl_matrix_free (m);
4839           free (iv0.indexes);
4840           free (iv1.indexes);
4841           return;
4842         }
4843       gsl_matrix_free (m);
4844
4845       if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4846         {
4847           msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4848           free (iv0.indexes);
4849           free (iv1.indexes);
4850           return;
4851         }
4852
4853       size[0] = d[0];
4854       size[1] = d[1];
4855     }
4856
4857   if (read->dst->n_indexes)
4858     {
4859       size_t submatrix_size[2];
4860       if (read->dst->n_indexes == 2)
4861         {
4862           submatrix_size[0] = iv0.n;
4863           submatrix_size[1] = iv1.n;
4864         }
4865       else if (read->dst->var->value->size1 == 1)
4866         {
4867           submatrix_size[0] = 1;
4868           submatrix_size[1] = iv0.n;
4869         }
4870       else
4871         {
4872           submatrix_size[0] = iv0.n;
4873           submatrix_size[1] = 1;
4874         }
4875
4876       if (read->size)
4877         {
4878           if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4879             {
4880               msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4881                          "(%zu,%zu)."),
4882                    size[0], size[1],
4883                    submatrix_size[0], submatrix_size[1]);
4884               free (iv0.indexes);
4885               free (iv1.indexes);
4886               return;
4887             }
4888         }
4889       else
4890         {
4891           size[0] = submatrix_size[0];
4892           size[1] = submatrix_size[1];
4893         }
4894     }
4895
4896   struct dfm_reader *reader = read_file_open (read->rf);
4897   if (read->reread)
4898     dfm_reread_record (reader, 1);
4899
4900   if (read->symmetric && size[0] != size[1])
4901     {
4902       msg (SE, _("Cannot read matrix with non-square dimensions (%zu,%zu) "
4903                  "using READ with MODE=SYMMETRIC."),
4904            size[0], size[1]);
4905       free (iv0.indexes);
4906       free (iv1.indexes);
4907       return;
4908     }
4909   gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4910   matrix_read (read, reader, tmp);
4911   matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4912 }
4913 \f
4914 static struct matrix_cmd *
4915 matrix_parse_write (struct matrix_state *s)
4916 {
4917   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4918   *cmd = (struct matrix_cmd) {
4919     .type = MCMD_WRITE,
4920     .write = { .format = FMT_F },
4921   };
4922
4923   struct write_command *write = &cmd->write;
4924   write->expression = matrix_parse_exp (s);
4925   if (!write->expression)
4926     goto error;
4927
4928   int by = 0;
4929   int repetitions = 0;
4930   int record_width = 0;
4931   bool seen_format = false;
4932   while (lex_match (s->lexer, T_SLASH))
4933     {
4934       if (lex_match_id (s->lexer, "OUTFILE"))
4935         {
4936           lex_match (s->lexer, T_EQUALS);
4937
4938           fh_unref (write->outfile);
4939           write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
4940           if (!write->outfile)
4941             goto error;
4942         }
4943       else if (lex_match_id (s->lexer, "ENCODING"))
4944         {
4945           lex_match (s->lexer, T_EQUALS);
4946           if (!lex_force_string (s->lexer))
4947             goto error;
4948
4949           free (write->encoding);
4950           write->encoding = ss_xstrdup (lex_tokss (s->lexer));
4951
4952           lex_get (s->lexer);
4953         }
4954       else if (lex_match_id (s->lexer, "FIELD"))
4955         {
4956           lex_match (s->lexer, T_EQUALS);
4957
4958           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4959             goto error;
4960           write->c1 = lex_integer (s->lexer);
4961           lex_get (s->lexer);
4962           if (!lex_force_match (s->lexer, T_TO)
4963               || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
4964             goto error;
4965           write->c2 = lex_integer (s->lexer) + 1;
4966           lex_get (s->lexer);
4967
4968           record_width = write->c2 - write->c1;
4969           if (lex_match (s->lexer, T_BY))
4970             {
4971               if (!lex_force_int_range (s->lexer, "BY", 1,
4972                                         write->c2 - write->c1))
4973                 goto error;
4974               by = lex_integer (s->lexer);
4975               lex_get (s->lexer);
4976
4977               if (record_width % by)
4978                 {
4979                   msg (SE, _("BY %d does not evenly divide record width %d."),
4980                        by, record_width);
4981                   goto error;
4982                 }
4983             }
4984           else
4985             by = 0;
4986         }
4987       else if (lex_match_id (s->lexer, "MODE"))
4988         {
4989           lex_match (s->lexer, T_EQUALS);
4990           if (lex_match_id (s->lexer, "RECTANGULAR"))
4991             write->triangular = false;
4992           else if (lex_match_id (s->lexer, "TRIANGULAR"))
4993             write->triangular = true;
4994           else
4995             {
4996               lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
4997               goto error;
4998             }
4999         }
5000       else if (lex_match_id (s->lexer, "HOLD"))
5001         write->hold = true;
5002       else if (lex_match_id (s->lexer, "FORMAT"))
5003         {
5004           if (seen_format)
5005             {
5006               lex_sbc_only_once ("FORMAT");
5007               goto error;
5008             }
5009           seen_format = true;
5010
5011           lex_match (s->lexer, T_EQUALS);
5012
5013           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5014             goto error;
5015
5016           const char *p = lex_tokcstr (s->lexer);
5017           if (c_isdigit (p[0]))
5018             {
5019               repetitions = atoi (p);
5020               p += strspn (p, "0123456789");
5021               if (!fmt_from_name (p, &write->format))
5022                 {
5023                   lex_error (s->lexer, _("Unknown format %s."), p);
5024                   goto error;
5025                 }
5026               lex_get (s->lexer);
5027             }
5028           else if (!fmt_from_name (p, &write->format))
5029             {
5030               struct fmt_spec format;
5031               if (!parse_format_specifier (s->lexer, &format))
5032                 goto error;
5033               write->format = format.type;
5034               write->w = format.w;
5035               write->d = format.d;
5036             }
5037         }
5038       else
5039         {
5040           lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5041                                "HOLD", "FORMAT");
5042           goto error;
5043         }
5044     }
5045
5046   if (!write->c1)
5047     {
5048       lex_sbc_missing ("FIELD");
5049       goto error;
5050     }
5051
5052   if (!write->outfile)
5053     {
5054       if (s->prev_write_outfile)
5055         write->outfile = fh_ref (s->prev_write_outfile);
5056       else
5057         {
5058           lex_sbc_missing ("OUTFILE");
5059           goto error;
5060         }
5061     }
5062   fh_unref (s->prev_write_outfile);
5063   s->prev_write_outfile = fh_ref (write->outfile);
5064
5065   /* Field width may be specified in multiple ways:
5066
5067      1. BY on FIELD.
5068      2. The format on FORMAT.
5069      3. The repetition factor on FORMAT.
5070
5071      (2) and (3) are mutually exclusive.
5072
5073      If more than one of these is present, they must agree.  If none of them is
5074      present, then free-field format is used.
5075    */
5076   if (repetitions > record_width)
5077     {
5078       msg (SE, _("%d repetitions cannot fit in record width %d."),
5079            repetitions, record_width);
5080       goto error;
5081     }
5082   int w = (repetitions ? record_width / repetitions
5083            : write->w ? write->w
5084            : by);
5085   if (by && w != by)
5086     {
5087       if (repetitions)
5088         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5089                    "which implies field width %d, "
5090                    "but BY specifies field width %d."),
5091              repetitions, record_width, w, by);
5092       else
5093         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5094              w, by);
5095       goto error;
5096     }
5097   write->w = w;
5098   return cmd;
5099
5100 error:
5101   matrix_cmd_destroy (cmd);
5102   return NULL;
5103 }
5104
5105 static void
5106 matrix_cmd_execute_write (struct write_command *write)
5107 {
5108   gsl_matrix *m = matrix_expr_evaluate (write->expression);
5109   if (!m)
5110     return;
5111
5112   if (write->triangular && m->size1 != m->size2)
5113     {
5114       msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5115                  "the matrix to be written has dimensions (%zu,%zu)."),
5116            m->size1, m->size2);
5117       gsl_matrix_free (m);
5118       return;
5119     }
5120
5121   struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
5122   if (!writer)
5123     {
5124       gsl_matrix_free (m);
5125       return;
5126     }
5127
5128   const struct fmt_settings *settings = settings_get_fmt_settings ();
5129   struct fmt_spec format = {
5130     .type = write->format,
5131     .w = write->w ? write->w : 40,
5132     .d = write->d
5133   };
5134   struct u8_line line = U8_LINE_INITIALIZER;
5135   for (size_t y = 0; y < m->size1; y++)
5136     {
5137       size_t nx = write->triangular ? y + 1 : m->size2;
5138       int x0 = write->c1;
5139       for (size_t x = 0; x < nx; x++)
5140         {
5141           /* XXX string values */
5142           union value v = { .f = gsl_matrix_get (m, y, x) };
5143           char *s = (write->w
5144                      ? data_out (&v, NULL, &format, settings)
5145                      : data_out_stretchy (&v, NULL, &format, settings, NULL));
5146           size_t len = strlen (s);
5147           int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5148           if (width + x0 > write->c2)
5149             {
5150               dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5151               u8_line_clear (&line);
5152               x0 = write->c1;
5153             }
5154           u8_line_put (&line, x0, x0 + width, s, len);
5155           free (s);
5156
5157           x0 += write->w ? write->w : width + 1;
5158         }
5159
5160       dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5161       u8_line_clear (&line);
5162     }
5163   u8_line_destroy (&line);
5164   dfm_close_writer (writer);
5165
5166   gsl_matrix_free (m);
5167 }
5168 \f
5169 static struct matrix_cmd *
5170 matrix_parse_get (struct matrix_state *s)
5171 {
5172   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5173   *cmd = (struct matrix_cmd) {
5174     .type = MCMD_GET,
5175     .get = {
5176       .user = { .treatment = MGET_ERROR },
5177       .system = { .treatment = MGET_ERROR },
5178     }
5179   };
5180
5181   struct get_command *get = &cmd->get;
5182   get->dst = matrix_lvalue_parse (s);
5183   if (!get->dst)
5184     goto error;
5185
5186   while (lex_match (s->lexer, T_SLASH))
5187     {
5188       if (lex_match_id (s->lexer, "FILE"))
5189         {
5190           if (get->variables.n)
5191             {
5192               lex_error (s->lexer, _("FILE must precede VARIABLES"));
5193               goto error;
5194             }
5195           lex_match (s->lexer, T_EQUALS);
5196
5197           fh_unref (get->file);
5198           get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5199           if (!get->file)
5200             goto error;
5201         }
5202       else if (lex_match_id (s->lexer, "ENCODING"))
5203         {
5204           if (get->variables.n)
5205             {
5206               lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5207               goto error;
5208             }
5209           lex_match (s->lexer, T_EQUALS);
5210           if (!lex_force_string (s->lexer))
5211             goto error;
5212
5213           free (get->encoding);
5214           get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5215
5216           lex_get (s->lexer);
5217         }
5218       else if (lex_match_id (s->lexer, "VARIABLES"))
5219         {
5220           lex_match (s->lexer, T_EQUALS);
5221
5222           struct dictionary *dict = NULL;
5223           if (!get->file)
5224             {
5225               dict = dataset_dict (s->dataset);
5226               if (dict_get_var_cnt (dict) == 0)
5227                 {
5228                   lex_error (s->lexer, _("GET cannot read empty active file."));
5229                   goto error;
5230                 }
5231             }
5232           else
5233             {
5234               struct casereader *reader = any_reader_open_and_decode (
5235                 get->file, get->encoding, &dict, NULL);
5236               if (!reader)
5237                 goto error;
5238               casereader_destroy (reader);
5239             }
5240
5241           struct variable **vars;
5242           size_t n_vars;
5243           bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5244                                      PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5245           if (!ok)
5246             {
5247               dict_unref (dict);
5248               goto error;
5249             }
5250
5251           string_array_clear (&get->variables);
5252           for (size_t i = 0; i < n_vars; i++)
5253             string_array_append (&get->variables, var_get_name (vars[i]));
5254           free (vars);
5255           dict_unref (dict);
5256         }
5257       else if (lex_match_id (s->lexer, "NAMES"))
5258         {
5259           lex_match (s->lexer, T_EQUALS);
5260           if (!lex_force_id (s->lexer))
5261             goto error;
5262
5263           struct substring name = lex_tokss (s->lexer);
5264           get->names = matrix_var_lookup (s, name);
5265           if (!get->names)
5266             get->names = matrix_var_create (s, name);
5267           lex_get (s->lexer);
5268         }
5269       else if (lex_match_id (s->lexer, "MISSING"))
5270         {
5271           lex_match (s->lexer, T_EQUALS);
5272           if (lex_match_id (s->lexer, "ACCEPT"))
5273             get->user.treatment = MGET_ACCEPT;
5274           else if (lex_match_id (s->lexer, "OMIT"))
5275             get->user.treatment = MGET_OMIT;
5276           else if (lex_is_number (s->lexer))
5277             {
5278               get->user.treatment = MGET_RECODE;
5279               get->user.substitute = lex_number (s->lexer);
5280               lex_get (s->lexer);
5281             }
5282           else
5283             {
5284               lex_error (s->lexer, NULL);
5285               goto error;
5286             }
5287         }
5288       else if (lex_match_id (s->lexer, "SYSMIS"))
5289         {
5290           lex_match (s->lexer, T_EQUALS);
5291           if (lex_match_id (s->lexer, "OMIT"))
5292             get->user.treatment = MGET_OMIT;
5293           else if (lex_is_number (s->lexer))
5294             {
5295               get->user.treatment = MGET_RECODE;
5296               get->user.substitute = lex_number (s->lexer);
5297               lex_get (s->lexer);
5298             }
5299           else
5300             {
5301               lex_error (s->lexer, NULL);
5302               goto error;
5303             }
5304         }
5305       else
5306         {
5307           lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5308                                "MISSING", "SYSMIS");
5309           goto error;
5310         }
5311     }
5312   return cmd;
5313
5314 error:
5315   matrix_cmd_destroy (cmd);
5316   return NULL;
5317 }
5318
5319 static void
5320 matrix_cmd_execute_get (struct get_command *get)
5321 {
5322   assert (get->file);           /* XXX */
5323
5324   struct dictionary *dict;
5325   struct casereader *reader = any_reader_open_and_decode (
5326     get->file, get->encoding, &dict, NULL);
5327   if (!reader)
5328     return;
5329
5330   const struct variable **vars = xnmalloc (
5331     get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5332     sizeof *vars);
5333   size_t n_vars = 0;
5334
5335   if (get->variables.n)
5336     {
5337       for (size_t i = 0; i < get->variables.n; i++)
5338         {
5339           const char *name = get->variables.strings[i];
5340           const struct variable *var = dict_lookup_var (dict, name);
5341           if (!var)
5342             {
5343               msg (SE, _("GET: Data file does not contain variable %s."),
5344                    name);
5345               dict_unref (dict);
5346               free (vars);
5347               return;
5348             }
5349           if (!var_is_numeric (var))
5350             {
5351               msg (SE, _("GET: Variable %s is not numeric."), name);
5352               dict_unref (dict);
5353               free (vars);
5354               return;
5355             }
5356           vars[n_vars++] = var;
5357         }
5358     }
5359   else
5360     {
5361       for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5362         {
5363           const struct variable *var = dict_get_var (dict, i);
5364           if (!var_is_numeric (var))
5365             {
5366               msg (SE, _("GET: Variable %s is not numeric."),
5367                    var_get_name (var));
5368               dict_unref (dict);
5369               free (vars);
5370               return;
5371             }
5372           vars[n_vars++] = var;
5373         }
5374     }
5375
5376   size_t n_rows = 0;
5377   gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5378   long long int casenum = 1;
5379   bool error = false;
5380   for (struct ccase *c = casereader_read (reader); c;
5381        c = casereader_read (reader), casenum++)
5382     {
5383       if (n_rows >= m->size1)
5384         {
5385           gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5386           for (size_t y = 0; y < n_rows; y++)
5387             for (size_t x = 0; x < n_vars; x++)
5388               gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5389           gsl_matrix_free (m);
5390           m = p;
5391         }
5392
5393       bool keep = true;
5394       for (size_t x = 0; x < n_vars; x++)
5395         {
5396           const struct variable *var = vars[x];
5397           double d = case_num (c, var);
5398           if (d == SYSMIS)
5399             {
5400               if (get->system.treatment == MGET_RECODE)
5401                 d = get->system.substitute;
5402               else if (get->system.treatment == MGET_OMIT)
5403                 keep = false;
5404               else
5405                 {
5406                   msg (SE, _("GET: Variable %s in case %lld "
5407                              "is system-missing."),
5408                        var_get_name (var), casenum);
5409                   error = true;
5410                 }
5411             }
5412           else if (var_is_num_missing (var, d, MV_USER))
5413             {
5414               if (get->user.treatment == MGET_RECODE)
5415                 d = get->user.substitute;
5416               else if (get->user.treatment == MGET_OMIT)
5417                 keep = false;
5418               else if (get->user.treatment != MGET_ACCEPT)
5419                 {
5420                   msg (SE, _("GET: Variable %s in case %lld has user-missing "
5421                              "value %g."),
5422                        var_get_name (var), casenum, d);
5423                   error = true;
5424                 }
5425             }
5426           gsl_matrix_set (m, n_rows, x, d);
5427         }
5428       case_unref (c);
5429       if (error)
5430         break;
5431       if (keep)
5432         n_rows++;
5433     }
5434   casereader_destroy (reader);
5435   if (!error)
5436     {
5437       m->size1 = n_rows;
5438       matrix_lvalue_evaluate_and_assign (get->dst, m);
5439     }
5440   else
5441     gsl_matrix_free (m);
5442   dict_unref (dict);
5443   free (vars);
5444 }
5445 \f
5446 static const char *
5447 match_rowtype (struct lexer *lexer)
5448 {
5449   static const char *rowtypes[] = {
5450     "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5451   };
5452   size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5453
5454   for (size_t i = 0; i < n_rowtypes; i++)
5455     if (lex_match_id (lexer, rowtypes[i]))
5456       return rowtypes[i];
5457
5458   lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5459   return NULL;
5460 }
5461
5462 static bool
5463 parse_var_names (struct lexer *lexer, struct string_array *sa)
5464 {
5465   lex_match (lexer, T_EQUALS);
5466
5467   string_array_clear (sa);
5468
5469   struct dictionary *dict = dict_create (get_default_encoding ());
5470   dict_create_var_assert (dict, "ROWTYPE_", 8);
5471   dict_create_var_assert (dict, "VARNAME_", 8);
5472   char **names;
5473   size_t n_names;
5474   bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5475                                   PV_NO_DUPLICATE | PV_NO_SCRATCH);
5476   dict_unref (dict);
5477
5478   if (ok)
5479     {
5480       string_array_clear (sa);
5481       sa->strings = names;
5482       sa->n = sa->allocated = n_names;
5483     }
5484   return ok;
5485 }
5486
5487 static void
5488 msave_common_uninit (struct msave_common *common)
5489 {
5490   if (common)
5491     {
5492       fh_unref (common->outfile);
5493       string_array_destroy (&common->variables);
5494       string_array_destroy (&common->fnames);
5495       string_array_destroy (&common->snames);
5496     }
5497 }
5498
5499 static bool
5500 compare_variables (const char *keyword,
5501                    const struct string_array *new,
5502                    const struct string_array *old)
5503 {
5504   if (new->n)
5505     {
5506       if (!old->n)
5507         {
5508           msg (SE, _("%s may only be specified on MSAVE if it was specified "
5509                      "on the first MSAVE within MATRIX."), keyword);
5510           return false;
5511         }
5512       else if (!string_array_equal_case (old, new))
5513         {
5514           msg (SE, _("%s must specify the same variables each time within "
5515                      "a given MATRIX."), keyword);
5516           return false;
5517         }
5518     }
5519   return true;
5520 }
5521
5522 static struct matrix_cmd *
5523 matrix_parse_msave (struct matrix_state *s)
5524 {
5525   struct msave_common common = { .outfile = NULL };
5526   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5527   *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5528
5529   struct msave_command *msave = &cmd->msave;
5530   if (lex_token (s->lexer) == T_ID)
5531     msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5532   msave->expr = matrix_parse_exp (s);
5533   if (!msave->expr)
5534     return NULL;
5535
5536   while (lex_match (s->lexer, T_SLASH))
5537     {
5538       if (lex_match_id (s->lexer, "TYPE"))
5539         {
5540           lex_match (s->lexer, T_EQUALS);
5541
5542           msave->rowtype = match_rowtype (s->lexer);
5543           if (!msave->rowtype)
5544             goto error;
5545         }
5546       else if (lex_match_id (s->lexer, "OUTFILE"))
5547         {
5548           lex_match (s->lexer, T_EQUALS);
5549
5550           fh_unref (common.outfile);
5551           common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5552           if (!common.outfile)
5553             goto error;
5554         }
5555       else if (lex_match_id (s->lexer, "VARIABLES"))
5556         {
5557           if (!parse_var_names (s->lexer, &common.variables))
5558             goto error;
5559         }
5560       else if (lex_match_id (s->lexer, "FNAMES"))
5561         {
5562           if (!parse_var_names (s->lexer, &common.fnames))
5563             goto error;
5564         }
5565       else if (lex_match_id (s->lexer, "SNAMES"))
5566         {
5567           if (!parse_var_names (s->lexer, &common.snames))
5568             goto error;
5569         }
5570       else if (lex_match_id (s->lexer, "SPLIT"))
5571         {
5572           lex_match (s->lexer, T_EQUALS);
5573
5574           matrix_expr_destroy (msave->splits);
5575           msave->splits = matrix_parse_exp (s);
5576           if (!msave->splits)
5577             goto error;
5578         }
5579       else if (lex_match_id (s->lexer, "FACTOR"))
5580         {
5581           lex_match (s->lexer, T_EQUALS);
5582
5583           matrix_expr_destroy (msave->factors);
5584           msave->factors = matrix_parse_exp (s);
5585           if (!msave->factors)
5586             goto error;
5587         }
5588       else
5589         {
5590           lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5591                                "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5592           goto error;
5593         }
5594     }
5595   if (!msave->rowtype)
5596     {
5597       lex_sbc_missing ("TYPE");
5598       goto error;
5599     }
5600   common.has_splits = msave->splits || common.snames.n;
5601   common.has_factors = msave->factors || common.fnames.n;
5602
5603   struct msave_common *c = s->common ? s->common : &common;
5604   if (c->fnames.n && !msave->factors)
5605     {
5606       msg (SE, _("FNAMES requires FACTOR."));
5607       goto error;
5608     }
5609   if (c->snames.n && !msave->splits)
5610     {
5611       msg (SE, _("SNAMES requires SPLIT."));
5612       goto error;
5613     }
5614   if (c->has_factors && !common.has_factors)
5615     {
5616       msg (SE, _("%s is required because it was present on the first "
5617                  "MSAVE in this MATRIX command."), "FACTOR");
5618       goto error;
5619     }
5620   if (c->has_splits && !common.has_splits)
5621     {
5622       msg (SE, _("%s is required because it was present on the first "
5623                  "MSAVE in this MATRIX command."), "SPLIT");
5624       goto error;
5625     }
5626
5627   if (!s->common)
5628     {
5629       if (!common.outfile)
5630         {
5631           lex_sbc_missing ("OUTFILE");
5632           goto error;
5633         }
5634       s->common = xmemdup (&common, sizeof common);
5635     }
5636   else
5637     {
5638       if (common.outfile)
5639         {
5640           bool same = common.outfile == s->common->outfile;
5641           fh_unref (common.outfile);
5642           if (!same)
5643             {
5644               msg (SE, _("OUTFILE must name the same file on each MSAVE "
5645                          "within a single MATRIX command."));
5646               goto error;
5647             }
5648         }
5649       if (!compare_variables ("VARIABLES",
5650                               &common.variables, &s->common->variables)
5651           || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5652           || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5653         goto error;
5654       msave_common_uninit (&common);
5655     }
5656   msave->common = s->common;
5657   if (!msave->varname_)
5658     msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5659   return cmd;
5660
5661 error:
5662   msave_common_uninit (&common);
5663   matrix_cmd_destroy (cmd);
5664   return NULL;
5665 }
5666
5667 static gsl_vector *
5668 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5669 {
5670   gsl_matrix *m = matrix_expr_evaluate (e);
5671   if (!m)
5672     return NULL;
5673
5674   if (!is_vector (m))
5675     {
5676       msg (SE, _("%s expression must evaluate to vector, not a matrix with "
5677                  "dimensions (%zu,%zu)."),
5678            name, m->size1, m->size2);
5679       gsl_matrix_free (m);
5680       return NULL;
5681     }
5682
5683   return matrix_to_vector (m);
5684 }
5685
5686 static const char *
5687 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5688 {
5689   for (size_t i = 0; i < vars->n; i++)
5690     if (!dict_create_var (d, vars->strings[i], 0))
5691       return vars->strings[i];
5692   return NULL;
5693 }
5694
5695 static struct dictionary *
5696 msave_create_dict (const struct msave_common *common)
5697 {
5698   struct dictionary *dict = dict_create (get_default_encoding ());
5699
5700   const char *dup_split = msave_add_vars (dict, &common->snames);
5701   if (dup_split)
5702     {
5703       msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5704       goto error;
5705     }
5706
5707   dict_create_var_assert (dict, "ROWTYPE_", 8);
5708
5709   const char *dup_factor = msave_add_vars (dict, &common->fnames);
5710   if (dup_factor)
5711     {
5712       msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5713       goto error;
5714     }
5715
5716   dict_create_var_assert (dict, "VARNAME_", 8);
5717
5718   const char *dup_var = msave_add_vars (dict, &common->variables);
5719   if (dup_var)
5720     {
5721       msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5722       goto error;
5723     }
5724
5725   return dict;
5726
5727 error:
5728   dict_unref (dict);
5729   return NULL;
5730 }
5731
5732 static void
5733 matrix_cmd_execute_msave (struct msave_command *msave)
5734 {
5735   struct msave_common *common = msave->common;
5736   gsl_matrix *m = NULL;
5737   gsl_vector *factors = NULL;
5738   gsl_vector *splits = NULL;
5739
5740   m = matrix_expr_evaluate (msave->expr);
5741   if (!m)
5742     goto error;
5743
5744   if (!common->variables.n)
5745     for (size_t i = 0; i < m->size2; i++)
5746       string_array_append_nocopy (&common->variables,
5747                                   xasprintf ("COL%zu", i + 1));
5748
5749   if (m->size2 != common->variables.n)
5750     {
5751       msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5752            m->size2, common->variables.n);
5753       goto error;
5754     }
5755
5756   if (msave->factors)
5757     {
5758       factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5759       if (!factors)
5760         goto error;
5761
5762       if (!common->fnames.n)
5763         for (size_t i = 0; i < factors->size; i++)
5764           string_array_append_nocopy (&common->fnames,
5765                                       xasprintf ("FAC%zu", i + 1));
5766
5767       if (factors->size != common->fnames.n)
5768         {
5769           msg (SE, _("There are %zu factor variables, "
5770                      "but %zu split values were supplied."),
5771                common->fnames.n, factors->size);
5772           goto error;
5773         }
5774     }
5775
5776   if (msave->splits)
5777     {
5778       splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5779       if (!splits)
5780         goto error;
5781
5782       if (!common->fnames.n)
5783         for (size_t i = 0; i < splits->size; i++)
5784           string_array_append_nocopy (&common->fnames,
5785                                       xasprintf ("SPL%zu", i + 1));
5786
5787       if (splits->size != common->fnames.n)
5788         {
5789           msg (SE, _("There are %zu split variables, "
5790                      "but %zu split values were supplied."),
5791                common->fnames.n, splits->size);
5792           goto error;
5793         }
5794     }
5795
5796   if (!common->writer)
5797     {
5798       struct dictionary *dict = msave_create_dict (common);
5799       if (!dict)
5800         goto error;
5801
5802       common->writer = any_writer_open (common->outfile, dict);
5803       if (!common->writer)
5804         {
5805           dict_unref (dict);
5806           goto error;
5807         }
5808
5809       common->dict = dict;
5810     }
5811
5812   for (size_t y = 0; y < m->size1; y++)
5813     {
5814       struct ccase *c = case_create (dict_get_proto (common->dict));
5815       size_t idx = 0;
5816
5817       /* Split variables */
5818       if (splits)
5819         for (size_t i = 0; i < splits->size; i++)
5820           case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5821
5822       /* ROWTYPE_. */
5823       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5824                          msave->rowtype, ' ');
5825
5826       /* Factors. */
5827       if (factors)
5828         for (size_t i = 0; i < factors->size; i++)
5829           *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5830
5831       /* VARNAME_. */
5832       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5833                          msave->varname_, ' ');
5834
5835       /* Continuous variables. */
5836       for (size_t x = 0; x < m->size2; x++)
5837         case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5838       casewriter_write (common->writer, c);
5839     }
5840
5841 error:
5842   gsl_matrix_free (m);
5843   gsl_vector_free (factors);
5844   gsl_vector_free (splits);
5845 }
5846 \f
5847 static struct matrix_cmd *
5848 matrix_parse_mget (struct matrix_state *s)
5849 {
5850   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5851   *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5852
5853   struct mget_command *mget = &cmd->mget;
5854
5855   for (;;)
5856     {
5857       if (lex_match_id (s->lexer, "FILE"))
5858         {
5859           lex_match (s->lexer, T_EQUALS);
5860
5861           fh_unref (mget->file);
5862           mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5863           if (!mget->file)
5864             goto error;
5865         }
5866       else if (lex_match_id (s->lexer, "TYPE"))
5867         {
5868           lex_match (s->lexer, T_EQUALS);
5869           while (lex_token (s->lexer) != T_SLASH
5870                  && lex_token (s->lexer) != T_ENDCMD)
5871             {
5872               const char *rowtype = match_rowtype (s->lexer);
5873               if (!rowtype)
5874                 goto error;
5875
5876               stringi_set_insert (&mget->rowtypes, rowtype);
5877             }
5878         }
5879       else
5880         {
5881           lex_error_expecting (s->lexer, "FILE", "TYPE");
5882           goto error;
5883         }
5884       if (lex_token (s->lexer) == T_ENDCMD)
5885         break;
5886
5887       if (!lex_force_match (s->lexer, T_SLASH))
5888         goto error;
5889     }
5890   return cmd;
5891
5892 error:
5893   matrix_cmd_destroy (cmd);
5894   return NULL;
5895 }
5896
5897 static const struct variable *
5898 get_a8_var (const struct dictionary *d, const char *name)
5899 {
5900   const struct variable *v = dict_lookup_var (d, name);
5901   if (!v)
5902     {
5903       msg (SE, _("Matrix data file lacks %s variable."), name);
5904       return NULL;
5905     }
5906   if (var_get_width (v) != 8)
5907     {
5908       msg (SE, _("%s variable in matrix data file must be "
5909                  "8-byte string, but it has width %d."),
5910            name, var_get_width (v));
5911       return NULL;
5912     }
5913   return v;
5914 }
5915
5916 static bool
5917 is_rowtype (const union value *v, const char *rowtype)
5918 {
5919   struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
5920   ss_rtrim (&vs, ss_cstr (" "));
5921   return ss_equals_case (vs, ss_cstr (rowtype));
5922 }
5923
5924 static void
5925 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
5926                         const struct dictionary *d,
5927                         const struct variable *rowtype_var,
5928                         struct matrix_state *s, size_t si, size_t fi,
5929                         size_t cs, size_t cn)
5930 {
5931   if (!n_rows)
5932     return;
5933
5934   const union value *rowtype_ = case_data (rows[0], rowtype_var);
5935   const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
5936                              : is_rowtype (rowtype_, "CORR") ? "CR"
5937                              : is_rowtype (rowtype_, "MEAN") ? "MN"
5938                              : is_rowtype (rowtype_, "STDDEV") ? "SD"
5939                              : is_rowtype (rowtype_, "N") ? "NC"
5940                              : "CN");
5941
5942   struct string name = DS_EMPTY_INITIALIZER;
5943   ds_put_cstr (&name, name_prefix);
5944   if (fi > 0)
5945     ds_put_format (&name, "F%zu", fi);
5946   if (si > 0)
5947     ds_put_format (&name, "S%zu", si);
5948
5949   struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
5950   if (!mv)
5951     mv = matrix_var_create (s, ds_ss (&name));
5952   else if (mv->value)
5953     {
5954       msg (SW, _("Matrix data file contains variable with existing name %s."),
5955            ds_cstr (&name));
5956       goto exit;
5957     }
5958
5959   gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
5960   size_t n_missing = 0;
5961   for (size_t y = 0; y < n_rows; y++)
5962     {
5963       for (size_t x = 0; x < cn; x++)
5964         {
5965           struct variable *var = dict_get_var (d, cs + x);
5966           double value = case_num (rows[y], var);
5967           if (var_is_num_missing (var, value, MV_ANY))
5968             {
5969               n_missing++;
5970               value = 0.0;
5971             }
5972           gsl_matrix_set (m, y, x, value);
5973         }
5974     }
5975
5976   if (n_missing)
5977     msg (SE, ngettext ("Matrix data file variable %s contains a missing "
5978                        "value, which was treated as zero.",
5979                        "Matrix data file variable %s contains %zu missing "
5980                        "values, which were treated as zero.", n_missing),
5981          ds_cstr (&name), n_missing);
5982   mv->value = m;
5983
5984 exit:
5985   ds_destroy (&name);
5986   for (size_t y = 0; y < n_rows; y++)
5987     case_unref (rows[y]);
5988 }
5989
5990 static bool
5991 var_changed (const struct ccase *ca, const struct ccase *cb,
5992              const struct variable *var)
5993 {
5994   return !value_equal (case_data (ca, var), case_data (cb, var),
5995                        var_get_width (var));
5996 }
5997
5998 static bool
5999 vars_changed (const struct ccase *ca, const struct ccase *cb,
6000               const struct dictionary *d,
6001               size_t first_var, size_t n_vars)
6002 {
6003   for (size_t i = 0; i < n_vars; i++)
6004     {
6005       const struct variable *v = dict_get_var (d, first_var + i);
6006       if (var_changed (ca, cb, v))
6007         return true;
6008     }
6009   return false;
6010 }
6011
6012 static void
6013 matrix_cmd_execute_mget (struct mget_command *mget)
6014 {
6015   struct dictionary *d;
6016   struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6017                                                      &d, NULL);
6018   if (!r)
6019     return;
6020
6021   const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6022   const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6023   if (!rowtype_ || !varname_)
6024     goto exit;
6025
6026   if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6027     {
6028       msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6029       goto exit;
6030     }
6031   if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6032     {
6033       msg (SE, _("Matrix data file contains no continuous variables."));
6034       goto exit;
6035     }
6036
6037   for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6038     {
6039       const struct variable *v = dict_get_var (d, i);
6040       if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6041         {
6042           msg (SE,
6043                _("Matrix data file contains unexpected string variable %s."),
6044                var_get_name (v));
6045           goto exit;
6046         }
6047     }
6048
6049   /* SPLIT variables. */
6050   size_t ss = 0;
6051   size_t sn = var_get_dict_index (rowtype_);
6052   struct ccase *sc = NULL;
6053   size_t si = 0;
6054
6055   /* FACTOR variables. */
6056   size_t fs = var_get_dict_index (rowtype_) + 1;
6057   size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6058   struct ccase *fc = NULL;
6059   size_t fi = 0;
6060
6061   /* Continuous variables. */
6062   size_t cs = var_get_dict_index (varname_) + 1;
6063   size_t cn = dict_get_var_cnt (d) - cs;
6064   struct ccase *cc = NULL;
6065
6066   /* Matrix. */
6067   struct ccase **rows = NULL;
6068   size_t allocated_rows = 0;
6069   size_t n_rows = 0;
6070
6071   struct ccase *c;
6072   while ((c = casereader_read (r)) != NULL)
6073     {
6074       bool sd = vars_changed (sc, c, d, ss, sn);
6075       bool fd = sd || vars_changed (fc, c, d, fs, fn);
6076       bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6077       if (sd)
6078         {
6079           si++;
6080           case_unref (sc);
6081           sc = case_ref (c);
6082         }
6083       if (fd)
6084         {
6085           fi++;
6086           case_unref (fc);
6087           fc = case_ref (c);
6088         }
6089       if (md)
6090         {
6091           matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6092                                   mget->state, si, fi, cs, cn);
6093           n_rows = 0;
6094           case_unref (cc);
6095           cc = case_ref (c);
6096         }
6097
6098       if (n_rows >= allocated_rows)
6099         rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6100       rows[n_rows++] = c;
6101     }
6102   matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6103                           mget->state, si, fi, cs, cn);
6104   free (rows);
6105
6106 exit:
6107   casereader_destroy (r);
6108 }
6109 \f
6110 static bool
6111 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6112 {
6113   if (!lex_force_id (s->lexer))
6114     return false;
6115
6116   *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6117   if (!*varp)
6118     *varp = matrix_var_create (s, lex_tokss (s->lexer));
6119   lex_get (s->lexer);
6120   return true;
6121 }
6122
6123 static struct matrix_cmd *
6124 matrix_parse_eigen (struct matrix_state *s)
6125 {
6126   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6127   *cmd = (struct matrix_cmd) {
6128     .type = MCMD_EIGEN,
6129     .eigen = { .expr = NULL }
6130   };
6131
6132   struct eigen_command *eigen = &cmd->eigen;
6133   if (!lex_force_match (s->lexer, T_LPAREN))
6134     goto error;
6135   eigen->expr = matrix_parse_expr (s);
6136   if (!eigen->expr
6137       || !lex_force_match (s->lexer, T_COMMA)
6138       || !matrix_parse_dst_var (s, &eigen->evec)
6139       || !lex_force_match (s->lexer, T_COMMA)
6140       || !matrix_parse_dst_var (s, &eigen->eval)
6141       || !lex_force_match (s->lexer, T_RPAREN))
6142     goto error;
6143
6144   return cmd;
6145
6146 error:
6147   matrix_cmd_destroy (cmd);
6148   return NULL;
6149 }
6150
6151 static void
6152 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6153 {
6154   gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6155   if (!A)
6156     return;
6157   if (!is_symmetric (A))
6158     {
6159       msg (SE, _("Argument of EIGEN must be symmetric."));
6160       gsl_matrix_free (A);
6161       return;
6162     }
6163
6164   gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6165   gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6166   gsl_vector v_eval = to_vector (eval);
6167   gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6168   gsl_eigen_symmv (A, &v_eval, evec, w);
6169   gsl_eigen_symmv_free (w);
6170
6171   gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6172
6173   gsl_matrix_free (A);
6174
6175   gsl_matrix_free (eigen->eval->value);
6176   eigen->eval->value = eval;
6177
6178   gsl_matrix_free (eigen->evec->value);
6179   eigen->evec->value = evec;
6180 }
6181 \f
6182 static struct matrix_cmd *
6183 matrix_parse_setdiag (struct matrix_state *s)
6184 {
6185   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6186   *cmd = (struct matrix_cmd) {
6187     .type = MCMD_SETDIAG,
6188     .setdiag = { .dst = NULL }
6189   };
6190
6191   struct setdiag_command *setdiag = &cmd->setdiag;
6192   if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6193     goto error;
6194
6195   setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6196   if (!setdiag->dst)
6197     {
6198       lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6199       goto error;
6200     }
6201   lex_get (s->lexer);
6202
6203   if (!lex_force_match (s->lexer, T_COMMA))
6204     goto error;
6205
6206   setdiag->expr = matrix_parse_expr (s);
6207   if (!setdiag->expr)
6208     goto error;
6209
6210   if (!lex_force_match (s->lexer, T_RPAREN))
6211     goto error;
6212
6213   return cmd;
6214
6215 error:
6216   matrix_cmd_destroy (cmd);
6217   return NULL;
6218 }
6219
6220 static void
6221 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6222 {
6223   gsl_matrix *dst = setdiag->dst->value;
6224   if (!dst)
6225     {
6226       msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6227            setdiag->dst->name);
6228       return;
6229     }
6230
6231   gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6232   if (!src)
6233     return;
6234
6235   size_t n = MIN (dst->size1, dst->size2);
6236   if (is_scalar (src))
6237     {
6238       double d = to_scalar (src);
6239       for (size_t i = 0; i < n; i++)
6240         gsl_matrix_set (dst, i, i, d);
6241     }
6242   else if (is_vector (src))
6243     {
6244       gsl_vector v = to_vector (src);
6245       for (size_t i = 0; i < n && i < v.size; i++)
6246         gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6247     }
6248   else
6249     msg (SE, _("SETDIAG argument 2 must be a scalar or a vector but it "
6250                "has dimensions (%zu,%zu)."),
6251          src->size1, src->size2);
6252   gsl_matrix_free (src);
6253 }
6254 \f
6255 static struct matrix_cmd *
6256 matrix_parse_svd (struct matrix_state *s)
6257 {
6258   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6259   *cmd = (struct matrix_cmd) {
6260     .type = MCMD_SVD,
6261     .svd = { .expr = NULL }
6262   };
6263
6264   struct svd_command *svd = &cmd->svd;
6265   if (!lex_force_match (s->lexer, T_LPAREN))
6266     goto error;
6267   svd->expr = matrix_parse_expr (s);
6268   if (!svd->expr
6269       || !lex_force_match (s->lexer, T_COMMA)
6270       || !matrix_parse_dst_var (s, &svd->u)
6271       || !lex_force_match (s->lexer, T_COMMA)
6272       || !matrix_parse_dst_var (s, &svd->s)
6273       || !lex_force_match (s->lexer, T_COMMA)
6274       || !matrix_parse_dst_var (s, &svd->v)
6275       || !lex_force_match (s->lexer, T_RPAREN))
6276     goto error;
6277
6278   return cmd;
6279
6280 error:
6281   matrix_cmd_destroy (cmd);
6282   return NULL;
6283 }
6284
6285 static void
6286 matrix_cmd_execute_svd (struct svd_command *svd)
6287 {
6288   gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6289   if (!m)
6290     return;
6291
6292   if (m->size1 >= m->size2)
6293     {
6294       gsl_matrix *A = m;
6295       gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6296       gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6297       gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6298       gsl_vector *work = gsl_vector_alloc (A->size2);
6299       gsl_linalg_SV_decomp (A, V, &Sv, work);
6300       gsl_vector_free (work);
6301
6302       matrix_var_set (svd->u, A);
6303       matrix_var_set (svd->s, S);
6304       matrix_var_set (svd->v, V);
6305     }
6306   else
6307     {
6308       gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6309       gsl_matrix_transpose_memcpy (At, m);
6310       gsl_matrix_free (m);
6311
6312       gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6313       gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6314       gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6315       gsl_vector *work = gsl_vector_alloc (At->size2);
6316       gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6317       gsl_vector_free (work);
6318
6319       matrix_var_set (svd->v, At);
6320       matrix_var_set (svd->s, St);
6321       matrix_var_set (svd->u, Vt);
6322     }
6323 }
6324 \f
6325 static bool
6326 matrix_cmd_execute (struct matrix_cmd *cmd)
6327 {
6328   switch (cmd->type)
6329     {
6330     case MCMD_COMPUTE:
6331       matrix_cmd_execute_compute (&cmd->compute);
6332       break;
6333
6334     case MCMD_PRINT:
6335       matrix_cmd_execute_print (&cmd->print);
6336       break;
6337
6338     case MCMD_DO_IF:
6339       return matrix_cmd_execute_do_if (&cmd->do_if);
6340
6341     case MCMD_LOOP:
6342       matrix_cmd_execute_loop (&cmd->loop);
6343       break;
6344
6345     case MCMD_BREAK:
6346       return false;
6347
6348     case MCMD_DISPLAY:
6349       matrix_cmd_execute_display (&cmd->display);
6350       break;
6351
6352     case MCMD_RELEASE:
6353       matrix_cmd_execute_release (&cmd->release);
6354       break;
6355
6356     case MCMD_SAVE:
6357       matrix_cmd_execute_save (&cmd->save);
6358       break;
6359
6360     case MCMD_READ:
6361       matrix_cmd_execute_read (&cmd->read);
6362       break;
6363
6364     case MCMD_WRITE:
6365       matrix_cmd_execute_write (&cmd->write);
6366       break;
6367
6368     case MCMD_GET:
6369       matrix_cmd_execute_get (&cmd->get);
6370       break;
6371
6372     case MCMD_MSAVE:
6373       matrix_cmd_execute_msave (&cmd->msave);
6374       break;
6375
6376     case MCMD_MGET:
6377       matrix_cmd_execute_mget (&cmd->mget);
6378       break;
6379
6380     case MCMD_EIGEN:
6381       matrix_cmd_execute_eigen (&cmd->eigen);
6382       break;
6383
6384     case MCMD_SETDIAG:
6385       matrix_cmd_execute_setdiag (&cmd->setdiag);
6386       break;
6387
6388     case MCMD_SVD:
6389       matrix_cmd_execute_svd (&cmd->svd);
6390       break;
6391     }
6392
6393   return true;
6394 }
6395
6396 static void
6397 matrix_cmds_uninit (struct matrix_cmds *cmds)
6398 {
6399   for (size_t i = 0; i < cmds->n; i++)
6400     matrix_cmd_destroy (cmds->commands[i]);
6401   free (cmds->commands);
6402 }
6403
6404 static void
6405 matrix_cmd_destroy (struct matrix_cmd *cmd)
6406 {
6407   if (!cmd)
6408     return;
6409
6410   switch (cmd->type)
6411     {
6412     case MCMD_COMPUTE:
6413       matrix_lvalue_destroy (cmd->compute.lvalue);
6414       matrix_expr_destroy (cmd->compute.rvalue);
6415       break;
6416
6417     case MCMD_PRINT:
6418       matrix_expr_destroy (cmd->print.expression);
6419       free (cmd->print.title);
6420       print_labels_destroy (cmd->print.rlabels);
6421       print_labels_destroy (cmd->print.clabels);
6422       break;
6423
6424     case MCMD_DO_IF:
6425       for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6426         {
6427           matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6428           matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6429         }
6430       free (cmd->do_if.clauses);
6431       break;
6432
6433     case MCMD_LOOP:
6434       matrix_expr_destroy (cmd->loop.start);
6435       matrix_expr_destroy (cmd->loop.end);
6436       matrix_expr_destroy (cmd->loop.increment);
6437       matrix_expr_destroy (cmd->loop.top_condition);
6438       matrix_expr_destroy (cmd->loop.bottom_condition);
6439       matrix_cmds_uninit (&cmd->loop.commands);
6440       break;
6441
6442     case MCMD_BREAK:
6443       break;
6444
6445     case MCMD_DISPLAY:
6446       break;
6447
6448     case MCMD_RELEASE:
6449       free (cmd->release.vars);
6450       break;
6451
6452     case MCMD_SAVE:
6453       matrix_expr_destroy (cmd->save.expression);
6454       fh_unref (cmd->save.outfile);
6455       string_array_destroy (cmd->save.variables);
6456       matrix_expr_destroy (cmd->save.names);
6457       stringi_set_destroy (&cmd->save.strings);
6458       break;
6459
6460     case MCMD_READ:
6461       matrix_lvalue_destroy (cmd->read.dst);
6462       matrix_expr_destroy (cmd->read.size);
6463       break;
6464
6465     case MCMD_WRITE:
6466       matrix_expr_destroy (cmd->write.expression);
6467       free (cmd->write.encoding);
6468       break;
6469
6470     case MCMD_GET:
6471       matrix_lvalue_destroy (cmd->get.dst);
6472       fh_unref (cmd->get.file);
6473       free (cmd->get.encoding);
6474       string_array_destroy (&cmd->get.variables);
6475       break;
6476
6477     case MCMD_MSAVE:
6478       free (cmd->msave.varname_);
6479       matrix_expr_destroy (cmd->msave.expr);
6480       matrix_expr_destroy (cmd->msave.factors);
6481       matrix_expr_destroy (cmd->msave.splits);
6482       break;
6483
6484     case MCMD_MGET:
6485       fh_unref (cmd->mget.file);
6486       stringi_set_destroy (&cmd->mget.rowtypes);
6487       break;
6488
6489     case MCMD_EIGEN:
6490       matrix_expr_destroy (cmd->eigen.expr);
6491       break;
6492
6493     case MCMD_SETDIAG:
6494       matrix_expr_destroy (cmd->setdiag.expr);
6495       break;
6496
6497     case MCMD_SVD:
6498       matrix_expr_destroy (cmd->svd.expr);
6499       break;
6500     }
6501   free (cmd);
6502 }
6503
6504 struct matrix_command_name
6505   {
6506     const char *name;
6507     struct matrix_cmd *(*parse) (struct matrix_state *);
6508   };
6509
6510 static const struct matrix_command_name *
6511 matrix_parse_command_name (struct lexer *lexer)
6512 {
6513   static const struct matrix_command_name commands[] = {
6514     { "COMPUTE", matrix_parse_compute },
6515     { "CALL EIGEN", matrix_parse_eigen },
6516     { "CALL SETDIAG", matrix_parse_setdiag },
6517     { "CALL SVD", matrix_parse_svd },
6518     { "PRINT", matrix_parse_print },
6519     { "DO IF", matrix_parse_do_if },
6520     { "LOOP", matrix_parse_loop },
6521     { "BREAK", matrix_parse_break },
6522     { "READ", matrix_parse_read },
6523     { "WRITE", matrix_parse_write },
6524     { "GET", matrix_parse_get },
6525     { "SAVE", matrix_parse_save },
6526     { "MGET", matrix_parse_mget },
6527     { "MSAVE", matrix_parse_msave },
6528     { "DISPLAY", matrix_parse_display },
6529     { "RELEASE", matrix_parse_release },
6530   };
6531   static size_t n = sizeof commands / sizeof *commands;
6532
6533   for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6534     if (lex_match_phrase (lexer, c->name))
6535       return c;
6536   return NULL;
6537 }
6538
6539 static struct matrix_cmd *
6540 matrix_parse_command (struct matrix_state *s)
6541 {
6542   size_t nesting_level = SIZE_MAX;
6543
6544   struct matrix_cmd *c = NULL;
6545   const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6546   if (!cmd)
6547     lex_error (s->lexer, _("Unknown matrix command."));
6548   else if (!cmd->parse)
6549     lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6550                cmd->name);
6551   else
6552     {
6553       nesting_level = output_open_group (
6554         group_item_create_nocopy (utf8_to_title (cmd->name),
6555                                   utf8_to_title (cmd->name)));
6556       c = cmd->parse (s);
6557     }
6558
6559   if (c)
6560     lex_end_of_command (s->lexer);
6561   lex_discard_rest_of_command (s->lexer);
6562   if (nesting_level != SIZE_MAX)
6563     output_close_groups (nesting_level);
6564
6565   return c;
6566 }
6567
6568 int
6569 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6570 {
6571   if (!lex_force_match (lexer, T_ENDCMD))
6572     return CMD_FAILURE;
6573
6574   struct matrix_state state = {
6575     .session = dataset_session (ds),
6576     .lexer = lexer,
6577     .vars = HMAP_INITIALIZER (state.vars),
6578   };
6579
6580   for (;;)
6581     {
6582       while (lex_match (lexer, T_ENDCMD))
6583         continue;
6584       if (lex_token (lexer) == T_STOP)
6585         {
6586           msg (SE, _("Unexpected end of input expecting matrix command."));
6587           break;
6588         }
6589
6590       if (lex_match_phrase (lexer, "END MATRIX"))
6591         break;
6592
6593       struct matrix_cmd *c = matrix_parse_command (&state);
6594       if (c)
6595         {
6596           matrix_cmd_execute (c);
6597           matrix_cmd_destroy (c);
6598         }
6599     }
6600
6601   struct matrix_var *var, *next;
6602   HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6603     {
6604       free (var->name);
6605       gsl_matrix_free (var->value);
6606       hmap_delete (&state.vars, &var->hmap_node);
6607       free (var);
6608     }
6609   hmap_destroy (&state.vars);
6610   fh_unref (state.prev_save_outfile);
6611   fh_unref (state.prev_write_outfile);
6612   if (state.common)
6613     {
6614       dict_unref (state.common->dict);
6615       casewriter_destroy (state.common->writer);
6616       free (state.common);
6617     }
6618   fh_unref (state.prev_read_file);
6619   for (size_t i = 0; i < state.n_read_files; i++)
6620     read_file_destroy (state.read_files[i]);
6621   free (state.read_files);
6622
6623   return CMD_SUCCESS;
6624 }