pintos-os.org Git - pspp/blob - src/language/stats/matrix.c

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