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