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