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