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

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