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