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 (fh_equal (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         }
4522       else
4523         {
4524           casewriter_destroy (sf->writer);
4525           dict_unref (sf->dict);
4526         }
4527       fh_unref (sf->file);
4528       string_array_destroy (&sf->variables);
4529       matrix_expr_destroy (sf->names);
4530       stringi_set_destroy (&sf->strings);
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   return cmd;
4637
4638 error:
4639   string_array_destroy (&variables);
4640   matrix_expr_destroy (names);
4641   stringi_set_destroy (&strings);
4642   fh_unref (fh);
4643   matrix_cmd_destroy (cmd);
4644   return NULL;
4645 }
4646
4647 static void
4648 matrix_cmd_execute_save (const struct save_command *save)
4649 {
4650   gsl_matrix *m = matrix_expr_evaluate (save->expression);
4651   if (!m)
4652     return;
4653
4654   struct casewriter *writer = save_file_open (save->sf, m);
4655   if (!writer)
4656     {
4657       gsl_matrix_free (m);
4658       return;
4659     }
4660
4661   const struct caseproto *proto = casewriter_get_proto (writer);
4662   for (size_t y = 0; y < m->size1; y++)
4663     {
4664       struct ccase *c = case_create (proto);
4665       for (size_t x = 0; x < m->size2; x++)
4666         {
4667           double d = gsl_matrix_get (m, y, x);
4668           int width = caseproto_get_width (proto, x);
4669           union value *value = case_data_rw_idx (c, x);
4670           if (width == 0)
4671             value->f = d;
4672           else
4673             memcpy (value->s, &d, width);
4674         }
4675       casewriter_write (writer, c);
4676     }
4677   gsl_matrix_free (m);
4678 }
4679 \f
4680 static struct matrix_cmd *
4681 matrix_parse_read (struct matrix_state *s)
4682 {
4683   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4684   *cmd = (struct matrix_cmd) {
4685     .type = MCMD_READ,
4686     .read = { .format = FMT_F },
4687   };
4688
4689   struct file_handle *fh = NULL;
4690   char *encoding = NULL;
4691   struct read_command *read = &cmd->read;
4692   read->dst = matrix_lvalue_parse (s);
4693   if (!read->dst)
4694     goto error;
4695
4696   int by = 0;
4697   int repetitions = 0;
4698   int record_width = 0;
4699   bool seen_format = false;
4700   while (lex_match (s->lexer, T_SLASH))
4701     {
4702       if (lex_match_id (s->lexer, "FILE"))
4703         {
4704           lex_match (s->lexer, T_EQUALS);
4705
4706           fh_unref (fh);
4707           fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4708           if (!fh)
4709             goto error;
4710         }
4711       else if (lex_match_id (s->lexer, "ENCODING"))
4712         {
4713           lex_match (s->lexer, T_EQUALS);
4714           if (!lex_force_string (s->lexer))
4715             goto error;
4716
4717           free (encoding);
4718           encoding = ss_xstrdup (lex_tokss (s->lexer));
4719
4720           lex_get (s->lexer);
4721         }
4722       else if (lex_match_id (s->lexer, "FIELD"))
4723         {
4724           lex_match (s->lexer, T_EQUALS);
4725
4726           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4727             goto error;
4728           read->c1 = lex_integer (s->lexer);
4729           lex_get (s->lexer);
4730           if (!lex_force_match (s->lexer, T_TO)
4731               || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4732             goto error;
4733           read->c2 = lex_integer (s->lexer) + 1;
4734           lex_get (s->lexer);
4735
4736           record_width = read->c2 - read->c1;
4737           if (lex_match (s->lexer, T_BY))
4738             {
4739               if (!lex_force_int_range (s->lexer, "BY", 1,
4740                                         read->c2 - read->c1))
4741                 goto error;
4742               by = lex_integer (s->lexer);
4743               lex_get (s->lexer);
4744
4745               if (record_width % by)
4746                 {
4747                   msg (SE, _("BY %d does not evenly divide record width %d."),
4748                        by, record_width);
4749                   goto error;
4750                 }
4751             }
4752           else
4753             by = 0;
4754         }
4755       else if (lex_match_id (s->lexer, "SIZE"))
4756         {
4757           lex_match (s->lexer, T_EQUALS);
4758           matrix_expr_destroy (read->size);
4759           read->size = matrix_parse_exp (s);
4760           if (!read->size)
4761             goto error;
4762         }
4763       else if (lex_match_id (s->lexer, "MODE"))
4764         {
4765           lex_match (s->lexer, T_EQUALS);
4766           if (lex_match_id (s->lexer, "RECTANGULAR"))
4767             read->symmetric = false;
4768           else if (lex_match_id (s->lexer, "SYMMETRIC"))
4769             read->symmetric = true;
4770           else
4771             {
4772               lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4773               goto error;
4774             }
4775         }
4776       else if (lex_match_id (s->lexer, "REREAD"))
4777         read->reread = true;
4778       else if (lex_match_id (s->lexer, "FORMAT"))
4779         {
4780           if (seen_format)
4781             {
4782               lex_sbc_only_once ("FORMAT");
4783               goto error;
4784             }
4785           seen_format = true;
4786
4787           lex_match (s->lexer, T_EQUALS);
4788
4789           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4790             goto error;
4791
4792           const char *p = lex_tokcstr (s->lexer);
4793           if (c_isdigit (p[0]))
4794             {
4795               repetitions = atoi (p);
4796               p += strspn (p, "0123456789");
4797               if (!fmt_from_name (p, &read->format))
4798                 {
4799                   lex_error (s->lexer, _("Unknown format %s."), p);
4800                   goto error;
4801                 }
4802               lex_get (s->lexer);
4803             }
4804           else if (fmt_from_name (p, &read->format))
4805             lex_get (s->lexer);
4806           else
4807             {
4808               struct fmt_spec format;
4809               if (!parse_format_specifier (s->lexer, &format))
4810                 goto error;
4811               read->format = format.type;
4812               read->w = format.w;
4813             }
4814         }
4815       else
4816         {
4817           lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4818                                "REREAD", "FORMAT");
4819           goto error;
4820         }
4821     }
4822
4823   if (!read->c1)
4824     {
4825       lex_sbc_missing ("FIELD");
4826       goto error;
4827     }
4828
4829   if (!read->dst->n_indexes && !read->size)
4830     {
4831       msg (SE, _("SIZE is required for reading data into a full matrix "
4832                  "(as opposed to a submatrix)."));
4833       goto error;
4834     }
4835
4836   if (!fh)
4837     {
4838       if (s->prev_read_file)
4839         fh = fh_ref (s->prev_read_file);
4840       else
4841         {
4842           lex_sbc_missing ("FILE");
4843           goto error;
4844         }
4845     }
4846   fh_unref (s->prev_read_file);
4847   s->prev_read_file = fh_ref (fh);
4848
4849   read->rf = read_file_create (s, fh);
4850   fh = NULL;
4851   if (encoding)
4852     {
4853       free (read->rf->encoding);
4854       read->rf->encoding = encoding;
4855       encoding = NULL;
4856     }
4857
4858   /* Field width may be specified in multiple ways:
4859
4860      1. BY on FIELD.
4861      2. The format on FORMAT.
4862      3. The repetition factor on FORMAT.
4863
4864      (2) and (3) are mutually exclusive.
4865
4866      If more than one of these is present, they must agree.  If none of them is
4867      present, then free-field format is used.
4868    */
4869   if (repetitions > record_width)
4870     {
4871       msg (SE, _("%d repetitions cannot fit in record width %d."),
4872            repetitions, record_width);
4873       goto error;
4874     }
4875   int w = (repetitions ? record_width / repetitions
4876            : read->w ? read->w
4877            : by);
4878   if (by && w != by)
4879     {
4880       if (repetitions)
4881         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4882                    "which implies field width %d, "
4883                    "but BY specifies field width %d."),
4884              repetitions, record_width, w, by);
4885       else
4886         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4887              w, by);
4888       goto error;
4889     }
4890   read->w = w;
4891   return cmd;
4892
4893 error:
4894   fh_unref (fh);
4895   matrix_cmd_destroy (cmd);
4896   free (encoding);
4897   return NULL;
4898 }
4899
4900 static void
4901 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4902              struct substring data, size_t y, size_t x,
4903              int first_column, int last_column, char *error)
4904 {
4905   int line_number = dfm_get_line_number (reader);
4906   struct msg_location *location = xmalloc (sizeof *location);
4907   *location = (struct msg_location) {
4908     .file_name = xstrdup (dfm_get_file_name (reader)),
4909     .first_line = line_number,
4910     .last_line = line_number + 1,
4911     .first_column = first_column,
4912     .last_column = last_column,
4913   };
4914   struct msg *m = xmalloc (sizeof *m);
4915   *m = (struct msg) {
4916     .category = MSG_C_DATA,
4917     .severity = MSG_S_WARNING,
4918     .location = location,
4919     .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4920                          "for matrix row %zu, column %zu: %s"),
4921                        (int) data.length, data.string, fmt_name (format),
4922                        y + 1, x + 1, error),
4923   };
4924   msg_emit (m);
4925
4926   free (error);
4927 }
4928
4929 static void
4930 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4931                        gsl_matrix *m, struct substring p, size_t y, size_t x,
4932                        const char *line_start)
4933 {
4934   const char *input_encoding = dfm_reader_get_encoding (reader);
4935   char *error;
4936   double f;
4937   if (fmt_is_numeric (read->format))
4938     {
4939       union value v;
4940       error = data_in (p, input_encoding, read->format,
4941                        settings_get_fmt_settings (), &v, 0, NULL);
4942       if (!error && v.f == SYSMIS)
4943         error = xstrdup (_("Matrix data may not contain missing value."));
4944       f = v.f;
4945     }
4946     else
4947       {
4948         uint8_t s[sizeof (double)];
4949         union value v = { .s = s };
4950         error = data_in (p, input_encoding, read->format,
4951                          settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4952         memcpy (&f, s, sizeof f);
4953       }
4954
4955   if (error)
4956     {
4957       int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4958       int c2 = c1 + ss_utf8_count_columns (p) - 1;
4959       parse_error (reader, read->format, p, y, x, c1, c2, error);
4960     }
4961   else
4962     {
4963       gsl_matrix_set (m, y, x, f);
4964       if (read->symmetric && x != y)
4965         gsl_matrix_set (m, x, y, f);
4966     }
4967 }
4968
4969 static bool
4970 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4971                   struct substring *line, const char **startp)
4972 {
4973   if (dfm_eof (reader))
4974     {
4975       msg (SE, _("Unexpected end of file reading matrix data."));
4976       return false;
4977     }
4978   dfm_expand_tabs (reader);
4979   struct substring record = dfm_get_record (reader);
4980   /* XXX need to recode record into UTF-8 */
4981   *startp = record.string;
4982   *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4983   return true;
4984 }
4985
4986 static void
4987 matrix_read (struct read_command *read, struct dfm_reader *reader,
4988              gsl_matrix *m)
4989 {
4990   for (size_t y = 0; y < m->size1; y++)
4991     {
4992       size_t nx = read->symmetric ? y + 1 : m->size2;
4993
4994       struct substring line = ss_empty ();
4995       const char *line_start = line.string;
4996       for (size_t x = 0; x < nx; x++)
4997         {
4998           struct substring p;
4999           if (!read->w)
5000             {
5001               for (;;)
5002                 {
5003                   ss_ltrim (&line, ss_cstr (" ,"));
5004                   if (!ss_is_empty (line))
5005                     break;
5006                   if (!matrix_read_line (read, reader, &line, &line_start))
5007                     return;
5008                   dfm_forward_record (reader);
5009                 }
5010
5011               ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5012             }
5013           else
5014             {
5015               if (!matrix_read_line (read, reader, &line, &line_start))
5016                 return;
5017               size_t fields_per_line = (read->c2 - read->c1) / read->w;
5018               int f = x % fields_per_line;
5019               if (f == fields_per_line - 1)
5020                 dfm_forward_record (reader);
5021
5022               p = ss_substr (line, read->w * f, read->w);
5023             }
5024
5025           matrix_read_set_field (read, reader, m, p, y, x, line_start);
5026         }
5027
5028       if (read->w)
5029         dfm_forward_record (reader);
5030       else
5031         {
5032           ss_ltrim (&line, ss_cstr (" ,"));
5033           if (!ss_is_empty (line))
5034             {
5035               /* XXX */
5036               msg (SW, _("Trailing garbage on line \"%.*s\""),
5037                    (int) line.length, line.string);
5038             }
5039         }
5040     }
5041 }
5042
5043 static void
5044 matrix_cmd_execute_read (struct read_command *read)
5045 {
5046   struct index_vector iv0, iv1;
5047   if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5048     return;
5049
5050   size_t size[2] = { SIZE_MAX, SIZE_MAX };
5051   if (read->size)
5052     {
5053       gsl_matrix *m = matrix_expr_evaluate (read->size);
5054       if (!m)
5055         return;
5056
5057       if (!is_vector (m))
5058         {
5059           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5060                      "not a %zu×%zu matrix."), m->size1, m->size2);
5061           gsl_matrix_free (m);
5062           free (iv0.indexes);
5063           free (iv1.indexes);
5064           return;
5065         }
5066
5067       gsl_vector v = to_vector (m);
5068       double d[2];
5069       if (v.size == 1)
5070         {
5071           d[0] = gsl_vector_get (&v, 0);
5072           d[1] = 1;
5073         }
5074       else if (v.size == 2)
5075         {
5076           d[0] = gsl_vector_get (&v, 0);
5077           d[1] = gsl_vector_get (&v, 1);
5078         }
5079       else
5080         {
5081           msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5082                      "not a %zu×%zu matrix."),
5083                m->size1, m->size2),
5084           gsl_matrix_free (m);
5085           free (iv0.indexes);
5086           free (iv1.indexes);
5087           return;
5088         }
5089       gsl_matrix_free (m);
5090
5091       if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5092         {
5093           msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5094                      "are outside valid range."),
5095                d[0], d[1]);
5096           free (iv0.indexes);
5097           free (iv1.indexes);
5098           return;
5099         }
5100
5101       size[0] = d[0];
5102       size[1] = d[1];
5103     }
5104
5105   if (read->dst->n_indexes)
5106     {
5107       size_t submatrix_size[2];
5108       if (read->dst->n_indexes == 2)
5109         {
5110           submatrix_size[0] = iv0.n;
5111           submatrix_size[1] = iv1.n;
5112         }
5113       else if (read->dst->var->value->size1 == 1)
5114         {
5115           submatrix_size[0] = 1;
5116           submatrix_size[1] = iv0.n;
5117         }
5118       else
5119         {
5120           submatrix_size[0] = iv0.n;
5121           submatrix_size[1] = 1;
5122         }
5123
5124       if (read->size)
5125         {
5126           if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5127             {
5128               msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5129                          "differ from submatrix dimensions %zu×%zu."),
5130                    size[0], size[1],
5131                    submatrix_size[0], submatrix_size[1]);
5132               free (iv0.indexes);
5133               free (iv1.indexes);
5134               return;
5135             }
5136         }
5137       else
5138         {
5139           size[0] = submatrix_size[0];
5140           size[1] = submatrix_size[1];
5141         }
5142     }
5143
5144   struct dfm_reader *reader = read_file_open (read->rf);
5145   if (read->reread)
5146     dfm_reread_record (reader, 1);
5147
5148   if (read->symmetric && size[0] != size[1])
5149     {
5150       msg (SE, _("Cannot read non-square %zu×%zu matrix "
5151                  "using READ with MODE=SYMMETRIC."),
5152            size[0], size[1]);
5153       free (iv0.indexes);
5154       free (iv1.indexes);
5155       return;
5156     }
5157   gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5158   matrix_read (read, reader, tmp);
5159   matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5160 }
5161 \f
5162 static struct matrix_cmd *
5163 matrix_parse_write (struct matrix_state *s)
5164 {
5165   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5166   *cmd = (struct matrix_cmd) {
5167     .type = MCMD_WRITE,
5168   };
5169
5170   struct file_handle *fh = NULL;
5171   char *encoding = NULL;
5172   struct write_command *write = &cmd->write;
5173   write->expression = matrix_parse_exp (s);
5174   if (!write->expression)
5175     goto error;
5176
5177   int by = 0;
5178   int repetitions = 0;
5179   int record_width = 0;
5180   enum fmt_type format = FMT_F;
5181   bool has_format = false;
5182   while (lex_match (s->lexer, T_SLASH))
5183     {
5184       if (lex_match_id (s->lexer, "OUTFILE"))
5185         {
5186           lex_match (s->lexer, T_EQUALS);
5187
5188           fh_unref (fh);
5189           fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5190           if (!fh)
5191             goto error;
5192         }
5193       else if (lex_match_id (s->lexer, "ENCODING"))
5194         {
5195           lex_match (s->lexer, T_EQUALS);
5196           if (!lex_force_string (s->lexer))
5197             goto error;
5198
5199           free (encoding);
5200           encoding = ss_xstrdup (lex_tokss (s->lexer));
5201
5202           lex_get (s->lexer);
5203         }
5204       else if (lex_match_id (s->lexer, "FIELD"))
5205         {
5206           lex_match (s->lexer, T_EQUALS);
5207
5208           if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5209             goto error;
5210           write->c1 = lex_integer (s->lexer);
5211           lex_get (s->lexer);
5212           if (!lex_force_match (s->lexer, T_TO)
5213               || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5214             goto error;
5215           write->c2 = lex_integer (s->lexer) + 1;
5216           lex_get (s->lexer);
5217
5218           record_width = write->c2 - write->c1;
5219           if (lex_match (s->lexer, T_BY))
5220             {
5221               if (!lex_force_int_range (s->lexer, "BY", 1,
5222                                         write->c2 - write->c1))
5223                 goto error;
5224               by = lex_integer (s->lexer);
5225               lex_get (s->lexer);
5226
5227               if (record_width % by)
5228                 {
5229                   msg (SE, _("BY %d does not evenly divide record width %d."),
5230                        by, record_width);
5231                   goto error;
5232                 }
5233             }
5234           else
5235             by = 0;
5236         }
5237       else if (lex_match_id (s->lexer, "MODE"))
5238         {
5239           lex_match (s->lexer, T_EQUALS);
5240           if (lex_match_id (s->lexer, "RECTANGULAR"))
5241             write->triangular = false;
5242           else if (lex_match_id (s->lexer, "TRIANGULAR"))
5243             write->triangular = true;
5244           else
5245             {
5246               lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5247               goto error;
5248             }
5249         }
5250       else if (lex_match_id (s->lexer, "HOLD"))
5251         write->hold = true;
5252       else if (lex_match_id (s->lexer, "FORMAT"))
5253         {
5254           if (has_format || write->format)
5255             {
5256               lex_sbc_only_once ("FORMAT");
5257               goto error;
5258             }
5259
5260           lex_match (s->lexer, T_EQUALS);
5261
5262           if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5263             goto error;
5264
5265           const char *p = lex_tokcstr (s->lexer);
5266           if (c_isdigit (p[0]))
5267             {
5268               repetitions = atoi (p);
5269               p += strspn (p, "0123456789");
5270               if (!fmt_from_name (p, &format))
5271                 {
5272                   lex_error (s->lexer, _("Unknown format %s."), p);
5273                   goto error;
5274                 }
5275               has_format = true;
5276               lex_get (s->lexer);
5277             }
5278           else if (fmt_from_name (p, &format))
5279             {
5280               has_format = true;
5281               lex_get (s->lexer);
5282             }
5283           else
5284             {
5285               struct fmt_spec spec;
5286               if (!parse_format_specifier (s->lexer, &spec))
5287                 goto error;
5288               write->format = xmemdup (&spec, sizeof spec);
5289             }
5290         }
5291       else
5292         {
5293           lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5294                                "HOLD", "FORMAT");
5295           goto error;
5296         }
5297     }
5298
5299   if (!write->c1)
5300     {
5301       lex_sbc_missing ("FIELD");
5302       goto error;
5303     }
5304
5305   if (!fh)
5306     {
5307       if (s->prev_write_file)
5308         fh = fh_ref (s->prev_write_file);
5309       else
5310         {
5311           lex_sbc_missing ("OUTFILE");
5312           goto error;
5313         }
5314     }
5315   fh_unref (s->prev_write_file);
5316   s->prev_write_file = fh_ref (fh);
5317
5318   write->wf = write_file_create (s, fh);
5319   fh = NULL;
5320   if (encoding)
5321     {
5322       free (write->wf->encoding);
5323       write->wf->encoding = encoding;
5324       encoding = NULL;
5325     }
5326
5327   /* Field width may be specified in multiple ways:
5328
5329      1. BY on FIELD.
5330      2. The format on FORMAT.
5331      3. The repetition factor on FORMAT.
5332
5333      (2) and (3) are mutually exclusive.
5334
5335      If more than one of these is present, they must agree.  If none of them is
5336      present, then free-field format is used.
5337    */
5338   if (repetitions > record_width)
5339     {
5340       msg (SE, _("%d repetitions cannot fit in record width %d."),
5341            repetitions, record_width);
5342       goto error;
5343     }
5344   int w = (repetitions ? record_width / repetitions
5345            : write->format ? write->format->w
5346            : by);
5347   if (by && w != by)
5348     {
5349       if (repetitions)
5350         msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5351                    "which implies field width %d, "
5352                    "but BY specifies field width %d."),
5353              repetitions, record_width, w, by);
5354       else
5355         msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5356              w, by);
5357       goto error;
5358     }
5359   if (w && !write->format)
5360     {
5361       write->format = xmalloc (sizeof *write->format);
5362       *write->format = (struct fmt_spec) { .type = format, .w = w };
5363
5364       if (!fmt_check_output (write->format))
5365         goto error;
5366     };
5367
5368   if (write->format && fmt_var_width (write->format) > sizeof (double))
5369     {
5370       char s[FMT_STRING_LEN_MAX + 1];
5371       fmt_to_string (write->format, s);
5372       msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5373            s, sizeof (double));
5374       goto error;
5375     }
5376
5377   return cmd;
5378
5379 error:
5380   fh_unref (fh);
5381   matrix_cmd_destroy (cmd);
5382   return NULL;
5383 }
5384
5385 static void
5386 matrix_cmd_execute_write (struct write_command *write)
5387 {
5388   gsl_matrix *m = matrix_expr_evaluate (write->expression);
5389   if (!m)
5390     return;
5391
5392   if (write->triangular && m->size1 != m->size2)
5393     {
5394       msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5395                  "the matrix to be written has dimensions %zu×%zu."),
5396            m->size1, m->size2);
5397       gsl_matrix_free (m);
5398       return;
5399     }
5400
5401   struct dfm_writer *writer = write_file_open (write->wf);
5402   if (!writer || !m->size1)
5403     {
5404       gsl_matrix_free (m);
5405       return;
5406     }
5407
5408   const struct fmt_settings *settings = settings_get_fmt_settings ();
5409   struct u8_line *line = write->wf->held;
5410   for (size_t y = 0; y < m->size1; y++)
5411     {
5412       if (!line)
5413         {
5414           line = xmalloc (sizeof *line);
5415           u8_line_init (line);
5416         }
5417       size_t nx = write->triangular ? y + 1 : m->size2;
5418       int x0 = write->c1;
5419       for (size_t x = 0; x < nx; x++)
5420         {
5421           char *s;
5422           double f = gsl_matrix_get (m, y, x);
5423           if (write->format)
5424             {
5425               union value v;
5426               if (fmt_is_numeric (write->format->type))
5427                 v.f = f;
5428               else
5429                 v.s = (uint8_t *) &f;
5430               s = data_out (&v, NULL, write->format, settings);
5431             }
5432           else
5433             {
5434               s = xmalloc (DBL_BUFSIZE_BOUND);
5435               if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5436                   >= DBL_BUFSIZE_BOUND)
5437                 abort ();
5438             }
5439           size_t len = strlen (s);
5440           int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5441           if (width + x0 > write->c2)
5442             {
5443               dfm_put_record_utf8 (writer, line->s.ss.string,
5444                                    line->s.ss.length);
5445               u8_line_clear (line);
5446               x0 = write->c1;
5447             }
5448           u8_line_put (line, x0, x0 + width, s, len);
5449           free (s);
5450
5451           x0 += write->format ? write->format->w : width + 1;
5452         }
5453
5454       if (y + 1 >= m->size1 && write->hold)
5455         break;
5456       dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5457       u8_line_clear (line);
5458     }
5459   if (!write->hold)
5460     {
5461       u8_line_destroy (line);
5462       free (line);
5463       line = NULL;
5464     }
5465   write->wf->held = line;
5466
5467   gsl_matrix_free (m);
5468 }
5469 \f
5470 static struct matrix_cmd *
5471 matrix_parse_get (struct matrix_state *s)
5472 {
5473   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5474   *cmd = (struct matrix_cmd) {
5475     .type = MCMD_GET,
5476     .get = {
5477       .dataset = s->dataset,
5478       .user = { .treatment = MGET_ERROR },
5479       .system = { .treatment = MGET_ERROR },
5480     }
5481   };
5482
5483   struct get_command *get = &cmd->get;
5484   get->dst = matrix_lvalue_parse (s);
5485   if (!get->dst)
5486     goto error;
5487
5488   while (lex_match (s->lexer, T_SLASH))
5489     {
5490       if (lex_match_id (s->lexer, "FILE"))
5491         {
5492           lex_match (s->lexer, T_EQUALS);
5493
5494           fh_unref (get->file);
5495           if (lex_match (s->lexer, T_ASTERISK))
5496             get->file = NULL;
5497           else
5498             {
5499               get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5500               if (!get->file)
5501                 goto error;
5502             }
5503         }
5504       else if (lex_match_id (s->lexer, "ENCODING"))
5505         {
5506           lex_match (s->lexer, T_EQUALS);
5507           if (!lex_force_string (s->lexer))
5508             goto error;
5509
5510           free (get->encoding);
5511           get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5512
5513           lex_get (s->lexer);
5514         }
5515       else if (lex_match_id (s->lexer, "VARIABLES"))
5516         {
5517           lex_match (s->lexer, T_EQUALS);
5518
5519           if (get->n_vars)
5520             {
5521               lex_sbc_only_once ("VARIABLES");
5522               goto error;
5523             }
5524
5525           if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5526             goto error;
5527         }
5528       else if (lex_match_id (s->lexer, "NAMES"))
5529         {
5530           lex_match (s->lexer, T_EQUALS);
5531           if (!lex_force_id (s->lexer))
5532             goto error;
5533
5534           struct substring name = lex_tokss (s->lexer);
5535           get->names = matrix_var_lookup (s, name);
5536           if (!get->names)
5537             get->names = matrix_var_create (s, name);
5538           lex_get (s->lexer);
5539         }
5540       else if (lex_match_id (s->lexer, "MISSING"))
5541         {
5542           lex_match (s->lexer, T_EQUALS);
5543           if (lex_match_id (s->lexer, "ACCEPT"))
5544             get->user.treatment = MGET_ACCEPT;
5545           else if (lex_match_id (s->lexer, "OMIT"))
5546             get->user.treatment = MGET_OMIT;
5547           else if (lex_is_number (s->lexer))
5548             {
5549               get->user.treatment = MGET_RECODE;
5550               get->user.substitute = lex_number (s->lexer);
5551               lex_get (s->lexer);
5552             }
5553           else
5554             {
5555               lex_error (s->lexer, NULL);
5556               goto error;
5557             }
5558         }
5559       else if (lex_match_id (s->lexer, "SYSMIS"))
5560         {
5561           lex_match (s->lexer, T_EQUALS);
5562           if (lex_match_id (s->lexer, "OMIT"))
5563             get->system.treatment = MGET_OMIT;
5564           else if (lex_is_number (s->lexer))
5565             {
5566               get->system.treatment = MGET_RECODE;
5567               get->system.substitute = lex_number (s->lexer);
5568               lex_get (s->lexer);
5569             }
5570           else
5571             {
5572               lex_error (s->lexer, NULL);
5573               goto error;
5574             }
5575         }
5576       else
5577         {
5578           lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5579                                "MISSING", "SYSMIS");
5580           goto error;
5581         }
5582     }
5583
5584   if (get->user.treatment != MGET_ACCEPT)
5585     get->system.treatment = MGET_ERROR;
5586
5587   return cmd;
5588
5589 error:
5590   matrix_cmd_destroy (cmd);
5591   return NULL;
5592 }
5593
5594 static void
5595 matrix_cmd_execute_get__ (struct get_command *get,
5596                           const struct dictionary *dict,
5597                           struct casereader *reader)
5598 {
5599   struct variable **vars;
5600   size_t n_vars = 0;
5601
5602   if (get->n_vars)
5603     {
5604       if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5605                                 &vars, &n_vars, PV_NUMERIC))
5606         return;
5607     }
5608   else
5609     {
5610       n_vars = dict_get_var_cnt (dict);
5611       vars = xnmalloc (n_vars, sizeof *vars);
5612       for (size_t i = 0; i < n_vars; i++)
5613         {
5614           struct variable *var = dict_get_var (dict, i);
5615           if (!var_is_numeric (var))
5616             {
5617               msg (SE, _("GET: Variable %s is not numeric."),
5618                    var_get_name (var));
5619               free (vars);
5620               return;
5621             }
5622           vars[i] = var;
5623         }
5624     }
5625
5626   if (get->names)
5627     {
5628       gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5629       for (size_t i = 0; i < n_vars; i++)
5630         {
5631           char s[sizeof (double)];
5632           double f;
5633           buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5634           memcpy (&f, s, sizeof f);
5635           gsl_matrix_set (names, i, 0, f);
5636         }
5637
5638       gsl_matrix_free (get->names->value);
5639       get->names->value = names;
5640     }
5641
5642   size_t n_rows = 0;
5643   gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5644   long long int casenum = 1;
5645   bool error = false;
5646   for (struct ccase *c = casereader_read (reader); c;
5647        c = casereader_read (reader), casenum++)
5648     {
5649       if (n_rows >= m->size1)
5650         {
5651           gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5652           for (size_t y = 0; y < n_rows; y++)
5653             for (size_t x = 0; x < n_vars; x++)
5654               gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5655           gsl_matrix_free (m);
5656           m = p;
5657         }
5658
5659       bool keep = true;
5660       for (size_t x = 0; x < n_vars; x++)
5661         {
5662           const struct variable *var = vars[x];
5663           double d = case_num (c, var);
5664           if (d == SYSMIS)
5665             {
5666               if (get->system.treatment == MGET_RECODE)
5667                 d = get->system.substitute;
5668               else if (get->system.treatment == MGET_OMIT)
5669                 keep = false;
5670               else
5671                 {
5672                   msg (SE, _("GET: Variable %s in case %lld "
5673                              "is system-missing."),
5674                        var_get_name (var), casenum);
5675                   error = true;
5676                 }
5677             }
5678           else if (var_is_num_missing (var, d, MV_USER))
5679             {
5680               if (get->user.treatment == MGET_RECODE)
5681                 d = get->user.substitute;
5682               else if (get->user.treatment == MGET_OMIT)
5683                 keep = false;
5684               else if (get->user.treatment != MGET_ACCEPT)
5685                 {
5686                   msg (SE, _("GET: Variable %s in case %lld has user-missing "
5687                              "value %g."),
5688                        var_get_name (var), casenum, d);
5689                   error = true;
5690                 }
5691             }
5692           gsl_matrix_set (m, n_rows, x, d);
5693         }
5694       case_unref (c);
5695       if (error)
5696         break;
5697       if (keep)
5698         n_rows++;
5699     }
5700   if (!error)
5701     {
5702       m->size1 = n_rows;
5703       matrix_lvalue_evaluate_and_assign (get->dst, m);
5704     }
5705   else
5706     gsl_matrix_free (m);
5707   free (vars);
5708 }
5709
5710 static void
5711 matrix_cmd_execute_get (struct get_command *get)
5712 {
5713   struct dictionary *dict;
5714   struct casereader *reader;
5715   if (get->file)
5716     {
5717        reader = any_reader_open_and_decode (get->file, get->encoding,
5718                                             &dict, NULL);
5719        if (!reader)
5720          return;
5721     }
5722   else
5723     {
5724       if (dict_get_var_cnt (dataset_dict (get->dataset)) == 0)
5725         {
5726           msg (ME, _("GET cannot read empty active file."));
5727           return;
5728         }
5729       reader = proc_open (get->dataset);
5730       dict = dict_ref (dataset_dict (get->dataset));
5731     }
5732
5733   matrix_cmd_execute_get__ (get, dict, reader);
5734
5735   dict_unref (dict);
5736   casereader_destroy (reader);
5737   if (!get->file)
5738     proc_commit (get->dataset);
5739 }
5740 \f
5741 static const char *
5742 match_rowtype (struct lexer *lexer)
5743 {
5744   static const char *rowtypes[] = {
5745     "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5746   };
5747   size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5748
5749   for (size_t i = 0; i < n_rowtypes; i++)
5750     if (lex_match_id (lexer, rowtypes[i]))
5751       return rowtypes[i];
5752
5753   lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5754   return NULL;
5755 }
5756
5757 static bool
5758 parse_var_names (struct lexer *lexer, struct string_array *sa)
5759 {
5760   lex_match (lexer, T_EQUALS);
5761
5762   string_array_clear (sa);
5763
5764   struct dictionary *dict = dict_create (get_default_encoding ());
5765   dict_create_var_assert (dict, "ROWTYPE_", 8);
5766   dict_create_var_assert (dict, "VARNAME_", 8);
5767   char **names;
5768   size_t n_names;
5769   bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5770                                   PV_NO_DUPLICATE | PV_NO_SCRATCH);
5771   dict_unref (dict);
5772
5773   if (ok)
5774     {
5775       string_array_clear (sa);
5776       sa->strings = names;
5777       sa->n = sa->allocated = n_names;
5778     }
5779   return ok;
5780 }
5781
5782 static void
5783 msave_common_uninit (struct msave_common *common)
5784 {
5785   if (common)
5786     {
5787       fh_unref (common->outfile);
5788       string_array_destroy (&common->variables);
5789       string_array_destroy (&common->fnames);
5790       string_array_destroy (&common->snames);
5791     }
5792 }
5793
5794 static bool
5795 compare_variables (const char *keyword,
5796                    const struct string_array *new,
5797                    const struct string_array *old)
5798 {
5799   if (new->n)
5800     {
5801       if (!old->n)
5802         {
5803           msg (SE, _("%s may only be specified on MSAVE if it was specified "
5804                      "on the first MSAVE within MATRIX."), keyword);
5805           return false;
5806         }
5807       else if (!string_array_equal_case (old, new))
5808         {
5809           msg (SE, _("%s must specify the same variables each time within "
5810                      "a given MATRIX."), keyword);
5811           return false;
5812         }
5813     }
5814   return true;
5815 }
5816
5817 static struct matrix_cmd *
5818 matrix_parse_msave (struct matrix_state *s)
5819 {
5820   struct msave_common common = { .outfile = NULL };
5821   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5822   *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5823
5824   struct msave_command *msave = &cmd->msave;
5825   if (lex_token (s->lexer) == T_ID)
5826     msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5827   msave->expr = matrix_parse_exp (s);
5828   if (!msave->expr)
5829     return NULL;
5830
5831   while (lex_match (s->lexer, T_SLASH))
5832     {
5833       if (lex_match_id (s->lexer, "TYPE"))
5834         {
5835           lex_match (s->lexer, T_EQUALS);
5836
5837           msave->rowtype = match_rowtype (s->lexer);
5838           if (!msave->rowtype)
5839             goto error;
5840         }
5841       else if (lex_match_id (s->lexer, "OUTFILE"))
5842         {
5843           lex_match (s->lexer, T_EQUALS);
5844
5845           fh_unref (common.outfile);
5846           common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5847           if (!common.outfile)
5848             goto error;
5849         }
5850       else if (lex_match_id (s->lexer, "VARIABLES"))
5851         {
5852           if (!parse_var_names (s->lexer, &common.variables))
5853             goto error;
5854         }
5855       else if (lex_match_id (s->lexer, "FNAMES"))
5856         {
5857           if (!parse_var_names (s->lexer, &common.fnames))
5858             goto error;
5859         }
5860       else if (lex_match_id (s->lexer, "SNAMES"))
5861         {
5862           if (!parse_var_names (s->lexer, &common.snames))
5863             goto error;
5864         }
5865       else if (lex_match_id (s->lexer, "SPLIT"))
5866         {
5867           lex_match (s->lexer, T_EQUALS);
5868
5869           matrix_expr_destroy (msave->splits);
5870           msave->splits = matrix_parse_exp (s);
5871           if (!msave->splits)
5872             goto error;
5873         }
5874       else if (lex_match_id (s->lexer, "FACTOR"))
5875         {
5876           lex_match (s->lexer, T_EQUALS);
5877
5878           matrix_expr_destroy (msave->factors);
5879           msave->factors = matrix_parse_exp (s);
5880           if (!msave->factors)
5881             goto error;
5882         }
5883       else
5884         {
5885           lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5886                                "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5887           goto error;
5888         }
5889     }
5890   if (!msave->rowtype)
5891     {
5892       lex_sbc_missing ("TYPE");
5893       goto error;
5894     }
5895   common.has_splits = msave->splits || common.snames.n;
5896   common.has_factors = msave->factors || common.fnames.n;
5897
5898   struct msave_common *c = s->common ? s->common : &common;
5899   if (c->fnames.n && !msave->factors)
5900     {
5901       msg (SE, _("FNAMES requires FACTOR."));
5902       goto error;
5903     }
5904   if (c->snames.n && !msave->splits)
5905     {
5906       msg (SE, _("SNAMES requires SPLIT."));
5907       goto error;
5908     }
5909   if (c->has_factors && !common.has_factors)
5910     {
5911       msg (SE, _("%s is required because it was present on the first "
5912                  "MSAVE in this MATRIX command."), "FACTOR");
5913       goto error;
5914     }
5915   if (c->has_splits && !common.has_splits)
5916     {
5917       msg (SE, _("%s is required because it was present on the first "
5918                  "MSAVE in this MATRIX command."), "SPLIT");
5919       goto error;
5920     }
5921
5922   if (!s->common)
5923     {
5924       if (!common.outfile)
5925         {
5926           lex_sbc_missing ("OUTFILE");
5927           goto error;
5928         }
5929       s->common = xmemdup (&common, sizeof common);
5930     }
5931   else
5932     {
5933       if (common.outfile)
5934         {
5935           bool same = common.outfile == s->common->outfile;
5936           fh_unref (common.outfile);
5937           if (!same)
5938             {
5939               msg (SE, _("OUTFILE must name the same file on each MSAVE "
5940                          "within a single MATRIX command."));
5941               goto error;
5942             }
5943         }
5944       if (!compare_variables ("VARIABLES",
5945                               &common.variables, &s->common->variables)
5946           || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5947           || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5948         goto error;
5949       msave_common_uninit (&common);
5950     }
5951   msave->common = s->common;
5952   if (!msave->varname_)
5953     msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5954   return cmd;
5955
5956 error:
5957   msave_common_uninit (&common);
5958   matrix_cmd_destroy (cmd);
5959   return NULL;
5960 }
5961
5962 static gsl_vector *
5963 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5964 {
5965   gsl_matrix *m = matrix_expr_evaluate (e);
5966   if (!m)
5967     return NULL;
5968
5969   if (!is_vector (m))
5970     {
5971       msg (SE, _("%s expression must evaluate to vector, "
5972                  "not a %zu×%zu matrix."),
5973            name, m->size1, m->size2);
5974       gsl_matrix_free (m);
5975       return NULL;
5976     }
5977
5978   return matrix_to_vector (m);
5979 }
5980
5981 static const char *
5982 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5983 {
5984   for (size_t i = 0; i < vars->n; i++)
5985     if (!dict_create_var (d, vars->strings[i], 0))
5986       return vars->strings[i];
5987   return NULL;
5988 }
5989
5990 static struct dictionary *
5991 msave_create_dict (const struct msave_common *common)
5992 {
5993   struct dictionary *dict = dict_create (get_default_encoding ());
5994
5995   const char *dup_split = msave_add_vars (dict, &common->snames);
5996   if (dup_split)
5997     {
5998       msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5999       goto error;
6000     }
6001
6002   dict_create_var_assert (dict, "ROWTYPE_", 8);
6003
6004   const char *dup_factor = msave_add_vars (dict, &common->fnames);
6005   if (dup_factor)
6006     {
6007       msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
6008       goto error;
6009     }
6010
6011   dict_create_var_assert (dict, "VARNAME_", 8);
6012
6013   const char *dup_var = msave_add_vars (dict, &common->variables);
6014   if (dup_var)
6015     {
6016       msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6017       goto error;
6018     }
6019
6020   return dict;
6021
6022 error:
6023   dict_unref (dict);
6024   return NULL;
6025 }
6026
6027 static void
6028 matrix_cmd_execute_msave (struct msave_command *msave)
6029 {
6030   struct msave_common *common = msave->common;
6031   gsl_matrix *m = NULL;
6032   gsl_vector *factors = NULL;
6033   gsl_vector *splits = NULL;
6034
6035   m = matrix_expr_evaluate (msave->expr);
6036   if (!m)
6037     goto error;
6038
6039   if (!common->variables.n)
6040     for (size_t i = 0; i < m->size2; i++)
6041       string_array_append_nocopy (&common->variables,
6042                                   xasprintf ("COL%zu", i + 1));
6043
6044   if (m->size2 != common->variables.n)
6045     {
6046       msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
6047            m->size2, common->variables.n);
6048       goto error;
6049     }
6050
6051   if (msave->factors)
6052     {
6053       factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6054       if (!factors)
6055         goto error;
6056
6057       if (!common->fnames.n)
6058         for (size_t i = 0; i < factors->size; i++)
6059           string_array_append_nocopy (&common->fnames,
6060                                       xasprintf ("FAC%zu", i + 1));
6061
6062       if (factors->size != common->fnames.n)
6063         {
6064           msg (SE, _("There are %zu factor variables, "
6065                      "but %zu split values were supplied."),
6066                common->fnames.n, factors->size);
6067           goto error;
6068         }
6069     }
6070
6071   if (msave->splits)
6072     {
6073       splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6074       if (!splits)
6075         goto error;
6076
6077       if (!common->fnames.n)
6078         for (size_t i = 0; i < splits->size; i++)
6079           string_array_append_nocopy (&common->fnames,
6080                                       xasprintf ("SPL%zu", i + 1));
6081
6082       if (splits->size != common->fnames.n)
6083         {
6084           msg (SE, _("There are %zu split variables, "
6085                      "but %zu split values were supplied."),
6086                common->fnames.n, splits->size);
6087           goto error;
6088         }
6089     }
6090
6091   if (!common->writer)
6092     {
6093       struct dictionary *dict = msave_create_dict (common);
6094       if (!dict)
6095         goto error;
6096
6097       common->writer = any_writer_open (common->outfile, dict);
6098       if (!common->writer)
6099         {
6100           dict_unref (dict);
6101           goto error;
6102         }
6103
6104       common->dict = dict;
6105     }
6106
6107   for (size_t y = 0; y < m->size1; y++)
6108     {
6109       struct ccase *c = case_create (dict_get_proto (common->dict));
6110       size_t idx = 0;
6111
6112       /* Split variables */
6113       if (splits)
6114         for (size_t i = 0; i < splits->size; i++)
6115           case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6116
6117       /* ROWTYPE_. */
6118       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6119                          msave->rowtype, ' ');
6120
6121       /* Factors. */
6122       if (factors)
6123         for (size_t i = 0; i < factors->size; i++)
6124           *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6125
6126       /* VARNAME_. */
6127       buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6128                          msave->varname_, ' ');
6129
6130       /* Continuous variables. */
6131       for (size_t x = 0; x < m->size2; x++)
6132         case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6133       casewriter_write (common->writer, c);
6134     }
6135
6136 error:
6137   gsl_matrix_free (m);
6138   gsl_vector_free (factors);
6139   gsl_vector_free (splits);
6140 }
6141 \f
6142 static struct matrix_cmd *
6143 matrix_parse_mget (struct matrix_state *s)
6144 {
6145   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6146   *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
6147
6148   struct mget_command *mget = &cmd->mget;
6149
6150   for (;;)
6151     {
6152       if (lex_match_id (s->lexer, "FILE"))
6153         {
6154           lex_match (s->lexer, T_EQUALS);
6155
6156           fh_unref (mget->file);
6157           mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6158           if (!mget->file)
6159             goto error;
6160         }
6161       else if (lex_match_id (s->lexer, "TYPE"))
6162         {
6163           lex_match (s->lexer, T_EQUALS);
6164           while (lex_token (s->lexer) != T_SLASH
6165                  && lex_token (s->lexer) != T_ENDCMD)
6166             {
6167               const char *rowtype = match_rowtype (s->lexer);
6168               if (!rowtype)
6169                 goto error;
6170
6171               stringi_set_insert (&mget->rowtypes, rowtype);
6172             }
6173         }
6174       else
6175         {
6176           lex_error_expecting (s->lexer, "FILE", "TYPE");
6177           goto error;
6178         }
6179       if (lex_token (s->lexer) == T_ENDCMD)
6180         break;
6181
6182       if (!lex_force_match (s->lexer, T_SLASH))
6183         goto error;
6184     }
6185   return cmd;
6186
6187 error:
6188   matrix_cmd_destroy (cmd);
6189   return NULL;
6190 }
6191
6192 static const struct variable *
6193 get_a8_var (const struct dictionary *d, const char *name)
6194 {
6195   const struct variable *v = dict_lookup_var (d, name);
6196   if (!v)
6197     {
6198       msg (SE, _("Matrix data file lacks %s variable."), name);
6199       return NULL;
6200     }
6201   if (var_get_width (v) != 8)
6202     {
6203       msg (SE, _("%s variable in matrix data file must be "
6204                  "8-byte string, but it has width %d."),
6205            name, var_get_width (v));
6206       return NULL;
6207     }
6208   return v;
6209 }
6210
6211 static bool
6212 is_rowtype (const union value *v, const char *rowtype)
6213 {
6214   struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6215   ss_rtrim (&vs, ss_cstr (" "));
6216   return ss_equals_case (vs, ss_cstr (rowtype));
6217 }
6218
6219 static void
6220 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6221                         const struct dictionary *d,
6222                         const struct variable *rowtype_var,
6223                         struct matrix_state *s, size_t si, size_t fi,
6224                         size_t cs, size_t cn)
6225 {
6226   if (!n_rows)
6227     return;
6228
6229   const union value *rowtype_ = case_data (rows[0], rowtype_var);
6230   const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6231                              : is_rowtype (rowtype_, "CORR") ? "CR"
6232                              : is_rowtype (rowtype_, "MEAN") ? "MN"
6233                              : is_rowtype (rowtype_, "STDDEV") ? "SD"
6234                              : is_rowtype (rowtype_, "N") ? "NC"
6235                              : "CN");
6236
6237   struct string name = DS_EMPTY_INITIALIZER;
6238   ds_put_cstr (&name, name_prefix);
6239   if (fi > 0)
6240     ds_put_format (&name, "F%zu", fi);
6241   if (si > 0)
6242     ds_put_format (&name, "S%zu", si);
6243
6244   struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6245   if (!mv)
6246     mv = matrix_var_create (s, ds_ss (&name));
6247   else if (mv->value)
6248     {
6249       msg (SW, _("Matrix data file contains variable with existing name %s."),
6250            ds_cstr (&name));
6251       goto exit;
6252     }
6253
6254   gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6255   size_t n_missing = 0;
6256   for (size_t y = 0; y < n_rows; y++)
6257     {
6258       for (size_t x = 0; x < cn; x++)
6259         {
6260           struct variable *var = dict_get_var (d, cs + x);
6261           double value = case_num (rows[y], var);
6262           if (var_is_num_missing (var, value, MV_ANY))
6263             {
6264               n_missing++;
6265               value = 0.0;
6266             }
6267           gsl_matrix_set (m, y, x, value);
6268         }
6269     }
6270
6271   if (n_missing)
6272     msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6273                        "value, which was treated as zero.",
6274                        "Matrix data file variable %s contains %zu missing "
6275                        "values, which were treated as zero.", n_missing),
6276          ds_cstr (&name), n_missing);
6277   mv->value = m;
6278
6279 exit:
6280   ds_destroy (&name);
6281   for (size_t y = 0; y < n_rows; y++)
6282     case_unref (rows[y]);
6283 }
6284
6285 static bool
6286 var_changed (const struct ccase *ca, const struct ccase *cb,
6287              const struct variable *var)
6288 {
6289   return !value_equal (case_data (ca, var), case_data (cb, var),
6290                        var_get_width (var));
6291 }
6292
6293 static bool
6294 vars_changed (const struct ccase *ca, const struct ccase *cb,
6295               const struct dictionary *d,
6296               size_t first_var, size_t n_vars)
6297 {
6298   for (size_t i = 0; i < n_vars; i++)
6299     {
6300       const struct variable *v = dict_get_var (d, first_var + i);
6301       if (var_changed (ca, cb, v))
6302         return true;
6303     }
6304   return false;
6305 }
6306
6307 static void
6308 matrix_cmd_execute_mget (struct mget_command *mget)
6309 {
6310   struct dictionary *d;
6311   struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6312                                                      &d, NULL);
6313   if (!r)
6314     return;
6315
6316   const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6317   const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6318   if (!rowtype_ || !varname_)
6319     goto exit;
6320
6321   if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6322     {
6323       msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6324       goto exit;
6325     }
6326   if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6327     {
6328       msg (SE, _("Matrix data file contains no continuous variables."));
6329       goto exit;
6330     }
6331
6332   for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6333     {
6334       const struct variable *v = dict_get_var (d, i);
6335       if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6336         {
6337           msg (SE,
6338                _("Matrix data file contains unexpected string variable %s."),
6339                var_get_name (v));
6340           goto exit;
6341         }
6342     }
6343
6344   /* SPLIT variables. */
6345   size_t ss = 0;
6346   size_t sn = var_get_dict_index (rowtype_);
6347   struct ccase *sc = NULL;
6348   size_t si = 0;
6349
6350   /* FACTOR variables. */
6351   size_t fs = var_get_dict_index (rowtype_) + 1;
6352   size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6353   struct ccase *fc = NULL;
6354   size_t fi = 0;
6355
6356   /* Continuous variables. */
6357   size_t cs = var_get_dict_index (varname_) + 1;
6358   size_t cn = dict_get_var_cnt (d) - cs;
6359   struct ccase *cc = NULL;
6360
6361   /* Matrix. */
6362   struct ccase **rows = NULL;
6363   size_t allocated_rows = 0;
6364   size_t n_rows = 0;
6365
6366   struct ccase *c;
6367   while ((c = casereader_read (r)) != NULL)
6368     {
6369       bool sd = vars_changed (sc, c, d, ss, sn);
6370       bool fd = sd || vars_changed (fc, c, d, fs, fn);
6371       bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6372       if (sd)
6373         {
6374           si++;
6375           case_unref (sc);
6376           sc = case_ref (c);
6377         }
6378       if (fd)
6379         {
6380           fi++;
6381           case_unref (fc);
6382           fc = case_ref (c);
6383         }
6384       if (md)
6385         {
6386           matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6387                                   mget->state, si, fi, cs, cn);
6388           n_rows = 0;
6389           case_unref (cc);
6390           cc = case_ref (c);
6391         }
6392
6393       if (n_rows >= allocated_rows)
6394         rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6395       rows[n_rows++] = c;
6396     }
6397   matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6398                           mget->state, si, fi, cs, cn);
6399   free (rows);
6400
6401 exit:
6402   casereader_destroy (r);
6403 }
6404 \f
6405 static bool
6406 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6407 {
6408   if (!lex_force_id (s->lexer))
6409     return false;
6410
6411   *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6412   if (!*varp)
6413     *varp = matrix_var_create (s, lex_tokss (s->lexer));
6414   lex_get (s->lexer);
6415   return true;
6416 }
6417
6418 static struct matrix_cmd *
6419 matrix_parse_eigen (struct matrix_state *s)
6420 {
6421   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6422   *cmd = (struct matrix_cmd) {
6423     .type = MCMD_EIGEN,
6424     .eigen = { .expr = NULL }
6425   };
6426
6427   struct eigen_command *eigen = &cmd->eigen;
6428   if (!lex_force_match (s->lexer, T_LPAREN))
6429     goto error;
6430   eigen->expr = matrix_parse_expr (s);
6431   if (!eigen->expr
6432       || !lex_force_match (s->lexer, T_COMMA)
6433       || !matrix_parse_dst_var (s, &eigen->evec)
6434       || !lex_force_match (s->lexer, T_COMMA)
6435       || !matrix_parse_dst_var (s, &eigen->eval)
6436       || !lex_force_match (s->lexer, T_RPAREN))
6437     goto error;
6438
6439   return cmd;
6440
6441 error:
6442   matrix_cmd_destroy (cmd);
6443   return NULL;
6444 }
6445
6446 static void
6447 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6448 {
6449   gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6450   if (!A)
6451     return;
6452   if (!is_symmetric (A))
6453     {
6454       msg (SE, _("Argument of EIGEN must be symmetric."));
6455       gsl_matrix_free (A);
6456       return;
6457     }
6458
6459   gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6460   gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6461   gsl_vector v_eval = to_vector (eval);
6462   gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6463   gsl_eigen_symmv (A, &v_eval, evec, w);
6464   gsl_eigen_symmv_free (w);
6465
6466   gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6467
6468   gsl_matrix_free (A);
6469
6470   gsl_matrix_free (eigen->eval->value);
6471   eigen->eval->value = eval;
6472
6473   gsl_matrix_free (eigen->evec->value);
6474   eigen->evec->value = evec;
6475 }
6476 \f
6477 static struct matrix_cmd *
6478 matrix_parse_setdiag (struct matrix_state *s)
6479 {
6480   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6481   *cmd = (struct matrix_cmd) {
6482     .type = MCMD_SETDIAG,
6483     .setdiag = { .dst = NULL }
6484   };
6485
6486   struct setdiag_command *setdiag = &cmd->setdiag;
6487   if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6488     goto error;
6489
6490   setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6491   if (!setdiag->dst)
6492     {
6493       lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6494       goto error;
6495     }
6496   lex_get (s->lexer);
6497
6498   if (!lex_force_match (s->lexer, T_COMMA))
6499     goto error;
6500
6501   setdiag->expr = matrix_parse_expr (s);
6502   if (!setdiag->expr)
6503     goto error;
6504
6505   if (!lex_force_match (s->lexer, T_RPAREN))
6506     goto error;
6507
6508   return cmd;
6509
6510 error:
6511   matrix_cmd_destroy (cmd);
6512   return NULL;
6513 }
6514
6515 static void
6516 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6517 {
6518   gsl_matrix *dst = setdiag->dst->value;
6519   if (!dst)
6520     {
6521       msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6522            setdiag->dst->name);
6523       return;
6524     }
6525
6526   gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6527   if (!src)
6528     return;
6529
6530   size_t n = MIN (dst->size1, dst->size2);
6531   if (is_scalar (src))
6532     {
6533       double d = to_scalar (src);
6534       for (size_t i = 0; i < n; i++)
6535         gsl_matrix_set (dst, i, i, d);
6536     }
6537   else if (is_vector (src))
6538     {
6539       gsl_vector v = to_vector (src);
6540       for (size_t i = 0; i < n && i < v.size; i++)
6541         gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6542     }
6543   else
6544     msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6545                "not a %zu×%zu matrix."),
6546          src->size1, src->size2);
6547   gsl_matrix_free (src);
6548 }
6549 \f
6550 static struct matrix_cmd *
6551 matrix_parse_svd (struct matrix_state *s)
6552 {
6553   struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6554   *cmd = (struct matrix_cmd) {
6555     .type = MCMD_SVD,
6556     .svd = { .expr = NULL }
6557   };
6558
6559   struct svd_command *svd = &cmd->svd;
6560   if (!lex_force_match (s->lexer, T_LPAREN))
6561     goto error;
6562   svd->expr = matrix_parse_expr (s);
6563   if (!svd->expr
6564       || !lex_force_match (s->lexer, T_COMMA)
6565       || !matrix_parse_dst_var (s, &svd->u)
6566       || !lex_force_match (s->lexer, T_COMMA)
6567       || !matrix_parse_dst_var (s, &svd->s)
6568       || !lex_force_match (s->lexer, T_COMMA)
6569       || !matrix_parse_dst_var (s, &svd->v)
6570       || !lex_force_match (s->lexer, T_RPAREN))
6571     goto error;
6572
6573   return cmd;
6574
6575 error:
6576   matrix_cmd_destroy (cmd);
6577   return NULL;
6578 }
6579
6580 static void
6581 matrix_cmd_execute_svd (struct svd_command *svd)
6582 {
6583   gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6584   if (!m)
6585     return;
6586
6587   if (m->size1 >= m->size2)
6588     {
6589       gsl_matrix *A = m;
6590       gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6591       gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6592       gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6593       gsl_vector *work = gsl_vector_alloc (A->size2);
6594       gsl_linalg_SV_decomp (A, V, &Sv, work);
6595       gsl_vector_free (work);
6596
6597       matrix_var_set (svd->u, A);
6598       matrix_var_set (svd->s, S);
6599       matrix_var_set (svd->v, V);
6600     }
6601   else
6602     {
6603       gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6604       gsl_matrix_transpose_memcpy (At, m);
6605       gsl_matrix_free (m);
6606
6607       gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6608       gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6609       gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6610       gsl_vector *work = gsl_vector_alloc (At->size2);
6611       gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6612       gsl_vector_free (work);
6613
6614       matrix_var_set (svd->v, At);
6615       matrix_var_set (svd->s, St);
6616       matrix_var_set (svd->u, Vt);
6617     }
6618 }
6619 \f
6620 static bool
6621 matrix_cmd_execute (struct matrix_cmd *cmd)
6622 {
6623   switch (cmd->type)
6624     {
6625     case MCMD_COMPUTE:
6626       matrix_cmd_execute_compute (&cmd->compute);
6627       break;
6628
6629     case MCMD_PRINT:
6630       matrix_cmd_execute_print (&cmd->print);
6631       break;
6632
6633     case MCMD_DO_IF:
6634       return matrix_cmd_execute_do_if (&cmd->do_if);
6635
6636     case MCMD_LOOP:
6637       matrix_cmd_execute_loop (&cmd->loop);
6638       break;
6639
6640     case MCMD_BREAK:
6641       return false;
6642
6643     case MCMD_DISPLAY:
6644       matrix_cmd_execute_display (&cmd->display);
6645       break;
6646
6647     case MCMD_RELEASE:
6648       matrix_cmd_execute_release (&cmd->release);
6649       break;
6650
6651     case MCMD_SAVE:
6652       matrix_cmd_execute_save (&cmd->save);
6653       break;
6654
6655     case MCMD_READ:
6656       matrix_cmd_execute_read (&cmd->read);
6657       break;
6658
6659     case MCMD_WRITE:
6660       matrix_cmd_execute_write (&cmd->write);
6661       break;
6662
6663     case MCMD_GET:
6664       matrix_cmd_execute_get (&cmd->get);
6665       break;
6666
6667     case MCMD_MSAVE:
6668       matrix_cmd_execute_msave (&cmd->msave);
6669       break;
6670
6671     case MCMD_MGET:
6672       matrix_cmd_execute_mget (&cmd->mget);
6673       break;
6674
6675     case MCMD_EIGEN:
6676       matrix_cmd_execute_eigen (&cmd->eigen);
6677       break;
6678
6679     case MCMD_SETDIAG:
6680       matrix_cmd_execute_setdiag (&cmd->setdiag);
6681       break;
6682
6683     case MCMD_SVD:
6684       matrix_cmd_execute_svd (&cmd->svd);
6685       break;
6686     }
6687
6688   return true;
6689 }
6690
6691 static void
6692 matrix_cmds_uninit (struct matrix_cmds *cmds)
6693 {
6694   for (size_t i = 0; i < cmds->n; i++)
6695     matrix_cmd_destroy (cmds->commands[i]);
6696   free (cmds->commands);
6697 }
6698
6699 static void
6700 matrix_cmd_destroy (struct matrix_cmd *cmd)
6701 {
6702   if (!cmd)
6703     return;
6704
6705   switch (cmd->type)
6706     {
6707     case MCMD_COMPUTE:
6708       matrix_lvalue_destroy (cmd->compute.lvalue);
6709       matrix_expr_destroy (cmd->compute.rvalue);
6710       break;
6711
6712     case MCMD_PRINT:
6713       matrix_expr_destroy (cmd->print.expression);
6714       free (cmd->print.title);
6715       print_labels_destroy (cmd->print.rlabels);
6716       print_labels_destroy (cmd->print.clabels);
6717       break;
6718
6719     case MCMD_DO_IF:
6720       for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6721         {
6722           matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6723           matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6724         }
6725       free (cmd->do_if.clauses);
6726       break;
6727
6728     case MCMD_LOOP:
6729       matrix_expr_destroy (cmd->loop.start);
6730       matrix_expr_destroy (cmd->loop.end);
6731       matrix_expr_destroy (cmd->loop.increment);
6732       matrix_expr_destroy (cmd->loop.top_condition);
6733       matrix_expr_destroy (cmd->loop.bottom_condition);
6734       matrix_cmds_uninit (&cmd->loop.commands);
6735       break;
6736
6737     case MCMD_BREAK:
6738       break;
6739
6740     case MCMD_DISPLAY:
6741       break;
6742
6743     case MCMD_RELEASE:
6744       free (cmd->release.vars);
6745       break;
6746
6747     case MCMD_SAVE:
6748       matrix_expr_destroy (cmd->save.expression);
6749       break;
6750
6751     case MCMD_READ:
6752       matrix_lvalue_destroy (cmd->read.dst);
6753       matrix_expr_destroy (cmd->read.size);
6754       break;
6755
6756     case MCMD_WRITE:
6757       matrix_expr_destroy (cmd->write.expression);
6758       free (cmd->write.format);
6759       break;
6760
6761     case MCMD_GET:
6762       matrix_lvalue_destroy (cmd->get.dst);
6763       fh_unref (cmd->get.file);
6764       free (cmd->get.encoding);
6765       var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6766       break;
6767
6768     case MCMD_MSAVE:
6769       free (cmd->msave.varname_);
6770       matrix_expr_destroy (cmd->msave.expr);
6771       matrix_expr_destroy (cmd->msave.factors);
6772       matrix_expr_destroy (cmd->msave.splits);
6773       break;
6774
6775     case MCMD_MGET:
6776       fh_unref (cmd->mget.file);
6777       stringi_set_destroy (&cmd->mget.rowtypes);
6778       break;
6779
6780     case MCMD_EIGEN:
6781       matrix_expr_destroy (cmd->eigen.expr);
6782       break;
6783
6784     case MCMD_SETDIAG:
6785       matrix_expr_destroy (cmd->setdiag.expr);
6786       break;
6787
6788     case MCMD_SVD:
6789       matrix_expr_destroy (cmd->svd.expr);
6790       break;
6791     }
6792   free (cmd);
6793 }
6794
6795 struct matrix_command_name
6796   {
6797     const char *name;
6798     struct matrix_cmd *(*parse) (struct matrix_state *);
6799   };
6800
6801 static const struct matrix_command_name *
6802 matrix_parse_command_name (struct lexer *lexer)
6803 {
6804   static const struct matrix_command_name commands[] = {
6805     { "COMPUTE", matrix_parse_compute },
6806     { "CALL EIGEN", matrix_parse_eigen },
6807     { "CALL SETDIAG", matrix_parse_setdiag },
6808     { "CALL SVD", matrix_parse_svd },
6809     { "PRINT", matrix_parse_print },
6810     { "DO IF", matrix_parse_do_if },
6811     { "LOOP", matrix_parse_loop },
6812     { "BREAK", matrix_parse_break },
6813     { "READ", matrix_parse_read },
6814     { "WRITE", matrix_parse_write },
6815     { "GET", matrix_parse_get },
6816     { "SAVE", matrix_parse_save },
6817     { "MGET", matrix_parse_mget },
6818     { "MSAVE", matrix_parse_msave },
6819     { "DISPLAY", matrix_parse_display },
6820     { "RELEASE", matrix_parse_release },
6821   };
6822   static size_t n = sizeof commands / sizeof *commands;
6823
6824   for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6825     if (lex_match_phrase (lexer, c->name))
6826       return c;
6827   return NULL;
6828 }
6829
6830 static struct matrix_cmd *
6831 matrix_parse_command (struct matrix_state *s)
6832 {
6833   size_t nesting_level = SIZE_MAX;
6834
6835   struct matrix_cmd *c = NULL;
6836   const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6837   if (!cmd)
6838     lex_error (s->lexer, _("Unknown matrix command."));
6839   else if (!cmd->parse)
6840     lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6841                cmd->name);
6842   else
6843     {
6844       nesting_level = output_open_group (
6845         group_item_create_nocopy (utf8_to_title (cmd->name),
6846                                   utf8_to_title (cmd->name)));
6847       c = cmd->parse (s);
6848     }
6849
6850   if (c)
6851     lex_end_of_command (s->lexer);
6852   lex_discard_rest_of_command (s->lexer);
6853   if (nesting_level != SIZE_MAX)
6854     output_close_groups (nesting_level);
6855
6856   return c;
6857 }
6858
6859 int
6860 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6861 {
6862   if (!lex_force_match (lexer, T_ENDCMD))
6863     return CMD_FAILURE;
6864
6865   struct matrix_state state = {
6866     .dataset = ds,
6867     .session = dataset_session (ds),
6868     .lexer = lexer,
6869     .vars = HMAP_INITIALIZER (state.vars),
6870   };
6871
6872   for (;;)
6873     {
6874       while (lex_match (lexer, T_ENDCMD))
6875         continue;
6876       if (lex_token (lexer) == T_STOP)
6877         {
6878           msg (SE, _("Unexpected end of input expecting matrix command."));
6879           break;
6880         }
6881
6882       if (lex_match_phrase (lexer, "END MATRIX"))
6883         break;
6884
6885       struct matrix_cmd *c = matrix_parse_command (&state);
6886       if (c)
6887         {
6888           matrix_cmd_execute (c);
6889           matrix_cmd_destroy (c);
6890         }
6891     }
6892
6893   struct matrix_var *var, *next;
6894   HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6895     {
6896       free (var->name);
6897       gsl_matrix_free (var->value);
6898       hmap_delete (&state.vars, &var->hmap_node);
6899       free (var);
6900     }
6901   hmap_destroy (&state.vars);
6902   if (state.common)
6903     {
6904       dict_unref (state.common->dict);
6905       casewriter_destroy (state.common->writer);
6906       free (state.common);
6907     }
6908   fh_unref (state.prev_read_file);
6909   for (size_t i = 0; i < state.n_read_files; i++)
6910     read_file_destroy (state.read_files[i]);
6911   free (state.read_files);
6912   fh_unref (state.prev_write_file);
6913   for (size_t i = 0; i < state.n_write_files; i++)
6914     write_file_destroy (state.write_files[i]);
6915   free (state.write_files);
6916   fh_unref (state.prev_save_file);
6917   for (size_t i = 0; i < state.n_save_files; i++)
6918     save_file_destroy (state.save_files[i]);
6919   free (state.save_files);
6920
6921   return CMD_SUCCESS;
6922 }