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

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