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