--- /dev/null
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2021 Free Software Foundation, Inc.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include <config.h>
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_cdf.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <limits.h>
+#include <math.h>
+#include <uniwidth.h>
+
+#include "data/any-reader.h"
+#include "data/any-writer.h"
+#include "data/casereader.h"
+#include "data/casewriter.h"
+#include "data/data-in.h"
+#include "data/data-out.h"
+#include "data/dataset.h"
+#include "data/dictionary.h"
+#include "data/file-handle-def.h"
+#include "language/command.h"
+#include "language/commands/data-reader.h"
+#include "language/commands/data-writer.h"
+#include "language/commands/file-handle.h"
+#include "language/lexer/format-parser.h"
+#include "language/lexer/lexer.h"
+#include "language/lexer/variable-parser.h"
+#include "libpspp/array.h"
+#include "libpspp/assertion.h"
+#include "libpspp/compiler.h"
+#include "libpspp/hmap.h"
+#include "libpspp/i18n.h"
+#include "libpspp/intern.h"
+#include "libpspp/misc.h"
+#include "libpspp/str.h"
+#include "libpspp/string-array.h"
+#include "libpspp/stringi-set.h"
+#include "libpspp/u8-line.h"
+#include "math/distributions.h"
+#include "math/random.h"
+#include "output/driver.h"
+#include "output/output-item.h"
+#include "output/pivot-table.h"
+
+#include "gl/c-ctype.h"
+#include "gl/c-strcase.h"
+#include "gl/ftoastr.h"
+#include "gl/minmax.h"
+#include "gl/xsize.h"
+
+#include "gettext.h"
+#define _(msgid) gettext (msgid)
+#define N_(msgid) (msgid)
+
+struct matrix_state;
+
+/* A variable in the matrix language. */
+struct matrix_var
+ {
+ struct hmap_node hmap_node; /* In matrix_state's 'vars' hmap. */
+ char *name; /* UTF-8. */
+ gsl_matrix *value; /* NULL, if the variable is uninitialized. */
+ };
+
+/* All the MSAVE commands within a matrix program share common configuration,
+ provided by the first MSAVE command within the program. This structure
+ encapsulates this configuration. */
+struct msave_common
+ {
+ /* Common configuration for all MSAVEs. */
+ struct msg_location *location; /* Range of lines for first MSAVE. */
+ struct file_handle *outfile; /* Output file for all the MSAVEs. */
+ struct msg_location *outfile_location;
+ struct string_array variables; /* VARIABLES subcommand. */
+ struct msg_location *variables_location;
+ struct string_array fnames; /* FNAMES subcommand. */
+ struct msg_location *fnames_location;
+ struct string_array snames; /* SNAMES subcommand. */
+ struct msg_location *snames_location;
+
+ /* Collects and owns factors and splits. The individual msave_command
+ structs point to these but do not own them. (This is because factors
+ and splits can be carried over from one MSAVE to the next, so it's
+ easiest to just take the most recent.) */
+ struct matrix_expr **factors;
+ size_t n_factors, allocated_factors;
+ struct matrix_expr **splits;
+ size_t n_splits, allocated_splits;
+
+ /* Execution state. */
+ struct dictionary *dict;
+ struct casewriter *writer;
+ };
+
+/* A file used by one or more READ commands. */
+struct read_file
+ {
+ /* Parse state. */
+ struct file_handle *file;
+
+ /* Execution state. */
+ struct dfm_reader *reader;
+ char *encoding;
+ };
+
+static struct read_file *read_file_create (struct matrix_state *,
+ struct file_handle *);
+static struct dfm_reader *read_file_open (struct read_file *);
+
+/* A file used by one or more WRITE comamnds. */
+struct write_file
+ {
+ /* Parse state. */
+ struct file_handle *file;
+
+ /* Execution state. */
+ struct dfm_writer *writer;
+ char *encoding;
+ struct u8_line *held; /* Output held by a previous WRITE with /HOLD. */
+ };
+
+static struct write_file *write_file_create (struct matrix_state *,
+ struct file_handle *);
+static struct dfm_writer *write_file_open (struct write_file *);
+static void write_file_destroy (struct write_file *);
+
+/* A file used by one or more SAVE commands. */
+struct save_file
+ {
+ /* Parse state. */
+ struct file_handle *file;
+ struct dataset *dataset;
+ struct string_array variables;
+ struct matrix_expr *names;
+ struct stringi_set strings;
+
+ /* Execution state. */
+ bool error;
+ struct casewriter *writer;
+ struct dictionary *dict;
+ struct msg_location *location;
+ };
+
+/* State of an entire matrix program. */
+struct matrix_state
+ {
+ /* State passed into MATRIX from outside. */
+ struct dataset *dataset;
+ struct session *session;
+ struct lexer *lexer;
+
+ /* Matrix program's own state. */
+ struct hmap vars; /* Dictionary of matrix variables. */
+ bool in_loop; /* True if parsing within a LOOP. */
+
+ /* MSAVE. */
+ struct msave_common *msave_common;
+
+ /* READ. */
+ struct file_handle *prev_read_file;
+ struct read_file **read_files;
+ size_t n_read_files;
+
+ /* WRITE. */
+ struct file_handle *prev_write_file;
+ struct write_file **write_files;
+ size_t n_write_files;
+
+ /* SAVE. */
+ struct file_handle *prev_save_file;
+ struct save_file **save_files;
+ size_t n_save_files;
+ };
+
+/* Finds and returns the variable with the given NAME (case-insensitive) within
+ S, if there is one, or a null pointer if there is not. */
+static struct matrix_var *
+matrix_var_lookup (struct matrix_state *s, struct substring name)
+{
+ struct matrix_var *var;
+
+ HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
+ utf8_hash_case_substring (name, 0), &s->vars)
+ if (!utf8_sscasecmp (ss_cstr (var->name), name))
+ return var;
+
+ return NULL;
+}
+
+/* Creates and returns a new variable named NAME within S. There must not
+ already be a variable with the same (case-insensitive) name. The variable
+ is created uninitialized. */
+static struct matrix_var *
+matrix_var_create (struct matrix_state *s, struct substring name)
+{
+ struct matrix_var *var = xmalloc (sizeof *var);
+ *var = (struct matrix_var) { .name = ss_xstrdup (name) };
+ hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
+ return var;
+}
+
+/* Replaces VAR's value by VALUE. Takes ownership of VALUE. */
+static void
+matrix_var_set (struct matrix_var *var, gsl_matrix *value)
+{
+ gsl_matrix_free (var->value);
+ var->value = value;
+}
+\f
+/* Matrix function catalog. */
+
+/* The third argument to F() is a "prototype". For most prototypes, the first
+ letter (before the _) represents the return type and each other letter
+ (after the _) is an argument type. The types are:
+
+ - "m": A matrix of unrestricted dimensions.
+
+ - "d": A scalar.
+
+ - "v": A row or column vector.
+
+ - "e": Primarily for the first argument, this is a matrix with
+ unrestricted dimensions treated elementwise. Each element in the matrix
+ is passed to the implementation function separately.
+
+ - "n": This gets passed the "const struct matrix_expr *" that represents
+ the expression. This allows the evaluation function to grab the source
+ location of arguments so that it can report accurate error locations.
+ This type doesn't correspond to an argument passed in by the user.
+
+ The fourth argument is an optional constraints string. For this purpose the
+ first argument is named "a", the second "b", and so on. The following kinds
+ of constraints are supported. For matrix arguments, the constraints are
+ applied to each value in the matrix separately:
+
+ - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
+ integer may substitute for 0 and 1. Half-open constraints (] and [) are
+ also supported.
+
+ - "ai": Restrict a to integer values.
+
+ - "a>0", "a<0", "a>=0", "a<=0", "a!=0".
+
+ - "a<b", "a>b", "a<=b", "a>=b", "b!=0".
+*/
+#define MATRIX_FUNCTIONS \
+ F(ABS, "ABS", m_e, NULL) \
+ F(ALL, "ALL", d_m, NULL) \
+ F(ANY, "ANY", d_m, NULL) \
+ F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
+ F(ARTAN, "ARTAN", m_e, NULL) \
+ F(BLOCK, "BLOCK", m_any, NULL) \
+ F(CHOL, "CHOL", m_mn, NULL) \
+ F(CMIN, "CMIN", m_m, NULL) \
+ F(CMAX, "CMAX", m_m, NULL) \
+ F(COS, "COS", m_e, NULL) \
+ F(CSSQ, "CSSQ", m_m, NULL) \
+ F(CSUM, "CSUM", m_m, NULL) \
+ F(DESIGN, "DESIGN", m_mn, NULL) \
+ F(DET, "DET", d_m, NULL) \
+ F(DIAG, "DIAG", m_m, NULL) \
+ F(EVAL, "EVAL", m_mn, NULL) \
+ F(EXP, "EXP", m_e, NULL) \
+ F(GINV, "GINV", m_m, NULL) \
+ F(GRADE, "GRADE", m_m, NULL) \
+ F(GSCH, "GSCH", m_mn, NULL) \
+ F(IDENT, "IDENT", IDENT, NULL) \
+ F(INV, "INV", m_m, NULL) \
+ F(KRONEKER, "KRONEKER", m_mm, NULL) \
+ F(LG10, "LG10", m_e, "a>0") \
+ F(LN, "LN", m_e, "a>0") \
+ F(MAGIC, "MAGIC", m_d, "ai>=3") \
+ F(MAKE, "MAKE", m_ddd, "ai>=0 bi>=0") \
+ F(MDIAG, "MDIAG", m_v, NULL) \
+ F(MMAX, "MMAX", d_m, NULL) \
+ F(MMIN, "MMIN", d_m, NULL) \
+ F(MOD, "MOD", m_md, "b!=0") \
+ F(MSSQ, "MSSQ", d_m, NULL) \
+ F(MSUM, "MSUM", d_m, NULL) \
+ F(NCOL, "NCOL", d_m, NULL) \
+ F(NROW, "NROW", d_m, NULL) \
+ F(RANK, "RANK", d_m, NULL) \
+ F(RESHAPE, "RESHAPE", m_mddn, NULL) \
+ F(RMAX, "RMAX", m_m, NULL) \
+ F(RMIN, "RMIN", m_m, NULL) \
+ F(RND, "RND", m_e, NULL) \
+ F(RNKORDER, "RNKORDER", m_m, NULL) \
+ F(RSSQ, "RSSQ", m_m, NULL) \
+ F(RSUM, "RSUM", m_m, NULL) \
+ F(SIN, "SIN", m_e, NULL) \
+ F(SOLVE, "SOLVE", m_mmn, NULL) \
+ F(SQRT, "SQRT", m_e, "a>=0") \
+ F(SSCP, "SSCP", m_m, NULL) \
+ F(SVAL, "SVAL", m_m, NULL) \
+ F(SWEEP, "SWEEP", m_mdn, NULL) \
+ F(T, "T", m_m, NULL) \
+ F(TRACE, "TRACE", d_m, NULL) \
+ F(TRANSPOS, "TRANSPOS", m_m, NULL) \
+ F(TRUNC, "TRUNC", m_e, NULL) \
+ F(UNIFORM, "UNIFORM", m_ddn, "ai>=0 bi>=0") \
+ F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
+ F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
+ F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
+ F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
+ F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
+ F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
+ F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
+ F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
+ F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
+ F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
+ F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
+ F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
+ F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
+ F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
+ F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
+ F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
+ F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
+ F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
+ F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
+ F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
+ F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
+ F(RV_EXP, "RV.EXP", d_d, "a>0") \
+ F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
+ F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
+ F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
+ F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
+ F(RV_F, "RV.F", d_dd, "a>0 b>0") \
+ F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
+ F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
+ F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
+ F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
+ F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
+ F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
+ F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
+ F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
+ F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
+ F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
+ F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
+ F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
+ F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
+ F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
+ F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
+ F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
+ F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
+ F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
+ F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
+ F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
+ F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
+ F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
+ F(CDFNORM, "CDFNORM", m_e, NULL) \
+ F(PROBIT, "PROBIT", m_e, "a(0,1)") \
+ F(NORMAL, "NORMAL", m_e, "a>0") \
+ F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
+ F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
+ F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
+ F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
+ F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
+ F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
+ F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
+ F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
+ F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
+ F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
+ F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
+ F(CDF_T, "CDF.T", m_ed, "b>0") \
+ F(TCDF, "TCDF", m_ed, "b>0") \
+ F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
+ F(PDF_T, "PDF.T", m_ed, "b>0") \
+ F(RV_T, "RV.T", d_d, "a>0") \
+ F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
+ F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
+ F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
+ F(RV_T1G, "RV.T1G", d_dd, NULL) \
+ F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
+ F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
+ F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
+ F(RV_T2G, "RV.T2G", d_dd, NULL) \
+ F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
+ F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
+ F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
+ F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
+ F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
+ F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
+ F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
+ F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
+ F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
+ F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
+ F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
+ F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
+ F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
+ F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
+ F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
+ F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
+ F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
+ F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
+ F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
+ F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
+ F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
+ F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
+ F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
+ F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
+ F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
+ F(RV_POISSON, "RV.POISSON", d_d, "a>0")
+
+/* Properties of a matrix function.
+
+ These come straight from the macro invocations above. */
+struct matrix_function_properties
+ {
+ const char *name;
+ const char *constraints;
+ };
+
+/* Minimum and maximum argument counts for each matrix function prototype. */
+enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
+enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
+enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
+enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
+enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
+enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
+enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
+enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
+enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
+enum { m_ddn_MIN_ARGS = 2, m_ddn_MAX_ARGS = 2 };
+enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
+enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
+enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
+enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
+enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
+enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
+enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
+enum { m_mddn_MIN_ARGS = 3, m_mddn_MAX_ARGS = 3 };
+enum { m_mdn_MIN_ARGS = 2, m_mdn_MAX_ARGS = 2 };
+enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
+enum { m_mmn_MIN_ARGS = 2, m_mmn_MAX_ARGS = 2 };
+enum { m_mn_MIN_ARGS = 1, m_mn_MAX_ARGS = 1 };
+enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
+
+/* C function prototype for each matrix function prototype. */
+typedef double matrix_proto_d_none (void);
+typedef double matrix_proto_d_d (double);
+typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_ddd (double, double, double);
+typedef gsl_matrix *matrix_proto_m_d (double);
+typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
+typedef gsl_matrix *matrix_proto_m_ddn (double, double,
+ const struct matrix_expr *);
+typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
+typedef gsl_matrix *matrix_proto_m_mn (gsl_matrix *,
+ const struct matrix_expr *);
+typedef double matrix_proto_m_e (double);
+typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
+typedef gsl_matrix *matrix_proto_m_mdn (gsl_matrix *, double,
+ const struct matrix_expr *);
+typedef double matrix_proto_m_ed (double, double);
+typedef gsl_matrix *matrix_proto_m_mddn (gsl_matrix *, double, double,
+ const struct matrix_expr *);
+typedef double matrix_proto_m_edd (double, double, double);
+typedef double matrix_proto_m_eddd (double, double, double, double);
+typedef double matrix_proto_m_eed (double, double, double);
+typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
+typedef gsl_matrix *matrix_proto_m_mmn (gsl_matrix *, gsl_matrix *,
+ const struct matrix_expr *);
+typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
+typedef double matrix_proto_d_m (gsl_matrix *);
+typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
+typedef gsl_matrix *matrix_proto_IDENT (double, double);
+
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+ static matrix_proto_##PROTO matrix_eval_##ENUM;
+MATRIX_FUNCTIONS
+#undef F
+\f
+/* Matrix expression data structure and parsing. */
+
+/* A node in a matrix expression. */
+struct matrix_expr
+ {
+ enum matrix_op
+ {
+ /* Functions. */
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
+ MATRIX_FUNCTIONS
+#undef F
+
+ /* Elementwise and scalar arithmetic. */
+ MOP_NEGATE, /* unary - */
+ MOP_ADD_ELEMS, /* + */
+ MOP_SUB_ELEMS, /* - */
+ MOP_MUL_ELEMS, /* &* */
+ MOP_DIV_ELEMS, /* / and &/ */
+ MOP_EXP_ELEMS, /* &** */
+ MOP_SEQ, /* a:b */
+ MOP_SEQ_BY, /* a:b:c */
+
+ /* Matrix arithmetic. */
+ MOP_MUL_MAT, /* * */
+ MOP_EXP_MAT, /* ** */
+
+ /* Relational. */
+ MOP_GT, /* > */
+ MOP_GE, /* >= */
+ MOP_LT, /* < */
+ MOP_LE, /* <= */
+ MOP_EQ, /* = */
+ MOP_NE, /* <> */
+
+ /* Logical. */
+ MOP_NOT, /* NOT */
+ MOP_AND, /* AND */
+ MOP_OR, /* OR */
+ MOP_XOR, /* XOR */
+
+ /* {}. */
+ MOP_PASTE_HORZ, /* a, b, c, ... */
+ MOP_PASTE_VERT, /* a; b; c; ... */
+ MOP_EMPTY, /* {} */
+
+ /* Sub-matrices. */
+ MOP_VEC_INDEX, /* x(y) */
+ MOP_VEC_ALL, /* x(:) */
+ MOP_MAT_INDEX, /* x(y,z) */
+ MOP_ROW_INDEX, /* x(y,:) */
+ MOP_COL_INDEX, /* x(:,z) */
+
+ /* Literals. */
+ MOP_NUMBER,
+ MOP_VARIABLE,
+
+ /* Oddball stuff. */
+ MOP_EOF, /* EOF('file') */
+ }
+ op;
+
+ union
+ {
+ /* Nonterminal expression nodes. */
+ struct
+ {
+ struct matrix_expr **subs;
+ size_t n_subs;
+ };
+
+ /* Terminal expression nodes. */
+ double number; /* MOP_NUMBER. */
+ struct matrix_var *variable; /* MOP_VARIABLE. */
+ struct read_file *eof; /* MOP_EOF. */
+ };
+
+ /* The syntax location corresponding to this expression node, for use in
+ error messages. This is always nonnull for terminal expression nodes.
+ For most others, it is null because it can be computed lazily if and
+ when it is needed.
+
+ Use matrix_expr_location() instead of using this member directly, so
+ that it gets computed lazily if needed. */
+ struct msg_location *location;
+ };
+
+static void
+matrix_expr_location__ (const struct matrix_expr *e,
+ const struct msg_location **minp,
+ const struct msg_location **maxp)
+{
+ struct msg_location *loc = e->location;
+ if (loc)
+ {
+ const struct msg_location *min = *minp;
+ if (loc->start.line
+ && (!min
+ || loc->start.line < min->start.line
+ || (loc->start.line == min->start.line
+ && loc->start.column < min->start.column)))
+ *minp = loc;
+
+ const struct msg_location *max = *maxp;
+ if (loc->end.line
+ && (!max
+ || loc->end.line > max->end.line
+ || (loc->end.line == max->end.line
+ && loc->end.column > max->end.column)))
+ *maxp = loc;
+
+ return;
+ }
+
+ assert (e->op != MOP_NUMBER && e->op != MOP_VARIABLE && e->op != MOP_EOF);
+ for (size_t i = 0; i < e->n_subs; i++)
+ matrix_expr_location__ (e->subs[i], minp, maxp);
+}
+
+/* Returns the source code location corresponding to expression E, computing it
+ lazily if needed. */
+static const struct msg_location *
+matrix_expr_location (const struct matrix_expr *e_)
+{
+ struct matrix_expr *e = CONST_CAST (struct matrix_expr *, e_);
+ if (!e)
+ return NULL;
+
+ if (!e->location)
+ {
+ const struct msg_location *min = NULL;
+ const struct msg_location *max = NULL;
+ matrix_expr_location__ (e, &min, &max);
+ if (min && max)
+ {
+ e->location = msg_location_dup (min);
+ e->location->end = max->end;
+ }
+ }
+ return e->location;
+}
+
+/* Sets e->location to the tokens in S's lexer from offset START_OFS to the
+ token before the current one. Has no effect if E already has a location or
+ if E is null. */
+static void
+matrix_expr_add_location (struct matrix_state *s, int start_ofs,
+ struct matrix_expr *e)
+{
+ if (e && !e->location)
+ e->location = lex_ofs_location (s->lexer, start_ofs,
+ lex_ofs (s->lexer) - 1);
+}
+
+/* Frees E and all the data and sub-expressions that it references. */
+static void
+matrix_expr_destroy (struct matrix_expr *e)
+{
+ if (!e)
+ return;
+
+ switch (e->op)
+ {
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
+MATRIX_FUNCTIONS
+#undef F
+ case MOP_NEGATE:
+ case MOP_ADD_ELEMS:
+ case MOP_SUB_ELEMS:
+ case MOP_MUL_ELEMS:
+ case MOP_DIV_ELEMS:
+ case MOP_EXP_ELEMS:
+ case MOP_SEQ:
+ case MOP_SEQ_BY:
+ case MOP_MUL_MAT:
+ case MOP_EXP_MAT:
+ case MOP_GT:
+ case MOP_GE:
+ case MOP_LT:
+ case MOP_LE:
+ case MOP_EQ:
+ case MOP_NE:
+ case MOP_NOT:
+ case MOP_AND:
+ case MOP_OR:
+ case MOP_XOR:
+ case MOP_EMPTY:
+ case MOP_PASTE_HORZ:
+ case MOP_PASTE_VERT:
+ case MOP_VEC_INDEX:
+ case MOP_VEC_ALL:
+ case MOP_MAT_INDEX:
+ case MOP_ROW_INDEX:
+ case MOP_COL_INDEX:
+ for (size_t i = 0; i < e->n_subs; i++)
+ matrix_expr_destroy (e->subs[i]);
+ free (e->subs);
+ break;
+
+ case MOP_NUMBER:
+ case MOP_VARIABLE:
+ case MOP_EOF:
+ break;
+ }
+ msg_location_destroy (e->location);
+ free (e);
+}
+
+/* Creates and returns a new matrix_expr with type OP, which must be a
+ nonterminal type. Initializes the new matrix_expr with the N_SUBS
+ expressions in SUBS as subexpressions. */
+static struct matrix_expr *
+matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
+ size_t n_subs)
+{
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) {
+ .op = op,
+ .subs = xmemdup (subs, n_subs * sizeof *subs),
+ .n_subs = n_subs
+ };
+ return e;
+}
+
+static struct matrix_expr *
+matrix_expr_create_0 (enum matrix_op op)
+{
+ struct matrix_expr *sub;
+ return matrix_expr_create_subs (op, &sub, 0);
+}
+
+static struct matrix_expr *
+matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
+{
+ return matrix_expr_create_subs (op, &sub, 1);
+}
+
+static struct matrix_expr *
+matrix_expr_create_2 (enum matrix_op op,
+ struct matrix_expr *sub0, struct matrix_expr *sub1)
+{
+ struct matrix_expr *subs[] = { sub0, sub1 };
+ return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
+}
+
+static struct matrix_expr *
+matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
+ struct matrix_expr *sub1, struct matrix_expr *sub2)
+{
+ struct matrix_expr *subs[] = { sub0, sub1, sub2 };
+ return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
+}
+
+/* Creates and returns a new MOP_NUMBER expression node to contain NUMBER. */
+static struct matrix_expr *
+matrix_expr_create_number (double number)
+{
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) {
+ .op = MOP_NUMBER,
+ .number = number,
+ };
+ return e;
+}
+
+static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
+
+/* A binary operator for matrix_parse_binary_operator(). */
+struct matrix_operator_syntax
+ {
+ /* Exactly one of these specifies the operator syntax. */
+ enum token_type token; /* A token, e.g. T_ASTERISK. */
+ const char *id; /* An identifier, e.g. "XOR". */
+ const char *phrase; /* A token phrase, e.g. "&**". */
+
+ /* The matrix operator corresponding to the syntax. */
+ enum matrix_op op;
+ };
+
+static bool
+matrix_operator_syntax_match (struct lexer *lexer,
+ const struct matrix_operator_syntax *syntax,
+ size_t n_syntax, enum matrix_op *op)
+{
+ const struct matrix_operator_syntax *end = &syntax[n_syntax];
+ for (const struct matrix_operator_syntax *syn = syntax; syn < end; syn++)
+ if (syn->id ? lex_match_id (lexer, syn->id)
+ : syn->phrase ? lex_match_phrase (lexer, syn->phrase)
+ : lex_match (lexer, syn->token))
+ {
+ *op = syn->op;
+ return true;
+ }
+ return false;
+}
+
+/* Parses a binary operator level in the recursive descent parser, returning a
+ matrix expression if successful or a null pointer otherwise. PARSE_NEXT
+ must be the function to parse the next level of precedence. The N_SYNTAX
+ elements of SYNTAX must specify the syntax and matrix_expr node type to
+ parse at this level. */
+static struct matrix_expr *
+matrix_parse_binary_operator (
+ struct matrix_state *s,
+ struct matrix_expr *(*parse_next) (struct matrix_state *),
+ const struct matrix_operator_syntax *syntax, size_t n_syntax)
+{
+ struct matrix_expr *lhs = parse_next (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (!matrix_operator_syntax_match (s->lexer, syntax, n_syntax, &op))
+ return lhs;
+
+ struct matrix_expr *rhs = parse_next (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+/* Parses a comma-separated list of expressions within {}, transforming them
+ into MOP_PASTE_HORZ operators. Returns the new expression or NULL on
+ error. */
+static struct matrix_expr *
+matrix_parse_curly_comma (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax op = {
+ .token = T_COMMA, .op = MOP_PASTE_HORZ
+ };
+ return matrix_parse_binary_operator (s, matrix_expr_parse, &op, 1);
+}
+
+/* Parses a semicolon-separated list of expressions within {}, transforming
+ them into MOP_PASTE_VERT operators. Returns the new expression or NULL on
+ error. */
+static struct matrix_expr *
+matrix_parse_curly_semi (struct matrix_state *s)
+{
+ if (lex_token (s->lexer) == T_RCURLY)
+ {
+ /* {} is a special case for a 0×0 matrix. */
+ return matrix_expr_create_0 (MOP_EMPTY);
+ }
+
+ static const struct matrix_operator_syntax op = {
+ .token = T_SEMICOLON, .op = MOP_PASTE_VERT
+ };
+ return matrix_parse_binary_operator (s, matrix_parse_curly_comma, &op, 1);
+}
+
+struct matrix_function
+ {
+ const char *name;
+ enum matrix_op op;
+ size_t min_args, max_args;
+ };
+
+static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
+
+static bool
+word_matches (const char **test, const char **name)
+{
+ size_t test_len = strcspn (*test, ".");
+ size_t name_len = strcspn (*name, ".");
+ if (test_len == name_len)
+ {
+ if (buf_compare_case (*test, *name, test_len))
+ return false;
+ }
+ else if (test_len < 3 || test_len > name_len)
+ return false;
+ else
+ {
+ if (buf_compare_case (*test, *name, test_len))
+ return false;
+ }
+
+ *test += test_len;
+ *name += name_len;
+ if (**test != **name)
+ return false;
+
+ if (**test == '.')
+ {
+ (*test)++;
+ (*name)++;
+ }
+ return true;
+}
+
+/* Returns 0 if TOKEN and FUNC do not match,
+ 1 if TOKEN is an acceptable abbreviation for FUNC,
+ 2 if TOKEN equals FUNC. */
+static int
+compare_function_names (const char *token_, const char *func_)
+{
+ const char *token = token_;
+ const char *func = func_;
+ while (*token || *func)
+ if (!word_matches (&token, &func))
+ return 0;
+ return !c_strcasecmp (token_, func_) ? 2 : 1;
+}
+
+static const struct matrix_function *
+matrix_parse_function_name (const char *token)
+{
+ static const struct matrix_function functions[] =
+ {
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+ { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
+ MATRIX_FUNCTIONS
+#undef F
+ };
+ enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
+
+ for (size_t i = 0; i < N_FUNCTIONS; i++)
+ {
+ if (compare_function_names (token, functions[i].name) > 0)
+ return &functions[i];
+ }
+ return NULL;
+}
+
+static bool
+matrix_parse_function (struct matrix_state *s, const char *token,
+ struct matrix_expr **exprp)
+{
+ *exprp = NULL;
+ if (lex_next_token (s->lexer, 1) != T_LPAREN)
+ return false;
+
+ int start_ofs = lex_ofs (s->lexer);
+ if (lex_match_id (s->lexer, "EOF"))
+ {
+ lex_get (s->lexer);
+ struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
+ if (!fh)
+ return true;
+
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ {
+ fh_unref (fh);
+ return true;
+ }
+
+ struct read_file *rf = read_file_create (s, fh);
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
+ matrix_expr_add_location (s, start_ofs, e);
+ *exprp = e;
+ return true;
+ }
+
+ const struct matrix_function *f = matrix_parse_function_name (token);
+ if (!f)
+ return false;
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = f->op };
+
+ lex_get_n (s->lexer, 2);
+ if (lex_token (s->lexer) != T_RPAREN)
+ {
+ size_t allocated_subs = 0;
+ do
+ {
+ struct matrix_expr *sub = matrix_expr_parse (s);
+ if (!sub)
+ goto error;
+
+ if (e->n_subs >= allocated_subs)
+ e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
+ e->subs[e->n_subs++] = sub;
+ }
+ while (lex_match (s->lexer, T_COMMA));
+ }
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ if (e->n_subs < f->min_args || e->n_subs > f->max_args)
+ {
+ if (f->min_args == f->max_args)
+ msg_at (SE, e->location,
+ ngettext ("Matrix function %s requires %zu argument.",
+ "Matrix function %s requires %zu arguments.",
+ f->min_args),
+ f->name, f->min_args);
+ else if (f->min_args == 1 && f->max_args == 2)
+ msg_at (SE, e->location,
+ ngettext ("Matrix function %s requires 1 or 2 arguments, "
+ "but %zu was provided.",
+ "Matrix function %s requires 1 or 2 arguments, "
+ "but %zu were provided.",
+ e->n_subs),
+ f->name, e->n_subs);
+ else if (f->min_args == 1 && f->max_args == INT_MAX)
+ msg_at (SE, e->location,
+ _("Matrix function %s requires at least one argument."),
+ f->name);
+ else
+ NOT_REACHED ();
+
+ goto error;
+ }
+
+ matrix_expr_add_location (s, start_ofs, e);
+
+ *exprp = e;
+ return true;
+
+error:
+ matrix_expr_destroy (e);
+ return true;
+}
+
+static struct matrix_expr *
+matrix_parse_primary__ (struct matrix_state *s)
+{
+ if (lex_is_number (s->lexer))
+ {
+ double number = lex_number (s->lexer);
+ lex_get (s->lexer);
+
+ return matrix_expr_create_number (number);
+ }
+ else if (lex_is_string (s->lexer))
+ {
+ char string[sizeof (double)];
+ buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
+ lex_get (s->lexer);
+
+ double number;
+ memcpy (&number, string, sizeof number);
+
+ return matrix_expr_create_number (number);
+ }
+ else if (lex_match (s->lexer, T_LPAREN))
+ {
+ struct matrix_expr *e = matrix_expr_parse (s);
+ if (!e || !lex_force_match (s->lexer, T_RPAREN))
+ {
+ matrix_expr_destroy (e);
+ return NULL;
+ }
+ return e;
+ }
+ else if (lex_match (s->lexer, T_LCURLY))
+ {
+ struct matrix_expr *e = matrix_parse_curly_semi (s);
+ if (!e || !lex_force_match (s->lexer, T_RCURLY))
+ {
+ matrix_expr_destroy (e);
+ return NULL;
+ }
+ return e;
+ }
+ else if (lex_token (s->lexer) == T_ID)
+ {
+ struct matrix_expr *retval;
+ if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
+ return retval;
+
+ struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (!var)
+ {
+ lex_error (s->lexer, _("Unknown variable %s."),
+ lex_tokcstr (s->lexer));
+ return NULL;
+ }
+ lex_get (s->lexer);
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
+ return e;
+ }
+ else if (lex_token (s->lexer) == T_ALL)
+ {
+ struct matrix_expr *retval;
+ if (matrix_parse_function (s, "ALL", &retval))
+ return retval;
+ }
+
+ lex_error (s->lexer, _("Syntax error expecting matrix expression."));
+ return NULL;
+}
+
+static struct matrix_expr *
+matrix_parse_primary (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+ struct matrix_expr *e = matrix_parse_primary__ (s);
+ matrix_expr_add_location (s, start_ofs, e);
+ return e;
+}
+
+static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
+
+static bool
+matrix_parse_index_expr (struct matrix_state *s,
+ struct matrix_expr **indexp,
+ struct msg_location **locationp)
+{
+ if (lex_match (s->lexer, T_COLON))
+ {
+ if (locationp)
+ *locationp = lex_get_location (s->lexer, -1, -1);
+ *indexp = NULL;
+ return true;
+ }
+ else
+ {
+ *indexp = matrix_expr_parse (s);
+ if (locationp && *indexp)
+ *locationp = msg_location_dup (matrix_expr_location (*indexp));
+ return *indexp != NULL;
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_postfix (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_primary (s);
+ if (!lhs || !lex_match (s->lexer, T_LPAREN))
+ return lhs;
+
+ struct matrix_expr *i0;
+ if (!matrix_parse_index_expr (s, &i0, NULL))
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ if (lex_match (s->lexer, T_RPAREN))
+ return (i0
+ ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
+ : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
+ else if (lex_match (s->lexer, T_COMMA))
+ {
+ struct matrix_expr *i1;
+ if (!matrix_parse_index_expr (s, &i1, NULL)
+ || !lex_force_match (s->lexer, T_RPAREN))
+ {
+ matrix_expr_destroy (lhs);
+ matrix_expr_destroy (i0);
+ matrix_expr_destroy (i1);
+ return NULL;
+ }
+ return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
+ : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
+ : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
+ : lhs);
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "`)'", "`,'");
+ return NULL;
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_unary (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+
+ struct matrix_expr *e;
+ if (lex_match (s->lexer, T_DASH))
+ {
+ struct matrix_expr *sub = matrix_parse_unary (s);
+ if (!sub)
+ return NULL;
+ e = matrix_expr_create_1 (MOP_NEGATE, sub);
+ }
+ else if (lex_match (s->lexer, T_PLUS))
+ {
+ e = matrix_parse_unary (s);
+ if (!e)
+ return NULL;
+ }
+ else
+ return matrix_parse_postfix (s);
+
+ matrix_expr_add_location (s, start_ofs, e);
+ e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
+ return e;
+}
+
+static struct matrix_expr *
+matrix_parse_seq (struct matrix_state *s)
+{
+ struct matrix_expr *start = matrix_parse_unary (s);
+ if (!start || !lex_match (s->lexer, T_COLON))
+ return start;
+
+ struct matrix_expr *end = matrix_parse_unary (s);
+ if (!end)
+ {
+ matrix_expr_destroy (start);
+ return NULL;
+ }
+
+ if (lex_match (s->lexer, T_COLON))
+ {
+ struct matrix_expr *increment = matrix_parse_unary (s);
+ if (!increment)
+ {
+ matrix_expr_destroy (start);
+ matrix_expr_destroy (end);
+ return NULL;
+ }
+ return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
+ }
+ else
+ return matrix_expr_create_2 (MOP_SEQ, start, end);
+}
+
+static struct matrix_expr *
+matrix_parse_exp (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax syntax[] = {
+ { .token = T_EXP, .op = MOP_EXP_MAT },
+ { .phrase = "&**", .op = MOP_EXP_ELEMS },
+ };
+ size_t n_syntax = sizeof syntax / sizeof *syntax;
+
+ return matrix_parse_binary_operator (s, matrix_parse_seq, syntax, n_syntax);
+}
+
+static struct matrix_expr *
+matrix_parse_mul_div (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax syntax[] = {
+ { .token = T_ASTERISK, .op = MOP_MUL_MAT },
+ { .token = T_SLASH, .op = MOP_DIV_ELEMS },
+ { .phrase = "&*", .op = MOP_MUL_ELEMS },
+ { .phrase = "&/", .op = MOP_DIV_ELEMS },
+ };
+ size_t n_syntax = sizeof syntax / sizeof *syntax;
+
+ return matrix_parse_binary_operator (s, matrix_parse_exp, syntax, n_syntax);
+}
+
+static struct matrix_expr *
+matrix_parse_add_sub (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_mul_div (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (lex_match (s->lexer, T_PLUS))
+ op = MOP_ADD_ELEMS;
+ else if (lex_match (s->lexer, T_DASH))
+ op = MOP_SUB_ELEMS;
+ else if (lex_token (s->lexer) == T_NEG_NUM)
+ op = MOP_ADD_ELEMS;
+ else
+ return lhs;
+
+ struct matrix_expr *rhs = matrix_parse_mul_div (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_relational (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax syntax[] = {
+ { .token = T_GT, .op = MOP_GT },
+ { .token = T_GE, .op = MOP_GE },
+ { .token = T_LT, .op = MOP_LT },
+ { .token = T_LE, .op = MOP_LE },
+ { .token = T_EQUALS, .op = MOP_EQ },
+ { .token = T_EQ, .op = MOP_EQ },
+ { .token = T_NE, .op = MOP_NE },
+ };
+ size_t n_syntax = sizeof syntax / sizeof *syntax;
+
+ return matrix_parse_binary_operator (s, matrix_parse_add_sub,
+ syntax, n_syntax);
+}
+
+static struct matrix_expr *
+matrix_parse_not (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+ if (lex_match (s->lexer, T_NOT))
+ {
+ struct matrix_expr *sub = matrix_parse_not (s);
+ if (!sub)
+ return NULL;
+
+ struct matrix_expr *e = matrix_expr_create_1 (MOP_NOT, sub);
+ matrix_expr_add_location (s, start_ofs, e);
+ e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
+ return e;
+ }
+ else
+ return matrix_parse_relational (s);
+}
+
+static struct matrix_expr *
+matrix_parse_and (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax op = {
+ .token = T_AND, .op = MOP_AND
+ };
+
+ return matrix_parse_binary_operator (s, matrix_parse_not, &op, 1);
+}
+
+static struct matrix_expr *
+matrix_expr_parse__ (struct matrix_state *s)
+{
+ static const struct matrix_operator_syntax syntax[] = {
+ { .token = T_OR, .op = MOP_OR },
+ { .id = "XOR", .op = MOP_XOR },
+ };
+ size_t n_syntax = sizeof syntax / sizeof *syntax;
+
+ return matrix_parse_binary_operator (s, matrix_parse_and, syntax, n_syntax);
+}
+
+static struct matrix_expr *
+matrix_expr_parse (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+ struct matrix_expr *e = matrix_expr_parse__ (s);
+ matrix_expr_add_location (s, start_ofs, e);
+ return e;
+}
+\f
+/* Matrix expression evaluation. */
+
+/* Iterates over all the elements in matrix M, setting Y and X to the row and
+ column indexes, respectively, and pointing D to the entry at each
+ position. */
+#define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
+ for (size_t Y = 0; Y < (M)->size1; Y++) \
+ for (size_t X = 0; X < (M)->size2; X++) \
+ for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
+
+static bool
+is_vector (const gsl_matrix *m)
+{
+ return m->size1 <= 1 || m->size2 <= 1;
+}
+
+static gsl_vector
+to_vector (gsl_matrix *m)
+{
+ return (m->size1 == 1
+ ? gsl_matrix_row (m, 0).vector
+ : gsl_matrix_column (m, 0).vector);
+}
+
+static double
+matrix_eval_ABS (double d)
+{
+ return fabs (d);
+}
+
+static double
+matrix_eval_ALL (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ if (*d == 0.0)
+ return 0.0;
+ return 1.0;
+}
+
+static double
+matrix_eval_ANY (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ if (*d != 0.0)
+ return 1.0;
+ return 0.0;
+}
+
+static double
+matrix_eval_ARSIN (double d)
+{
+ return asin (d);
+}
+
+static double
+matrix_eval_ARTAN (double d)
+{
+ return atan (d);
+}
+
+static gsl_matrix *
+matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
+{
+ size_t r = 0;
+ size_t c = 0;
+ for (size_t i = 0; i < n; i++)
+ {
+ r += m[i]->size1;
+ c += m[i]->size2;
+ }
+ gsl_matrix *block = gsl_matrix_calloc (r, c);
+ r = c = 0;
+ for (size_t i = 0; i < n; i++)
+ {
+ for (size_t y = 0; y < m[i]->size1; y++)
+ for (size_t x = 0; x < m[i]->size2; x++)
+ gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
+ r += m[i]->size1;
+ c += m[i]->size2;
+ }
+ return block;
+}
+
+static gsl_matrix *
+matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e)
+{
+ if (!gsl_linalg_cholesky_decomp1 (m))
+ {
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = y + 1; x < m->size2; x++)
+ gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
+
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = 0; x < y; x++)
+ gsl_matrix_set (m, y, x, 0);
+ return m;
+ }
+ else
+ {
+ msg_at (SE, e->subs[0]->location,
+ _("Input to CHOL function is not positive-definite."));
+ return NULL;
+ }
+}
+
+static gsl_matrix *
+matrix_eval_col_extremum (gsl_matrix *m, bool min)
+{
+ if (m->size1 <= 1)
+ return m;
+ else if (!m->size2)
+ return gsl_matrix_alloc (1, 0);
+
+ gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double ext = gsl_matrix_get (m, 0, x);
+ for (size_t y = 1; y < m->size1; y++)
+ {
+ double value = gsl_matrix_get (m, y, x);
+ if (min ? value < ext : value > ext)
+ ext = value;
+ }
+ gsl_matrix_set (cext, 0, x, ext);
+ }
+ return cext;
+}
+
+static gsl_matrix *
+matrix_eval_CMAX (gsl_matrix *m)
+{
+ return matrix_eval_col_extremum (m, false);
+}
+
+static gsl_matrix *
+matrix_eval_CMIN (gsl_matrix *m)
+{
+ return matrix_eval_col_extremum (m, true);
+}
+
+static double
+matrix_eval_COS (double d)
+{
+ return cos (d);
+}
+
+static gsl_matrix *
+matrix_eval_col_sum (gsl_matrix *m, bool square)
+{
+ if (m->size1 == 0)
+ return m;
+ else if (!m->size2)
+ return gsl_matrix_alloc (1, 0);
+
+ gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double sum = 0;
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ double d = gsl_matrix_get (m, y, x);
+ sum += square ? pow2 (d) : d;
+ }
+ gsl_matrix_set (result, 0, x, sum);
+ }
+ return result;
+}
+
+static gsl_matrix *
+matrix_eval_CSSQ (gsl_matrix *m)
+{
+ return matrix_eval_col_sum (m, true);
+}
+
+static gsl_matrix *
+matrix_eval_CSUM (gsl_matrix *m)
+{
+ return matrix_eval_col_sum (m, false);
+}
+
+static int
+compare_double_3way (const void *a_, const void *b_)
+{
+ const double *a = a_;
+ const double *b = b_;
+ return *a < *b ? -1 : *a > *b;
+}
+
+static gsl_matrix *
+matrix_eval_DESIGN (gsl_matrix *m, const struct matrix_expr *e)
+{
+ double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
+ gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
+ gsl_matrix_transpose_memcpy (&m2, m);
+
+ for (size_t y = 0; y < m2.size1; y++)
+ qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
+
+ size_t *n = xcalloc (m2.size1, sizeof *n);
+ size_t n_total = 0;
+ for (size_t i = 0; i < m2.size1; i++)
+ {
+ double *row = tmp + m2.size2 * i;
+ for (size_t j = 0; j < m2.size2; )
+ {
+ size_t k;
+ for (k = j + 1; k < m2.size2; k++)
+ if (row[j] != row[k])
+ break;
+ row[n[i]++] = row[j];
+ j = k;
+ }
+
+ if (n[i] <= 1)
+ msg_at (MW, e->subs[0]->location,
+ _("Column %zu in DESIGN argument has constant value."), i + 1);
+ else
+ n_total += n[i];
+ }
+
+ gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
+ size_t x = 0;
+ for (size_t i = 0; i < m->size2; i++)
+ {
+ if (n[i] <= 1)
+ continue;
+
+ const double *unique = tmp + m2.size2 * i;
+ for (size_t j = 0; j < n[i]; j++, x++)
+ {
+ double value = unique[j];
+ for (size_t y = 0; y < m->size1; y++)
+ gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
+ }
+ }
+
+ free (n);
+ free (tmp);
+
+ return result;
+}
+
+static double
+matrix_eval_DET (gsl_matrix *m)
+{
+ gsl_permutation *p = gsl_permutation_alloc (m->size1);
+ int signum;
+ gsl_linalg_LU_decomp (m, p, &signum);
+ gsl_permutation_free (p);
+ return gsl_linalg_LU_det (m, signum);
+}
+
+static gsl_matrix *
+matrix_eval_DIAG (gsl_matrix *m)
+{
+ gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
+ for (size_t i = 0; i < diag->size1; i++)
+ gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
+ return diag;
+}
+
+static bool
+is_symmetric (const gsl_matrix *m)
+{
+ if (m->size1 != m->size2)
+ return false;
+
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = 0; x < y; x++)
+ if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
+ return false;
+
+ return true;
+}
+
+static int
+compare_double_desc (const void *a_, const void *b_)
+{
+ const double *a = a_;
+ const double *b = b_;
+ return *a > *b ? -1 : *a < *b;
+}
+
+static gsl_matrix *
+matrix_eval_EVAL (gsl_matrix *m, const struct matrix_expr *e)
+{
+ if (!is_symmetric (m))
+ {
+ msg_at (SE, e->subs[0]->location,
+ _("Argument of EVAL must be symmetric."));
+ return NULL;
+ }
+
+ gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
+ gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
+ gsl_vector v_eval = to_vector (eval);
+ gsl_eigen_symm (m, &v_eval, w);
+ gsl_eigen_symm_free (w);
+
+ assert (v_eval.stride == 1);
+ qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
+
+ return eval;
+}
+
+static double
+matrix_eval_EXP (double d)
+{
+ return exp (d);
+}
+
+/* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
+ marked as:
+
+ Charl Linssen <charl@itfromb.it>
+ Feb 2016
+ PUBLIC DOMAIN */
+static gsl_matrix *
+matrix_eval_GINV (gsl_matrix *A)
+{
+ size_t n = A->size1;
+ size_t m = A->size2;
+ bool swap = m > n;
+ gsl_matrix *tmp_mat = NULL;
+ if (swap)
+ {
+ /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
+ tmp_mat = gsl_matrix_alloc (m, n);
+ gsl_matrix_transpose_memcpy (tmp_mat, A);
+ A = tmp_mat;
+ size_t i = m;
+ m = n;
+ n = i;
+ }
+
+ /* Do SVD. */
+ gsl_matrix *V = gsl_matrix_alloc (m, m);
+ gsl_vector *u = gsl_vector_alloc (m);
+
+ gsl_vector *tmp_vec = gsl_vector_alloc (m);
+ gsl_linalg_SV_decomp (A, V, u, tmp_vec);
+ gsl_vector_free (tmp_vec);
+
+ /* Compute Σ⁻¹. */
+ gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
+ gsl_matrix_set_zero (Sigma_pinv);
+ double cutoff = 1e-15 * gsl_vector_max (u);
+
+ for (size_t i = 0; i < m; ++i)
+ {
+ double x = gsl_vector_get (u, i);
+ gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
+ }
+
+ /* libgsl SVD yields "thin" SVD. Pad to full matrix by adding zeros. */
+ gsl_matrix *U = gsl_matrix_calloc (n, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
+
+ /* Two dot products to obtain pseudoinverse. */
+ gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
+ gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
+
+ gsl_matrix *A_pinv;
+ if (swap)
+ {
+ A_pinv = gsl_matrix_alloc (n, m);
+ gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
+ }
+ else
+ {
+ A_pinv = gsl_matrix_alloc (m, n);
+ gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
+ }
+
+ gsl_matrix_free (tmp_mat);
+ gsl_matrix_free (tmp_mat2);
+ gsl_matrix_free (U);
+ gsl_matrix_free (Sigma_pinv);
+ gsl_vector_free (u);
+ gsl_matrix_free (V);
+
+ return A_pinv;
+}
+
+struct grade
+ {
+ size_t y, x;
+ double value;
+ };
+
+static int
+grade_compare_3way (const void *a_, const void *b_)
+{
+ const struct grade *a = a_;
+ const struct grade *b = b_;
+
+ return (a->value < b->value ? -1
+ : a->value > b->value ? 1
+ : a->y < b->y ? -1
+ : a->y > b->y ? 1
+ : a->x < b->x ? -1
+ : a->x > b->x);
+}
+
+static gsl_matrix *
+matrix_eval_GRADE (gsl_matrix *m)
+{
+ size_t n = m->size1 * m->size2;
+ struct grade *grades = xmalloc (n * sizeof *grades);
+
+ size_t i = 0;
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
+ qsort (grades, n, sizeof *grades, grade_compare_3way);
+
+ for (size_t i = 0; i < n; i++)
+ gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
+
+ free (grades);
+
+ return m;
+}
+
+static double
+dot (gsl_vector *a, gsl_vector *b)
+{
+ double result = 0.0;
+ for (size_t i = 0; i < a->size; i++)
+ result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
+ return result;
+}
+
+static double
+norm2 (gsl_vector *v)
+{
+ double result = 0.0;
+ for (size_t i = 0; i < v->size; i++)
+ result += pow2 (gsl_vector_get (v, i));
+ return result;
+}
+
+static double
+norm (gsl_vector *v)
+{
+ return sqrt (norm2 (v));
+}
+
+static gsl_matrix *
+matrix_eval_GSCH (gsl_matrix *v, const struct matrix_expr *e)
+{
+ if (v->size2 < v->size1)
+ {
+ msg_at (SE, e->subs[0]->location,
+ _("GSCH requires its argument to have at least as many columns "
+ "as rows, but it has dimensions %zu×%zu."),
+ v->size1, v->size2);
+ return NULL;
+ }
+ if (!v->size1 || !v->size2)
+ return v;
+
+ gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
+ size_t ux = 0;
+ for (size_t vx = 0; vx < v->size2; vx++)
+ {
+ gsl_vector u_i = gsl_matrix_column (u, ux).vector;
+ gsl_vector v_i = gsl_matrix_column (v, vx).vector;
+
+ gsl_vector_memcpy (&u_i, &v_i);
+ for (size_t j = 0; j < ux; j++)
+ {
+ gsl_vector u_j = gsl_matrix_column (u, j).vector;
+ double scale = dot (&u_j, &u_i) / norm2 (&u_j);
+ for (size_t k = 0; k < u_i.size; k++)
+ gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
+ - scale * gsl_vector_get (&u_j, k)));
+ }
+
+ double len = norm (&u_i);
+ if (len > 1e-15)
+ {
+ gsl_vector_scale (&u_i, 1.0 / len);
+ if (++ux >= v->size1)
+ break;
+ }
+ }
+
+ if (ux < v->size1)
+ {
+ msg_at (SE, e->subs[0]->location,
+ _("%zu×%zu argument to GSCH contains only "
+ "%zu linearly independent columns."),
+ v->size1, v->size2, ux);
+ gsl_matrix_free (u);
+ return NULL;
+ }
+
+ u->size2 = v->size1;
+ return u;
+}
+
+static gsl_matrix *
+matrix_eval_IDENT (double s1, double s2)
+{
+ gsl_matrix *m = gsl_matrix_alloc (s1, s2);
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = x == y;
+ return m;
+}
+
+/* Inverts X, storing the inverse into INVERSE. As a side effect, replaces X
+ by its LU decomposition. */
+static void
+invert_matrix (gsl_matrix *x, gsl_matrix *inverse)
+{
+ gsl_permutation *p = gsl_permutation_alloc (x->size1);
+ int signum;
+ gsl_linalg_LU_decomp (x, p, &signum);
+ gsl_linalg_LU_invert (x, p, inverse);
+ gsl_permutation_free (p);
+}
+
+static gsl_matrix *
+matrix_eval_INV (gsl_matrix *src)
+{
+ gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
+ invert_matrix (src, dst);
+ return dst;
+}
+
+static gsl_matrix *
+matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
+{
+ gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
+ a->size2 * b->size2);
+ size_t y = 0;
+ for (size_t ar = 0; ar < a->size1; ar++)
+ for (size_t br = 0; br < b->size1; br++, y++)
+ {
+ size_t x = 0;
+ for (size_t ac = 0; ac < a->size2; ac++)
+ for (size_t bc = 0; bc < b->size2; bc++, x++)
+ {
+ double av = gsl_matrix_get (a, ar, ac);
+ double bv = gsl_matrix_get (b, br, bc);
+ gsl_matrix_set (k, y, x, av * bv);
+ }
+ }
+ return k;
+}
+
+static double
+matrix_eval_LG10 (double d)
+{
+ return log10 (d);
+}
+
+static double
+matrix_eval_LN (double d)
+{
+ return log (d);
+}
+
+static void
+matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
+{
+ /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
+ size_t y = 0;
+ size_t x = n / 2;
+ for (size_t i = 1; i <= n * n; i++)
+ {
+ gsl_matrix_set (m, y, x, i);
+
+ size_t y1 = !y ? n - 1 : y - 1;
+ size_t x1 = x + 1 >= n ? 0 : x + 1;
+ if (gsl_matrix_get (m, y1, x1) == 0)
+ {
+ y = y1;
+ x = x1;
+ }
+ else
+ y = y + 1 >= n ? 0 : y + 1;
+ }
+}
+
+static void
+magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
+{
+ double a = gsl_matrix_get (m, y1, x1);
+ double b = gsl_matrix_get (m, y2, x2);
+ gsl_matrix_set (m, y1, x1, b);
+ gsl_matrix_set (m, y2, x2, a);
+}
+
+static void
+matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
+{
+ size_t x, y;
+
+ /* A. Umar, "On the Construction of Even Order Magic Squares",
+ https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
+ x = y = 0;
+ for (size_t i = 1; i <= n * n / 2; i++)
+ {
+ gsl_matrix_set (m, y, x, i);
+ if (++y >= n)
+ {
+ y = 0;
+ x++;
+ }
+ }
+
+ x = n - 1;
+ y = 0;
+ for (size_t i = n * n; i > n * n / 2; i--)
+ {
+ gsl_matrix_set (m, y, x, i);
+ if (++y >= n)
+ {
+ y = 0;
+ x--;
+ }
+ }
+
+ for (size_t y = 0; y < n; y++)
+ for (size_t x = 0; x < n / 2; x++)
+ {
+ unsigned int d = gsl_matrix_get (m, y, x);
+ if (d % 2 != (y < n / 2))
+ magic_exchange (m, y, x, y, n - x - 1);
+ }
+
+ size_t y1 = n / 2;
+ size_t y2 = n - 1;
+ size_t x1 = n / 2 - 1;
+ size_t x2 = n / 2;
+ magic_exchange (m, y1, x1, y2, x1);
+ magic_exchange (m, y1, x2, y2, x2);
+}
+
+static void
+matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
+{
+ /* A. Umar, "On the Construction of Even Order Magic Squares",
+ https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
+ size_t x, y;
+
+ x = y = 0;
+ for (size_t i = 1; ; i++)
+ {
+ gsl_matrix_set (m, y, x, i);
+ if (++y == n / 2 - 1)
+ y += 2;
+ else if (y >= n)
+ {
+ y = 0;
+ if (++x >= n / 2)
+ break;
+ }
+ }
+
+ x = n - 1;
+ y = 0;
+ for (size_t i = n * n; ; i--)
+ {
+ gsl_matrix_set (m, y, x, i);
+ if (++y == n / 2 - 1)
+ y += 2;
+ else if (y >= n)
+ {
+ y = 0;
+ if (--x < n / 2)
+ break;
+ }
+ }
+ for (size_t y = 0; y < n; y++)
+ if (y != n / 2 - 1 && y != n / 2)
+ for (size_t x = 0; x < n / 2; x++)
+ {
+ unsigned int d = gsl_matrix_get (m, y, x);
+ if (d % 2 != (y < n / 2))
+ magic_exchange (m, y, x, y, n - x - 1);
+ }
+
+ size_t a0 = (n * n - 2 * n) / 2 + 1;
+ for (size_t i = 0; i < n / 2; i++)
+ {
+ size_t a = a0 + i;
+ gsl_matrix_set (m, n / 2, i, a);
+ gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
+ }
+ for (size_t i = 0; i < n / 2; i++)
+ {
+ size_t a = a0 + i + n / 2;
+ gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
+ gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
+ }
+ for (size_t x = 1; x < n / 2; x += 2)
+ magic_exchange (m, n / 2, x, n / 2 - 1, x);
+ for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
+ magic_exchange (m, n / 2, x, n / 2 - 1, x);
+ size_t x1 = n / 2 - 2;
+ size_t x2 = n / 2 + 1;
+ size_t y1 = n / 2 - 2;
+ size_t y2 = n / 2 + 1;
+ magic_exchange (m, y1, x1, y2, x1);
+ magic_exchange (m, y1, x2, y2, x2);
+}
+
+static gsl_matrix *
+matrix_eval_MAGIC (double n_)
+{
+ size_t n = n_;
+
+ gsl_matrix *m = gsl_matrix_calloc (n, n);
+ if (n % 2)
+ matrix_eval_MAGIC_odd (m, n);
+ else if (n % 4)
+ matrix_eval_MAGIC_singly_even (m, n);
+ else
+ matrix_eval_MAGIC_doubly_even (m, n);
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_MAKE (double r, double c, double value)
+{
+ gsl_matrix *m = gsl_matrix_alloc (r, c);
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = value;
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_MDIAG (gsl_vector *v)
+{
+ gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
+ gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
+ gsl_vector_memcpy (&diagonal, v);
+ return m;
+}
+
+static double
+matrix_eval_MMAX (gsl_matrix *m)
+{
+ return gsl_matrix_max (m);
+}
+
+static double
+matrix_eval_MMIN (gsl_matrix *m)
+{
+ return gsl_matrix_min (m);
+}
+
+static gsl_matrix *
+matrix_eval_MOD (gsl_matrix *m, double divisor)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = fmod (*d, divisor);
+ return m;
+}
+
+static double
+matrix_eval_MSSQ (gsl_matrix *m)
+{
+ double mssq = 0.0;
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ mssq += *d * *d;
+ return mssq;
+}
+
+static double
+matrix_eval_MSUM (gsl_matrix *m)
+{
+ double msum = 0.0;
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ msum += *d;
+ return msum;
+}
+
+static double
+matrix_eval_NCOL (gsl_matrix *m)
+{
+ return m->size2;
+}
+
+static double
+matrix_eval_NROW (gsl_matrix *m)
+{
+ return m->size1;
+}
+
+static double
+matrix_eval_RANK (gsl_matrix *m)
+{
+ gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
+ gsl_linalg_QR_decomp (m, tau);
+ gsl_vector_free (tau);
+
+ return gsl_linalg_QRPT_rank (m, -1);
+}
+
+static gsl_matrix *
+matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_,
+ const struct matrix_expr *e)
+{
+ bool r_ok = r_ >= 0 && r_ < SIZE_MAX;
+ bool c_ok = c_ >= 0 && c_ < SIZE_MAX;
+ if (!r_ok || !c_ok)
+ {
+ msg_at (SE,
+ !r_ok ? e->subs[1]->location : e->subs[2]->location,
+ _("Arguments 2 and 3 to RESHAPE must be integers."));
+ return NULL;
+ }
+ size_t r = r_;
+ size_t c = c_;
+ if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
+ {
+ struct msg_location *loc = msg_location_dup (e->subs[1]->location);
+ loc->end = e->subs[2]->location->end;
+ msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) "
+ "differs from product of matrix dimensions "
+ "(%zu×%zu = %zu)."),
+ r, c, r * c,
+ m->size1, m->size2, m->size1 * m->size2);
+ msg_location_destroy (loc);
+ return NULL;
+ }
+
+ gsl_matrix *dst = gsl_matrix_alloc (r, c);
+ size_t y1 = 0;
+ size_t x1 = 0;
+ MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
+ {
+ gsl_matrix_set (dst, y1, x1, *d);
+ if (++x1 >= c)
+ {
+ x1 = 0;
+ y1++;
+ }
+ }
+ return dst;
+}
+
+static gsl_matrix *
+matrix_eval_row_extremum (gsl_matrix *m, bool min)
+{
+ if (m->size2 <= 1)
+ return m;
+ else if (!m->size1)
+ return gsl_matrix_alloc (0, 1);
+
+ gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ double ext = gsl_matrix_get (m, y, 0);
+ for (size_t x = 1; x < m->size2; x++)
+ {
+ double value = gsl_matrix_get (m, y, x);
+ if (min ? value < ext : value > ext)
+ ext = value;
+ }
+ gsl_matrix_set (rext, y, 0, ext);
+ }
+ return rext;
+}
+
+static gsl_matrix *
+matrix_eval_RMAX (gsl_matrix *m)
+{
+ return matrix_eval_row_extremum (m, false);
+}
+
+static gsl_matrix *
+matrix_eval_RMIN (gsl_matrix *m)
+{
+ return matrix_eval_row_extremum (m, true);
+}
+
+static double
+matrix_eval_RND (double d)
+{
+ return rint (d);
+}
+
+struct rank
+ {
+ size_t y, x;
+ double value;
+ };
+
+static int
+rank_compare_3way (const void *a_, const void *b_)
+{
+ const struct rank *a = a_;
+ const struct rank *b = b_;
+
+ return a->value < b->value ? -1 : a->value > b->value;
+}
+
+static gsl_matrix *
+matrix_eval_RNKORDER (gsl_matrix *m)
+{
+ size_t n = m->size1 * m->size2;
+ struct rank *ranks = xmalloc (n * sizeof *ranks);
+ size_t i = 0;
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
+ qsort (ranks, n, sizeof *ranks, rank_compare_3way);
+
+ for (size_t i = 0; i < n; )
+ {
+ size_t j;
+ for (j = i + 1; j < n; j++)
+ if (ranks[i].value != ranks[j].value)
+ break;
+
+ double rank = (i + j + 1.0) / 2.0;
+ for (size_t k = i; k < j; k++)
+ gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
+
+ i = j;
+ }
+
+ free (ranks);
+
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_row_sum (gsl_matrix *m, bool square)
+{
+ if (m->size1 == 0)
+ return m;
+ else if (!m->size1)
+ return gsl_matrix_alloc (0, 1);
+
+ gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ double sum = 0;
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double d = gsl_matrix_get (m, y, x);
+ sum += square ? pow2 (d) : d;
+ }
+ gsl_matrix_set (result, y, 0, sum);
+ }
+ return result;
+}
+
+static gsl_matrix *
+matrix_eval_RSSQ (gsl_matrix *m)
+{
+ return matrix_eval_row_sum (m, true);
+}
+
+static gsl_matrix *
+matrix_eval_RSUM (gsl_matrix *m)
+{
+ return matrix_eval_row_sum (m, false);
+}
+
+static double
+matrix_eval_SIN (double d)
+{
+ return sin (d);
+}
+
+static gsl_matrix *
+matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e)
+{
+ if (m1->size1 != m2->size1)
+ {
+ struct msg_location *loc = msg_location_dup (e->subs[0]->location);
+ loc->end = e->subs[1]->location->end;
+
+ msg_at (SE, e->location,
+ _("SOLVE arguments must have the same number of rows."));
+ msg_at (SN, e->subs[0]->location,
+ _("Argument 1 has dimensions %zu×%zu."), m1->size1, m1->size2);
+ msg_at (SN, e->subs[1]->location,
+ _("Argument 2 has dimensions %zu×%zu."), m2->size1, m2->size2);
+
+ msg_location_destroy (loc);
+ return NULL;
+ }
+
+ gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
+ gsl_permutation *p = gsl_permutation_alloc (m1->size1);
+ int signum;
+ gsl_linalg_LU_decomp (m1, p, &signum);
+ for (size_t i = 0; i < m2->size2; i++)
+ {
+ gsl_vector bi = gsl_matrix_column (m2, i).vector;
+ gsl_vector xi = gsl_matrix_column (x, i).vector;
+ gsl_linalg_LU_solve (m1, p, &bi, &xi);
+ }
+ gsl_permutation_free (p);
+ return x;
+}
+
+static double
+matrix_eval_SQRT (double d)
+{
+ return sqrt (d);
+}
+
+static gsl_matrix *
+matrix_eval_SSCP (gsl_matrix *m)
+{
+ gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
+ gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
+ return sscp;
+}
+
+static gsl_matrix *
+matrix_eval_SVAL (gsl_matrix *m)
+{
+ gsl_matrix *tmp_mat = NULL;
+ if (m->size2 > m->size1)
+ {
+ tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
+ gsl_matrix_transpose_memcpy (tmp_mat, m);
+ m = tmp_mat;
+ }
+
+ /* Do SVD. */
+ gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
+ gsl_vector *S = gsl_vector_alloc (m->size2);
+ gsl_vector *work = gsl_vector_alloc (m->size2);
+ gsl_linalg_SV_decomp (m, V, S, work);
+
+ gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
+ for (size_t i = 0; i < m->size2; i++)
+ gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
+
+ gsl_matrix_free (V);
+ gsl_vector_free (S);
+ gsl_vector_free (work);
+ gsl_matrix_free (tmp_mat);
+
+ return vals;
+}
+
+static gsl_matrix *
+matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e)
+{
+ if (d < 1 || d > SIZE_MAX)
+ {
+ msg_at (SE, e->subs[1]->location,
+ _("Scalar argument to SWEEP must be integer."));
+ return NULL;
+ }
+ size_t k = d - 1;
+ if (k >= MIN (m->size1, m->size2))
+ {
+ msg_at (SE, e->subs[1]->location,
+ _("Scalar argument to SWEEP must be integer less than or "
+ "equal to the smaller of the matrix argument's rows and "
+ "columns."));
+ return NULL;
+ }
+
+ double m_kk = gsl_matrix_get (m, k, k);
+ if (fabs (m_kk) > 1e-19)
+ {
+ gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
+ MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
+ {
+ double m_ij = gsl_matrix_get (m, i, j);
+ double m_ik = gsl_matrix_get (m, i, k);
+ double m_kj = gsl_matrix_get (m, k, j);
+ *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
+ : i != k ? -m_ik
+ : j != k ? m_kj
+ : 1.0) / m_kk;
+ }
+ return a;
+ }
+ else
+ {
+ for (size_t i = 0; i < m->size1; i++)
+ {
+ gsl_matrix_set (m, i, k, 0);
+ gsl_matrix_set (m, k, i, 0);
+ }
+ return m;
+ }
+}
+
+static double
+matrix_eval_TRACE (gsl_matrix *m)
+{
+ double sum = 0;
+ size_t n = MIN (m->size1, m->size2);
+ for (size_t i = 0; i < n; i++)
+ sum += gsl_matrix_get (m, i, i);
+ return sum;
+}
+
+static gsl_matrix *
+matrix_eval_T (gsl_matrix *m)
+{
+ return matrix_eval_TRANSPOS (m);
+}
+
+static gsl_matrix *
+matrix_eval_TRANSPOS (gsl_matrix *m)
+{
+ if (m->size1 == m->size2)
+ {
+ gsl_matrix_transpose (m);
+ return m;
+ }
+ else
+ {
+ gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
+ gsl_matrix_transpose_memcpy (t, m);
+ return t;
+ }
+}
+
+static double
+matrix_eval_TRUNC (double d)
+{
+ return trunc (d);
+}
+
+static gsl_matrix *
+matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e)
+{
+ size_t r = r_;
+ size_t c = c_;
+ if (size_overflow_p (xtimes (r, xmax (c, 1))))
+ {
+ struct msg_location *loc = msg_location_dup (e->subs[0]->location);
+ loc->end = e->subs[1]->location->end;
+
+ msg_at (SE, loc,
+ _("Product of arguments to UNIFORM exceeds memory size."));
+
+ msg_location_destroy (loc);
+ return NULL;
+ }
+
+ gsl_matrix *m = gsl_matrix_alloc (r, c);
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = gsl_ran_flat (get_rng (), 0, 1);
+ return m;
+}
+
+static double
+matrix_eval_PDF_BETA (double x, double a, double b)
+{
+ return gsl_ran_beta_pdf (x, a, b);
+}
+
+static double
+matrix_eval_CDF_BETA (double x, double a, double b)
+{
+ return gsl_cdf_beta_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_BETA (double P, double a, double b)
+{
+ return gsl_cdf_beta_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_RV_BETA (double a, double b)
+{
+ return gsl_ran_beta (get_rng (), a, b);
+}
+
+static double
+matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
+{
+ return ncdf_beta (x, a, b, lambda);
+}
+
+static double
+matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
+{
+ return npdf_beta (x, a, b, lambda);
+}
+
+static double
+matrix_eval_CDF_BVNOR (double x0, double x1, double r)
+{
+ return cdf_bvnor (x0, x1, r);
+}
+
+static double
+matrix_eval_PDF_BVNOR (double x0, double x1, double r)
+{
+ return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
+}
+
+static double
+matrix_eval_CDF_CAUCHY (double x, double a, double b)
+{
+ return gsl_cdf_cauchy_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_CAUCHY (double P, double a, double b)
+{
+ return a + b * gsl_cdf_cauchy_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_CAUCHY (double x, double a, double b)
+{
+ return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
+}
+
+static double
+matrix_eval_RV_CAUCHY (double a, double b)
+{
+ return a + b * gsl_ran_cauchy (get_rng (), 1);
+}
+
+static double
+matrix_eval_CDF_CHISQ (double x, double df)
+{
+ return gsl_cdf_chisq_P (x, df);
+}
+
+static double
+matrix_eval_CHICDF (double x, double df)
+{
+ return matrix_eval_CDF_CHISQ (x, df);
+}
+
+static double
+matrix_eval_IDF_CHISQ (double P, double df)
+{
+ return gsl_cdf_chisq_Pinv (P, df);
+}
+
+static double
+matrix_eval_PDF_CHISQ (double x, double df)
+{
+ return gsl_ran_chisq_pdf (x, df);
+}
+
+static double
+matrix_eval_RV_CHISQ (double df)
+{
+ return gsl_ran_chisq (get_rng (), df);
+}
+
+static double
+matrix_eval_SIG_CHISQ (double x, double df)
+{
+ return gsl_cdf_chisq_Q (x, df);
+}
+
+static double
+matrix_eval_CDF_EXP (double x, double a)
+{
+ return gsl_cdf_exponential_P (x, 1. / a);
+}
+
+static double
+matrix_eval_IDF_EXP (double P, double a)
+{
+ return gsl_cdf_exponential_Pinv (P, 1. / a);
+}
+
+static double
+matrix_eval_PDF_EXP (double x, double a)
+{
+ return gsl_ran_exponential_pdf (x, 1. / a);
+}
+
+static double
+matrix_eval_RV_EXP (double a)
+{
+ return gsl_ran_exponential (get_rng (), 1. / a);
+}
+
+static double
+matrix_eval_PDF_XPOWER (double x, double a, double b)
+{
+ return gsl_ran_exppow_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_XPOWER (double a, double b)
+{
+ return gsl_ran_exppow (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_F (double x, double df1, double df2)
+{
+ return gsl_cdf_fdist_P (x, df1, df2);
+}
+
+static double
+matrix_eval_FCDF (double x, double df1, double df2)
+{
+ return matrix_eval_CDF_F (x, df1, df2);
+}
+
+static double
+matrix_eval_IDF_F (double P, double df1, double df2)
+{
+ return idf_fdist (P, df1, df2);
+}
+
+static double
+matrix_eval_RV_F (double df1, double df2)
+{
+ return gsl_ran_fdist (get_rng (), df1, df2);
+}
+
+static double
+matrix_eval_PDF_F (double x, double df1, double df2)
+{
+ return gsl_ran_fdist_pdf (x, df1, df2);
+}
+
+static double
+matrix_eval_SIG_F (double x, double df1, double df2)
+{
+ return gsl_cdf_fdist_Q (x, df1, df2);
+}
+
+static double
+matrix_eval_CDF_GAMMA (double x, double a, double b)
+{
+ return gsl_cdf_gamma_P (x, a, 1. / b);
+}
+
+static double
+matrix_eval_IDF_GAMMA (double P, double a, double b)
+{
+ return gsl_cdf_gamma_Pinv (P, a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_GAMMA (double x, double a, double b)
+{
+ return gsl_ran_gamma_pdf (x, a, 1. / b);
+}
+
+static double
+matrix_eval_RV_GAMMA (double a, double b)
+{
+ return gsl_ran_gamma (get_rng (), a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_LANDAU (double x)
+{
+ return gsl_ran_landau_pdf (x);
+}
+
+static double
+matrix_eval_RV_LANDAU (void)
+{
+ return gsl_ran_landau (get_rng ());
+}
+
+static double
+matrix_eval_CDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_cdf_laplace_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LAPLACE (double P, double a, double b)
+{
+ return a + b * gsl_cdf_laplace_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_ran_laplace_pdf ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_RV_LAPLACE (double a, double b)
+{
+ return a + b * gsl_ran_laplace (get_rng (), 1);
+}
+
+static double
+matrix_eval_RV_LEVY (double c, double alpha)
+{
+ return gsl_ran_levy (get_rng (), c, alpha);
+}
+
+static double
+matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
+{
+ return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
+}
+
+static double
+matrix_eval_CDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_cdf_logistic_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LOGISTIC (double P, double a, double b)
+{
+ return a + b * gsl_cdf_logistic_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
+}
+
+static double
+matrix_eval_RV_LOGISTIC (double a, double b)
+{
+ return a + b * gsl_ran_logistic (get_rng (), 1);
+}
+
+static double
+matrix_eval_CDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_cdf_lognormal_P (x, log (m), s);
+}
+
+static double
+matrix_eval_IDF_LNORMAL (double P, double m, double s)
+{
+ return gsl_cdf_lognormal_Pinv (P, log (m), s);;
+}
+
+static double
+matrix_eval_PDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_ran_lognormal_pdf (x, log (m), s);
+}
+
+static double
+matrix_eval_RV_LNORMAL (double m, double s)
+{
+ return gsl_ran_lognormal (get_rng (), log (m), s);
+}
+
+static double
+matrix_eval_CDF_NORMAL (double x, double u, double s)
+{
+ return gsl_cdf_gaussian_P (x - u, s);
+}
+
+static double
+matrix_eval_IDF_NORMAL (double P, double u, double s)
+{
+ return u + gsl_cdf_gaussian_Pinv (P, s);
+}
+
+static double
+matrix_eval_PDF_NORMAL (double x, double u, double s)
+{
+ return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
+}
+
+static double
+matrix_eval_RV_NORMAL (double u, double s)
+{
+ return u + gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_CDFNORM (double x)
+{
+ return gsl_cdf_ugaussian_P (x);
+}
+
+static double
+matrix_eval_PROBIT (double P)
+{
+ return gsl_cdf_ugaussian_Pinv (P);
+}
+
+static double
+matrix_eval_NORMAL (double s)
+{
+ return gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_PDF_NTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
+}
+
+static double
+matrix_eval_RV_NTAIL (double a, double sigma)
+{
+ return gsl_ran_gaussian_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_PARETO (double x, double a, double b)
+{
+ return gsl_cdf_pareto_P (x, b, a);
+}
+
+static double
+matrix_eval_IDF_PARETO (double P, double a, double b)
+{
+ return gsl_cdf_pareto_Pinv (P, b, a);
+}
+
+static double
+matrix_eval_PDF_PARETO (double x, double a, double b)
+{
+ return gsl_ran_pareto_pdf (x, b, a);
+}
+
+static double
+matrix_eval_RV_PARETO (double a, double b)
+{
+ return gsl_ran_pareto (get_rng (), b, a);
+}
+
+static double
+matrix_eval_CDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_cdf_rayleigh_P (x, sigma);
+}
+
+static double
+matrix_eval_IDF_RAYLEIGH (double P, double sigma)
+{
+ return gsl_cdf_rayleigh_Pinv (P, sigma);
+}
+
+static double
+matrix_eval_PDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_ran_rayleigh_pdf (x, sigma);
+}
+
+static double
+matrix_eval_RV_RAYLEIGH (double sigma)
+{
+ return gsl_ran_rayleigh (get_rng (), sigma);
+}
+
+static double
+matrix_eval_PDF_RTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
+}
+
+static double
+matrix_eval_RV_RTAIL (double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_T (double x, double df)
+{
+ return gsl_cdf_tdist_P (x, df);
+}
+
+static double
+matrix_eval_TCDF (double x, double df)
+{
+ return matrix_eval_CDF_T (x, df);
+}
+
+static double
+matrix_eval_IDF_T (double P, double df)
+{
+ return gsl_cdf_tdist_Pinv (P, df);
+}
+
+static double
+matrix_eval_PDF_T (double x, double df)
+{
+ return gsl_ran_tdist_pdf (x, df);
+}
+
+static double
+matrix_eval_RV_T (double df)
+{
+ return gsl_ran_tdist (get_rng (), df);
+}
+
+static double
+matrix_eval_CDF_T1G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T1G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T1G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T1G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_T2G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T2G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T2G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T2G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_cdf_flat_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_UNIFORM (double P, double a, double b)
+{
+ return gsl_cdf_flat_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_ran_flat_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_UNIFORM (double a, double b)
+{
+ return gsl_ran_flat (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_cdf_weibull_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_WEIBULL (double P, double a, double b)
+{
+ return gsl_cdf_weibull_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_ran_weibull_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_WEIBULL (double a, double b)
+{
+ return gsl_ran_weibull (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_BERNOULLI (double k, double p)
+{
+ return k ? 1 : 1 - p;
+}
+
+static double
+matrix_eval_PDF_BERNOULLI (double k, double p)
+{
+ return gsl_ran_bernoulli_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_BERNOULLI (double p)
+{
+ return gsl_ran_bernoulli (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_BINOM (double k, double n, double p)
+{
+ return gsl_cdf_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_BINOM (double k, double n, double p)
+{
+ return gsl_ran_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_BINOM (double n, double p)
+{
+ return gsl_ran_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_GEOM (double k, double p)
+{
+ return gsl_cdf_geometric_P (k, p);
+}
+
+static double
+matrix_eval_PDF_GEOM (double k, double p)
+{
+ return gsl_ran_geometric_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_GEOM (double p)
+{
+ return gsl_ran_geometric (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_cdf_hypergeometric_P (k, c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
+}
+
+static double
+matrix_eval_RV_HYPER (double a, double b, double c)
+{
+ return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_LOG (double k, double p)
+{
+ return gsl_ran_logarithmic_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_LOG (double p)
+{
+ return gsl_ran_logarithmic (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_cdf_negative_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_ran_negative_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_NEGBIN (double n, double p)
+{
+ return gsl_ran_negative_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_POISSON (double k, double mu)
+{
+ return gsl_cdf_poisson_P (k, mu);
+}
+
+static double
+matrix_eval_PDF_POISSON (double k, double mu)
+{
+ return gsl_ran_poisson_pdf (k, mu);
+}
+
+static double
+matrix_eval_RV_POISSON (double mu)
+{
+ return gsl_ran_poisson (get_rng (), mu);
+}
+
+static double
+matrix_op_eval (enum matrix_op op, double a, double b)
+{
+ switch (op)
+ {
+ case MOP_ADD_ELEMS: return a + b;
+ case MOP_SUB_ELEMS: return a - b;
+ case MOP_MUL_ELEMS: return a * b;
+ case MOP_DIV_ELEMS: return a / b;
+ case MOP_EXP_ELEMS: return pow (a, b);
+ case MOP_GT: return a > b;
+ case MOP_GE: return a >= b;
+ case MOP_LT: return a < b;
+ case MOP_LE: return a <= b;
+ case MOP_EQ: return a == b;
+ case MOP_NE: return a != b;
+ case MOP_AND: return (a > 0) && (b > 0);
+ case MOP_OR: return (a > 0) || (b > 0);
+ case MOP_XOR: return (a > 0) != (b > 0);
+
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
+ MATRIX_FUNCTIONS
+#undef F
+ case MOP_NEGATE:
+ case MOP_SEQ:
+ case MOP_SEQ_BY:
+ case MOP_MUL_MAT:
+ case MOP_EXP_MAT:
+ case MOP_NOT:
+ case MOP_PASTE_HORZ:
+ case MOP_PASTE_VERT:
+ case MOP_EMPTY:
+ case MOP_VEC_INDEX:
+ case MOP_VEC_ALL:
+ case MOP_MAT_INDEX:
+ case MOP_ROW_INDEX:
+ case MOP_COL_INDEX:
+ case MOP_NUMBER:
+ case MOP_VARIABLE:
+ case MOP_EOF:
+ NOT_REACHED ();
+ }
+ NOT_REACHED ();
+}
+
+static const char *
+matrix_op_name (enum matrix_op op)
+{
+ switch (op)
+ {
+ case MOP_ADD_ELEMS: return "+";
+ case MOP_SUB_ELEMS: return "-";
+ case MOP_MUL_ELEMS: return "&*";
+ case MOP_DIV_ELEMS: return "&/";
+ case MOP_EXP_ELEMS: return "&**";
+ case MOP_GT: return ">";
+ case MOP_GE: return ">=";
+ case MOP_LT: return "<";
+ case MOP_LE: return "<=";
+ case MOP_EQ: return "=";
+ case MOP_NE: return "<>";
+ case MOP_AND: return "AND";
+ case MOP_OR: return "OR";
+ case MOP_XOR: return "XOR";
+
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
+ MATRIX_FUNCTIONS
+#undef F
+ case MOP_NEGATE:
+ case MOP_SEQ:
+ case MOP_SEQ_BY:
+ case MOP_MUL_MAT:
+ case MOP_EXP_MAT:
+ case MOP_NOT:
+ case MOP_PASTE_HORZ:
+ case MOP_PASTE_VERT:
+ case MOP_EMPTY:
+ case MOP_VEC_INDEX:
+ case MOP_VEC_ALL:
+ case MOP_MAT_INDEX:
+ case MOP_ROW_INDEX:
+ case MOP_COL_INDEX:
+ case MOP_NUMBER:
+ case MOP_VARIABLE:
+ case MOP_EOF:
+ NOT_REACHED ();
+ }
+ NOT_REACHED ();
+}
+
+static bool
+is_scalar (const gsl_matrix *m)
+{
+ return m->size1 == 1 && m->size2 == 1;
+}
+
+static double
+to_scalar (const gsl_matrix *m)
+{
+ assert (is_scalar (m));
+ return gsl_matrix_get (m, 0, 0);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_elementwise (const struct matrix_expr *e,
+ enum matrix_op op,
+ gsl_matrix *a, gsl_matrix *b)
+{
+ if (is_scalar (b))
+ {
+ double be = to_scalar (b);
+ for (size_t r = 0; r < a->size1; r++)
+ for (size_t c = 0; c < a->size2; c++)
+ {
+ double *ae = gsl_matrix_ptr (a, r, c);
+ *ae = matrix_op_eval (op, *ae, be);
+ }
+ return a;
+ }
+ else if (is_scalar (a))
+ {
+ double ae = to_scalar (a);
+ for (size_t r = 0; r < b->size1; r++)
+ for (size_t c = 0; c < b->size2; c++)
+ {
+ double *be = gsl_matrix_ptr (b, r, c);
+ *be = matrix_op_eval (op, ae, *be);
+ }
+ return b;
+ }
+ else if (a->size1 == b->size1 && a->size2 == b->size2)
+ {
+ for (size_t r = 0; r < a->size1; r++)
+ for (size_t c = 0; c < a->size2; c++)
+ {
+ double *ae = gsl_matrix_ptr (a, r, c);
+ double be = gsl_matrix_get (b, r, c);
+ *ae = matrix_op_eval (op, *ae, be);
+ }
+ return a;
+ }
+ else
+ {
+ msg_at (SE, matrix_expr_location (e),
+ _("The operands of %s must have the same dimensions or one "
+ "must be a scalar."),
+ matrix_op_name (op));
+ msg_at (SN, matrix_expr_location (e->subs[0]),
+ _("The left-hand operand is a %zu×%zu matrix."),
+ a->size1, a->size2);
+ msg_at (SN, matrix_expr_location (e->subs[1]),
+ _("The right-hand operand is a %zu×%zu matrix."),
+ b->size1, b->size2);
+ return NULL;
+ }
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_mul_mat (const struct matrix_expr *e,
+ gsl_matrix *a, gsl_matrix *b)
+{
+ if (is_scalar (a) || is_scalar (b))
+ return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b);
+
+ if (a->size2 != b->size1)
+ {
+ msg_at (SE, e->location,
+ _("Matrices not conformable for multiplication."));
+ msg_at (SN, matrix_expr_location (e->subs[0]),
+ _("The left-hand operand is a %zu×%zu matrix."),
+ a->size1, a->size2);
+ msg_at (SN, matrix_expr_location (e->subs[1]),
+ _("The right-hand operand is a %zu×%zu matrix."),
+ b->size1, b->size2);
+ return NULL;
+ }
+
+ gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
+ gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
+ return c;
+}
+
+static void
+swap_matrix (gsl_matrix **a, gsl_matrix **b)
+{
+ gsl_matrix *tmp = *a;
+ *a = *b;
+ *b = tmp;
+}
+
+static void
+mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
+ gsl_matrix **tmp)
+{
+ gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
+ swap_matrix (z, tmp);
+}
+
+static void
+square_matrix (gsl_matrix **x, gsl_matrix **tmp)
+{
+ mul_matrix (x, *x, *x, tmp);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
+ gsl_matrix *x_, gsl_matrix *b)
+{
+ gsl_matrix *x = x_;
+ if (x->size1 != x->size2)
+ {
+ msg_at (SE, matrix_expr_location (e->subs[0]),
+ _("Matrix exponentation with ** requires a square matrix on "
+ "the left-hand size, not one with dimensions %zu×%zu."),
+ x->size1, x->size2);
+ return NULL;
+ }
+ if (!is_scalar (b))
+ {
+ msg_at (SE, matrix_expr_location (e->subs[1]),
+ _("Matrix exponentiation with ** requires a scalar on the "
+ "right-hand side, not a matrix with dimensions %zu×%zu."),
+ b->size1, b->size2);
+ return NULL;
+ }
+ double bf = to_scalar (b);
+ if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
+ {
+ msg_at (SE, matrix_expr_location (e->subs[1]),
+ _("Exponent %.1f in matrix exponentiation is non-integer "
+ "or outside the valid range."), bf);
+ return NULL;
+ }
+ long int bl = bf;
+
+ gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
+ gsl_matrix *y = y_;
+ gsl_matrix_set_identity (y);
+ if (bl == 0)
+ return y;
+
+ gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
+ gsl_matrix *t = t_;
+ for (unsigned long int n = labs (bl); n > 1; n /= 2)
+ if (n & 1)
+ {
+ mul_matrix (&y, x, y, &t);
+ square_matrix (&x, &t);
+ }
+ else
+ square_matrix (&x, &t);
+
+ mul_matrix (&y, x, y, &t);
+ if (bf < 0)
+ {
+ invert_matrix (y, x);
+ swap_matrix (&x, &y);
+ }
+
+ /* Garbage collection.
+
+ There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
+ a permutation of them. We are returning one of them; that one must not be
+ destroyed. We must not destroy 'x_' because the caller owns it. */
+ if (y != y_)
+ gsl_matrix_free (y_);
+ if (y != t_)
+ gsl_matrix_free (t_);
+
+ return y;
+}
+
+static void
+note_operand_size (const gsl_matrix *m, const struct matrix_expr *e)
+{
+ msg_at (SN, matrix_expr_location (e),
+ _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
+}
+
+static void
+note_nonscalar (const gsl_matrix *m, const struct matrix_expr *e)
+{
+ if (!is_scalar (m))
+ note_operand_size (m, e);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_seq (const struct matrix_expr *e,
+ gsl_matrix *start_, gsl_matrix *end_,
+ gsl_matrix *by_)
+{
+ if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
+ {
+ msg_at (SE, matrix_expr_location (e),
+ _("All operands of : operator must be scalars."));
+
+ note_nonscalar (start_, e->subs[0]);
+ note_nonscalar (end_, e->subs[1]);
+ if (by_)
+ note_nonscalar (by_, e->subs[2]);
+ return NULL;
+ }
+
+ long int start = to_scalar (start_);
+ long int end = to_scalar (end_);
+ long int by = by_ ? to_scalar (by_) : 1;
+
+ if (!by)
+ {
+ msg_at (SE, matrix_expr_location (e->subs[2]),
+ _("The increment operand to : must be nonzero."));
+ return NULL;
+ }
+
+ long int n = (end >= start && by > 0 ? (end - start + by) / by
+ : end <= start && by < 0 ? (start - end - by) / -by
+ : 0);
+ gsl_matrix *m = gsl_matrix_alloc (1, n);
+ for (long int i = 0; i < n; i++)
+ gsl_matrix_set (m, 0, i, start + i * by);
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_not (gsl_matrix *a)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, a)
+ *d = !(*d > 0);
+ return a;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_paste_horz (const struct matrix_expr *e,
+ gsl_matrix *a, gsl_matrix *b)
+{
+ if (a->size1 != b->size1)
+ {
+ if (!a->size1 || !a->size2)
+ return b;
+ else if (!b->size1 || !b->size2)
+ return a;
+
+ msg_at (SE, matrix_expr_location (e),
+ _("This expression tries to horizontally join matrices with "
+ "differing numbers of rows."));
+ note_operand_size (a, e->subs[0]);
+ note_operand_size (b, e->subs[1]);
+ return NULL;
+ }
+
+ gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
+ for (size_t y = 0; y < a->size1; y++)
+ {
+ for (size_t x = 0; x < a->size2; x++)
+ gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
+ for (size_t x = 0; x < b->size2; x++)
+ gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
+ }
+ return c;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_paste_vert (const struct matrix_expr *e,
+ gsl_matrix *a, gsl_matrix *b)
+{
+ if (a->size2 != b->size2)
+ {
+ if (!a->size1 || !a->size2)
+ return b;
+ else if (!b->size1 || !b->size2)
+ return a;
+
+ msg_at (SE, matrix_expr_location (e),
+ _("This expression tries to vertically join matrices with "
+ "differing numbers of columns."));
+ note_operand_size (a, e->subs[0]);
+ note_operand_size (b, e->subs[1]);
+ return NULL;
+ }
+
+ gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
+ for (size_t x = 0; x < a->size2; x++)
+ {
+ for (size_t y = 0; y < a->size1; y++)
+ gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
+ for (size_t y = 0; y < b->size1; y++)
+ gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
+ }
+ return c;
+}
+
+static gsl_vector *
+matrix_to_vector (gsl_matrix *m)
+{
+ assert (m->owner);
+ gsl_vector v = to_vector (m);
+ assert (v.block == m->block || !v.block);
+ assert (!v.owner);
+ v.owner = 1;
+ m->owner = 0;
+ gsl_matrix_free (m);
+ return xmemdup (&v, sizeof v);
+}
+
+enum index_type {
+ IV_ROW,
+ IV_COLUMN,
+ IV_VECTOR
+};
+
+struct index_vector
+ {
+ size_t *indexes;
+ size_t n;
+ };
+#define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
+
+static void
+index_vector_uninit (struct index_vector *iv)
+{
+ if (iv)
+ free (iv->indexes);
+}
+
+static bool
+matrix_normalize_index_vector (const gsl_matrix *m,
+ const struct matrix_expr *me, size_t size,
+ enum index_type index_type, size_t other_size,
+ struct index_vector *iv)
+{
+ if (m)
+ {
+ if (!is_vector (m))
+ {
+ switch (index_type)
+ {
+ case IV_VECTOR:
+ msg_at (SE, matrix_expr_location (me),
+ _("Vector index must be scalar or vector, not a "
+ "%zu×%zu matrix."),
+ m->size1, m->size2);
+ break;
+
+ case IV_ROW:
+ msg_at (SE, matrix_expr_location (me),
+ _("Matrix row index must be scalar or vector, not a "
+ "%zu×%zu matrix."),
+ m->size1, m->size2);
+ break;
+
+ case IV_COLUMN:
+ msg_at (SE, matrix_expr_location (me),
+ _("Matrix column index must be scalar or vector, not a "
+ "%zu×%zu matrix."),
+ m->size1, m->size2);
+ break;
+ }
+ return false;
+ }
+
+ gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
+ *iv = (struct index_vector) {
+ .indexes = xnmalloc (v.size, sizeof *iv->indexes),
+ .n = v.size,
+ };
+ for (size_t i = 0; i < v.size; i++)
+ {
+ double index = gsl_vector_get (&v, i);
+ if (index < 1 || index >= size + 1)
+ {
+ switch (index_type)
+ {
+ case IV_VECTOR:
+ msg_at (SE, matrix_expr_location (me),
+ _("Index %g is out of range for vector "
+ "with %zu elements."), index, size);
+ break;
+
+ case IV_ROW:
+ msg_at (SE, matrix_expr_location (me),
+ _("%g is not a valid row index for "
+ "a %zu×%zu matrix."),
+ index, size, other_size);
+ break;
+
+ case IV_COLUMN:
+ msg_at (SE, matrix_expr_location (me),
+ _("%g is not a valid column index for "
+ "a %zu×%zu matrix."),
+ index, other_size, size);
+ break;
+ }
+
+ index_vector_uninit (iv);
+ return false;
+ }
+ iv->indexes[i] = index - 1;
+ }
+ return true;
+ }
+ else
+ {
+ *iv = (struct index_vector) {
+ .indexes = xnmalloc (size, sizeof *iv->indexes),
+ .n = size,
+ };
+ for (size_t i = 0; i < size; i++)
+ iv->indexes[i] = i;
+ return true;
+ }
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_vec_all (const struct matrix_expr *e,
+ gsl_matrix *sm)
+{
+ if (!is_vector (sm))
+ {
+ msg_at (SE, matrix_expr_location (e->subs[0]),
+ _("Vector index operator may not be applied to "
+ "a %zu×%zu matrix."),
+ sm->size1, sm->size2);
+ return NULL;
+ }
+
+ return sm;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_vec_index (const struct matrix_expr *e,
+ gsl_matrix *sm, gsl_matrix *im)
+{
+ if (!matrix_expr_evaluate_vec_all (e, sm))
+ return NULL;
+
+ gsl_vector sv = to_vector (sm);
+ struct index_vector iv;
+ if (!matrix_normalize_index_vector (im, e->subs[1],
+ sv.size, IV_VECTOR, 0, &iv))
+ return NULL;
+
+ gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
+ sm->size1 == 1 ? iv.n : 1);
+ gsl_vector dv = to_vector (dm);
+ for (size_t dx = 0; dx < iv.n; dx++)
+ {
+ size_t sx = iv.indexes[dx];
+ gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
+ }
+ index_vector_uninit (&iv);
+
+ return dm;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_mat_index (gsl_matrix *sm,
+ gsl_matrix *im0, const struct matrix_expr *eim0,
+ gsl_matrix *im1, const struct matrix_expr *eim1)
+{
+ struct index_vector iv0;
+ if (!matrix_normalize_index_vector (im0, eim0, sm->size1,
+ IV_ROW, sm->size2, &iv0))
+ return NULL;
+
+ struct index_vector iv1;
+ if (!matrix_normalize_index_vector (im1, eim1, sm->size2,
+ IV_COLUMN, sm->size1, &iv1))
+ {
+ index_vector_uninit (&iv0);
+ return NULL;
+ }
+
+ gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
+ for (size_t dy = 0; dy < iv0.n; dy++)
+ {
+ size_t sy = iv0.indexes[dy];
+
+ for (size_t dx = 0; dx < iv1.n; dx++)
+ {
+ size_t sx = iv1.indexes[dx];
+ gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
+ }
+ }
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return dm;
+}
+
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+ static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
+ const struct matrix_function_properties *, gsl_matrix *[], \
+ const struct matrix_expr *, matrix_proto_##PROTO *);
+MATRIX_FUNCTIONS
+#undef F
+
+static bool
+check_scalar_arg (const char *name, gsl_matrix *subs[],
+ const struct matrix_expr *e, size_t index)
+{
+ if (!is_scalar (subs[index]))
+ {
+ msg_at (SE, matrix_expr_location (e->subs[index]),
+ _("Function %s argument %zu must be a scalar, "
+ "not a %zu×%zu matrix."),
+ name, index + 1, subs[index]->size1, subs[index]->size2);
+ return false;
+ }
+ return true;
+}
+
+static bool
+check_vector_arg (const char *name, gsl_matrix *subs[],
+ const struct matrix_expr *e, size_t index)
+{
+ if (!is_vector (subs[index]))
+ {
+ msg_at (SE, matrix_expr_location (e->subs[index]),
+ _("Function %s argument %zu must be a vector, "
+ "not a %zu×%zu matrix."),
+ name, index + 1, subs[index]->size1, subs[index]->size2);
+ return false;
+ }
+ return true;
+}
+
+static bool
+to_scalar_args (const char *name, gsl_matrix *subs[],
+ const struct matrix_expr *e, double d[])
+{
+ for (size_t i = 0; i < e->n_subs; i++)
+ {
+ if (!check_scalar_arg (name, subs, e, i))
+ return false;
+ d[i] = to_scalar (subs[i]);
+ }
+ return true;
+}
+
+static int
+parse_constraint_value (const char **constraintsp)
+{
+ char *tail;
+ long retval = strtol (*constraintsp, &tail, 10);
+ assert (tail > *constraintsp);
+ *constraintsp = tail;
+ return retval;
+}
+
+enum matrix_argument_relop
+ {
+ MRR_GT, /* > */
+ MRR_GE, /* >= */
+ MRR_LT, /* < */
+ MRR_LE, /* <= */
+ MRR_NE, /* <> */
+ };
+
+static void
+argument_inequality_error (
+ const struct matrix_function_properties *props, const struct matrix_expr *e,
+ size_t ai, gsl_matrix *a, size_t y, size_t x,
+ size_t bi, double b,
+ enum matrix_argument_relop relop)
+{
+ const struct msg_location *loc = matrix_expr_location (e);
+ switch (relop)
+ {
+ case MRR_GE:
+ msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
+ "than or equal to argument %zu."),
+ ai + 1, props->name, bi + 1);
+ break;
+
+ case MRR_GT:
+ msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
+ "than argument %zu."),
+ ai + 1, props->name, bi + 1);
+ break;
+
+ case MRR_LE:
+ msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
+ "or equal to argument %zu."),
+ ai + 1, props->name, bi + 1);
+ break;
+
+ case MRR_LT:
+ msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
+ "argument %zu."),
+ ai + 1, props->name, bi + 1);
+ break;
+
+ case MRR_NE:
+ msg_at (ME, loc, _("Argument %zu to matrix function %s must not be equal "
+ "to argument %zu."),
+ ai + 1, props->name, bi + 1);
+ break;
+ }
+
+ const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
+ if (is_scalar (a))
+ msg_at (SN, a_loc, _("Argument %zu is %g."),
+ ai + 1, gsl_matrix_get (a, y, x));
+ else
+ msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
+ y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
+
+ msg_at (SN, matrix_expr_location (e->subs[bi]),
+ _("Argument %zu is %g."), bi + 1, b);
+}
+
+static void
+argument_value_error (
+ const struct matrix_function_properties *props, const struct matrix_expr *e,
+ size_t ai, gsl_matrix *a, size_t y, size_t x,
+ double b,
+ enum matrix_argument_relop relop)
+{
+ const struct msg_location *loc = matrix_expr_location (e);
+ switch (relop)
+ {
+ case MRR_GE:
+ msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
+ "than or equal to %g."),
+ ai + 1, props->name, b);
+ break;
+
+ case MRR_GT:
+ msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
+ "than %g."),
+ ai + 1, props->name, b);
+ break;
+
+ case MRR_LE:
+ msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
+ "or equal to %g."),
+ ai + 1, props->name, b);
+ break;
+
+ case MRR_LT:
+ msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
+ "%g."),
+ ai + 1, props->name, b);
+ break;
+
+ case MRR_NE:
+ msg_at (SE, loc, _("Argument %zu to matrix function %s must not be equal "
+ "to %g."),
+ ai + 1, props->name, b);
+ break;
+ }
+
+ const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
+ if (is_scalar (a))
+ {
+ if (relop != MRR_NE)
+ msg_at (SN, a_loc, _("Argument %zu is %g."),
+ ai + 1, gsl_matrix_get (a, y, x));
+ }
+ else
+ msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
+ y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
+}
+
+static bool
+matrix_argument_relop_is_satisfied (double a, double b,
+ enum matrix_argument_relop relop)
+{
+ switch (relop)
+ {
+ case MRR_GE: return a >= b;
+ case MRR_GT: return a > b;
+ case MRR_LE: return a <= b;
+ case MRR_LT: return a < b;
+ case MRR_NE: return a != b;
+ }
+
+ NOT_REACHED ();
+}
+
+static enum matrix_argument_relop
+matrix_argument_relop_flip (enum matrix_argument_relop relop)
+{
+ switch (relop)
+ {
+ case MRR_GE: return MRR_LE;
+ case MRR_GT: return MRR_LT;
+ case MRR_LE: return MRR_GE;
+ case MRR_LT: return MRR_GT;
+ case MRR_NE: return MRR_NE;
+ }
+
+ NOT_REACHED ();
+}
+
+static bool
+check_constraints (const struct matrix_function_properties *props,
+ gsl_matrix *args[], const struct matrix_expr *e)
+{
+ size_t n_args = e->n_subs;
+ const char *constraints = props->constraints;
+ if (!constraints)
+ return true;
+
+ size_t arg_index = SIZE_MAX;
+ while (*constraints)
+ {
+ if (*constraints >= 'a' && *constraints <= 'd')
+ {
+ arg_index = *constraints++ - 'a';
+ assert (arg_index < n_args);
+ }
+ else if (*constraints == '[' || *constraints == '(')
+ {
+ assert (arg_index < n_args);
+ bool open_lower = *constraints++ == '(';
+ int minimum = parse_constraint_value (&constraints);
+ assert (*constraints == ',');
+ constraints++;
+ int maximum = parse_constraint_value (&constraints);
+ assert (*constraints == ']' || *constraints == ')');
+ bool open_upper = *constraints++ == ')';
+
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+ if ((open_lower ? *d <= minimum : *d < minimum)
+ || (open_upper ? *d >= maximum : *d > maximum))
+ {
+ if (!is_scalar (args[arg_index]))
+ msg_at (SE, matrix_expr_location (e->subs[arg_index]),
+ _("Row %zu, column %zu of argument %zu to matrix "
+ "function %s is %g, which is outside "
+ "the valid range %c%d,%d%c."),
+ y + 1, x + 1, arg_index + 1, props->name, *d,
+ open_lower ? '(' : '[',
+ minimum, maximum,
+ open_upper ? ')' : ']');
+ else
+ msg_at (SE, matrix_expr_location (e->subs[arg_index]),
+ _("Argument %zu to matrix function %s is %g, "
+ "which is outside the valid range %c%d,%d%c."),
+ arg_index + 1, props->name, *d,
+ open_lower ? '(' : '[',
+ minimum, maximum,
+ open_upper ? ')' : ']');
+ return false;
+ }
+ }
+ else if (*constraints == 'i')
+ {
+ constraints++;
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+ if (*d != floor (*d))
+ {
+ if (!is_scalar (args[arg_index]))
+ msg_at (SE, matrix_expr_location (e->subs[arg_index]),
+ _("Argument %zu to matrix function %s, which must be "
+ "integer, contains non-integer value %g in "
+ "row %zu, column %zu."),
+ arg_index + 1, props->name, *d, y + 1, x + 1);
+ else
+ msg_at (SE, matrix_expr_location (e->subs[arg_index]),
+ _("Argument %zu to matrix function %s, which must be "
+ "integer, has non-integer value %g."),
+ arg_index + 1, props->name, *d);
+ return false;
+ }
+ }
+ else if (*constraints == '>'
+ || *constraints == '<'
+ || *constraints == '!')
+ {
+ enum matrix_argument_relop relop;
+ switch (*constraints++)
+ {
+ case '>':
+ if (*constraints == '=')
+ {
+ constraints++;
+ relop = MRR_GE;
+ }
+ else
+ relop = MRR_GT;
+ break;
+
+ case '<':
+ if (*constraints == '=')
+ {
+ constraints++;
+ relop = MRR_LE;
+ }
+ else
+ relop = MRR_LT;
+ break;
+
+ case '!':
+ assert (*constraints == '=');
+ constraints++;
+ relop = MRR_NE;
+ break;
+
+ default:
+ NOT_REACHED ();
+ }
+
+ if (*constraints >= 'a' && *constraints <= 'd')
+ {
+ size_t a_index = arg_index;
+ size_t b_index = *constraints - 'a';
+ assert (a_index < n_args);
+ assert (b_index < n_args);
+
+ /* We only support one of the two arguments being non-scalar.
+ It's easier to support only the first one being non-scalar, so
+ flip things around if it's the other way. */
+ if (!is_scalar (args[b_index]))
+ {
+ assert (is_scalar (args[a_index]));
+ size_t tmp_index = a_index;
+ a_index = b_index;
+ b_index = tmp_index;
+ relop = matrix_argument_relop_flip (relop);
+ }
+
+ double b = to_scalar (args[b_index]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
+ if (!matrix_argument_relop_is_satisfied (*a, b, relop))
+ {
+ argument_inequality_error (
+ props, e,
+ a_index, args[a_index], y, x,
+ b_index, b,
+ relop);
+ return false;
+ }
+ }
+ else
+ {
+ int comparand = parse_constraint_value (&constraints);
+
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+ if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
+ {
+ argument_value_error (
+ props, e,
+ arg_index, args[arg_index], y, x,
+ comparand,
+ relop);
+ return false;
+ }
+ }
+ }
+ else
+ {
+ assert (*constraints == ' ');
+ constraints++;
+ arg_index = SIZE_MAX;
+ }
+ }
+ return true;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_none (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_d_none *f)
+{
+ assert (e->n_subs == 0);
+
+ if (!check_constraints (props, subs, e))
+ return NULL;
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f ());
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_d_d *f)
+{
+ assert (e->n_subs == 1);
+
+ double d;
+ if (!to_scalar_args (props->name, subs, e, &d)
+ || !check_constraints (props, subs, e))
+ return NULL;
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f (d));
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_d_dd *f)
+{
+ assert (e->n_subs == 2);
+
+ double d[2];
+ if (!to_scalar_args (props->name, subs, e, d)
+ && !check_constraints (props, subs, e))
+ return NULL;
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_d_ddd *f)
+{
+ assert (e->n_subs == 3);
+
+ double d[3];
+ if (!to_scalar_args (props->name, subs, e, d)
+ || !check_constraints (props, subs, e))
+ return NULL;
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_d *f)
+{
+ assert (e->n_subs == 1);
+
+ double d;
+ return (to_scalar_args (props->name, subs, e, &d)
+ && check_constraints (props, subs, e)
+ ? f(d)
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_ddd *f)
+{
+ assert (e->n_subs == 3);
+
+ double d[3];
+ return (to_scalar_args (props->name, subs, e, d)
+ && check_constraints (props, subs, e)
+ ? f(d[0], d[1], d[2])
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_ddn *f)
+{
+ assert (e->n_subs == 2);
+
+ double d[2];
+ return (to_scalar_args (props->name, subs, e, d)
+ && check_constraints (props, subs, e)
+ ? f(d[0], d[1], e)
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_m (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_m *f)
+{
+ assert (e->n_subs == 1);
+ return check_constraints (props, subs, e) ? f (subs[0]) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_mn *f)
+{
+ assert (e->n_subs == 1);
+ return check_constraints (props, subs, e) ? f (subs[0], e) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_e *f)
+{
+ assert (e->n_subs == 1);
+
+ if (!check_constraints (props, subs, e))
+ return NULL;
+
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a);
+ return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_md (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_md *f)
+{
+ assert (e->n_subs == 2);
+ return (check_scalar_arg (props->name, subs, e, 1)
+ && check_constraints (props, subs, e)
+ ? f (subs[0], to_scalar (subs[1]))
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_mdn *f)
+{
+ assert (e->n_subs == 2);
+ return (check_scalar_arg (props->name, subs, e, 1)
+ && check_constraints (props, subs, e)
+ ? f (subs[0], to_scalar (subs[1]), e)
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_ed *f)
+{
+ assert (e->n_subs == 2);
+ if (!check_scalar_arg (props->name, subs, e, 1)
+ || !check_constraints (props, subs, e))
+ return NULL;
+
+ double b = to_scalar (subs[1]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a, b);
+ return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_mddn *f)
+{
+ assert (e->n_subs == 3);
+ if (!check_scalar_arg (props->name, subs, e, 1)
+ || !check_scalar_arg (props->name, subs, e, 2)
+ || !check_constraints (props, subs, e))
+ return NULL;
+ return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_edd *f)
+{
+ assert (e->n_subs == 3);
+ if (!check_scalar_arg (props->name, subs, e, 1)
+ || !check_scalar_arg (props->name, subs, e, 2)
+ || !check_constraints (props, subs, e))
+ return NULL;
+
+ double b = to_scalar (subs[1]);
+ double c = to_scalar (subs[2]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a, b, c);
+ return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_eddd *f)
+{
+ assert (e->n_subs == 4);
+ for (size_t i = 1; i < 4; i++)
+ if (!check_scalar_arg (props->name, subs, e, i))
+ return NULL;
+
+ if (!check_constraints (props, subs, e))
+ return NULL;
+
+ double b = to_scalar (subs[1]);
+ double c = to_scalar (subs[2]);
+ double d = to_scalar (subs[3]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a, b, c, d);
+ return subs[0];
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_eed *f)
+{
+ assert (e->n_subs == 3);
+ if (!check_scalar_arg (props->name, subs, e, 2))
+ return NULL;
+
+ if (!is_scalar (subs[0]) && !is_scalar (subs[1])
+ && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
+ {
+ struct msg_location *loc = msg_location_dup (e->subs[0]->location);
+ loc->end = e->subs[1]->location->end;
+
+ msg_at (ME, loc,
+ _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
+ "%zu×%zu, but %s requires these arguments either to have "
+ "the same dimensions or for one of them to be a scalar."),
+ props->name,
+ subs[0]->size1, subs[0]->size2,
+ subs[1]->size1, subs[1]->size2,
+ props->name);
+
+ msg_location_destroy (loc);
+ return NULL;
+ }
+
+ if (!check_constraints (props, subs, e))
+ return NULL;
+
+ double c = to_scalar (subs[2]);
+
+ if (is_scalar (subs[0]))
+ {
+ double a = to_scalar (subs[0]);
+ MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
+ *b = f (a, *b, c);
+ return subs[1];
+ }
+ else
+ {
+ double b = to_scalar (subs[1]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a, b, c);
+ return subs[0];
+ }
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_mm *f)
+{
+ assert (e->n_subs == 2);
+ return check_constraints (props, subs, e) ? f (subs[0], subs[1]) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_mmn *f)
+{
+ assert (e->n_subs == 2);
+ return check_constraints (props, subs, e) ? f (subs[0], subs[1], e) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_v *f)
+{
+ assert (e->n_subs == 1);
+ if (!check_vector_arg (props->name, subs, e, 0)
+ || !check_constraints (props, subs, e))
+ return NULL;
+ gsl_vector v = to_vector (subs[0]);
+ return f (&v);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_m (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_d_m *f)
+{
+ assert (e->n_subs == 1);
+
+ if (!check_constraints (props, subs, e))
+ return NULL;
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f (subs[0]));
+ return m;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_any (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_m_any *f)
+{
+ return check_constraints (props, subs, e) ? f (subs, e->n_subs) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props_ UNUSED,
+ gsl_matrix *subs[], const struct matrix_expr *e,
+ matrix_proto_IDENT *f)
+{
+ static const struct matrix_function_properties p1 = {
+ .name = "IDENT",
+ .constraints = "ai>=0"
+ };
+ static const struct matrix_function_properties p2 = {
+ .name = "IDENT",
+ .constraints = "ai>=0 bi>=0"
+ };
+ const struct matrix_function_properties *props = e->n_subs == 1 ? &p1 : &p2;
+
+ assert (e->n_subs <= 2);
+
+ double d[2];
+ return (to_scalar_args (props->name, subs, e, d)
+ && check_constraints (props, subs, e)
+ ? f (d[0], d[e->n_subs - 1])
+ : NULL);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate (const struct matrix_expr *e)
+{
+ if (e->op == MOP_NUMBER)
+ {
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, e->number);
+ return m;
+ }
+ else if (e->op == MOP_VARIABLE)
+ {
+ const gsl_matrix *src = e->variable->value;
+ if (!src)
+ {
+ msg_at (SE, e->location,
+ _("Uninitialized variable %s used in expression."),
+ e->variable->name);
+ return NULL;
+ }
+
+ gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
+ gsl_matrix_memcpy (dst, src);
+ return dst;
+ }
+ else if (e->op == MOP_EOF)
+ {
+ struct dfm_reader *reader = read_file_open (e->eof);
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
+ return m;
+ }
+
+ enum { N_LOCAL = 3 };
+ gsl_matrix *local_subs[N_LOCAL];
+ gsl_matrix **subs = (e->n_subs < N_LOCAL
+ ? local_subs
+ : xmalloc (e->n_subs * sizeof *subs));
+
+ for (size_t i = 0; i < e->n_subs; i++)
+ {
+ subs[i] = matrix_expr_evaluate (e->subs[i]);
+ if (!subs[i])
+ {
+ for (size_t j = 0; j < i; j++)
+ gsl_matrix_free (subs[j]);
+ if (subs != local_subs)
+ free (subs);
+ return NULL;
+ }
+ }
+
+ gsl_matrix *result = NULL;
+ switch (e->op)
+ {
+#define F(ENUM, STRING, PROTO, CONSTRAINTS) \
+ case MOP_F_##ENUM: \
+ { \
+ static const struct matrix_function_properties props = { \
+ .name = STRING, \
+ .constraints = CONSTRAINTS, \
+ }; \
+ result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
+ matrix_eval_##ENUM); \
+ } \
+ break;
+ MATRIX_FUNCTIONS
+#undef F
+
+ case MOP_NEGATE:
+ gsl_matrix_scale (subs[0], -1.0);
+ result = subs[0];
+ break;
+
+ case MOP_ADD_ELEMS:
+ case MOP_SUB_ELEMS:
+ case MOP_MUL_ELEMS:
+ case MOP_DIV_ELEMS:
+ case MOP_EXP_ELEMS:
+ case MOP_GT:
+ case MOP_GE:
+ case MOP_LT:
+ case MOP_LE:
+ case MOP_EQ:
+ case MOP_NE:
+ case MOP_AND:
+ case MOP_OR:
+ case MOP_XOR:
+ result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
+ break;
+
+ case MOP_NOT:
+ result = matrix_expr_evaluate_not (subs[0]);
+ break;
+
+ case MOP_SEQ:
+ result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
+ break;
+
+ case MOP_SEQ_BY:
+ result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
+ break;
+
+ case MOP_MUL_MAT:
+ result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
+ break;
+
+ case MOP_EXP_MAT:
+ result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
+ break;
+
+ case MOP_PASTE_HORZ:
+ result = matrix_expr_evaluate_paste_horz (e, subs[0], subs[1]);
+ break;
+
+ case MOP_PASTE_VERT:
+ result = matrix_expr_evaluate_paste_vert (e, subs[0], subs[1]);
+ break;
+
+ case MOP_EMPTY:
+ result = gsl_matrix_alloc (0, 0);
+ break;
+
+ case MOP_VEC_INDEX:
+ result = matrix_expr_evaluate_vec_index (e, subs[0], subs[1]);
+ break;
+
+ case MOP_VEC_ALL:
+ result = matrix_expr_evaluate_vec_all (e, subs[0]);
+ break;
+
+ case MOP_MAT_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0],
+ subs[1], e->subs[1],
+ subs[2], e->subs[2]);
+ break;
+
+ case MOP_ROW_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0],
+ subs[1], e->subs[1],
+ NULL, NULL);
+ break;
+
+ case MOP_COL_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0],
+ NULL, NULL,
+ subs[1], e->subs[1]);
+ break;
+
+ case MOP_NUMBER:
+ case MOP_VARIABLE:
+ case MOP_EOF:
+ NOT_REACHED ();
+ }
+
+ for (size_t i = 0; i < e->n_subs; i++)
+ if (subs[i] != result)
+ gsl_matrix_free (subs[i]);
+ if (subs != local_subs)
+ free (subs);
+ return result;
+}
+
+static bool
+matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
+ double *d)
+{
+ gsl_matrix *m = matrix_expr_evaluate (e);
+ if (!m)
+ return false;
+
+ if (!is_scalar (m))
+ {
+ msg_at (SE, matrix_expr_location (e),
+ _("Expression for %s must evaluate to scalar, "
+ "not a %zu×%zu matrix."),
+ context, m->size1, m->size2);
+ gsl_matrix_free (m);
+ return false;
+ }
+
+ *d = to_scalar (m);
+ gsl_matrix_free (m);
+ return true;
+}
+
+static bool
+matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
+ long int *integer)
+{
+ double d;
+ if (!matrix_expr_evaluate_scalar (e, context, &d))
+ return false;
+
+ d = trunc (d);
+ if (d < LONG_MIN || d > LONG_MAX)
+ {
+ msg_at (SE, matrix_expr_location (e),
+ _("Expression for %s is outside the integer range."), context);
+ return false;
+ }
+ *integer = d;
+ return true;
+}
+\f
+/* Matrix lvalues.
+
+ An lvalue is an expression that can appear on the left side of a COMPUTE
+ command and in other contexts that assign values.
+
+ An lvalue is parsed once, with matrix_lvalue_parse(). It can then be
+ evaluated (with matrix_lvalue_evaluate()) and assigned (with
+ matrix_lvalue_assign()).
+
+ There are three kinds of lvalues:
+
+ - A variable name. A variable used as an lvalue need not be initialized,
+ since the assignment will initialize.
+
+ - A subvector, e.g. "var(index0)". The variable must be initialized and
+ must have the form of a vector (it must have 1 column or 1 row).
+
+ - A submatrix, e.g. "var(index0, index1)". The variable must be
+ initialized. */
+struct matrix_lvalue
+ {
+ struct matrix_var *var; /* Destination variable. */
+ struct matrix_expr *indexes[2]; /* Index expressions, if any. */
+ size_t n_indexes; /* Number of indexes. */
+
+ struct msg_location *var_location; /* Variable name. */
+ struct msg_location *full_location; /* Variable name plus indexing. */
+ struct msg_location *index_locations[2]; /* Index expressions. */
+ };
+
+/* Frees LVALUE. */
+static void
+matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
+{
+ if (lvalue)
+ {
+ msg_location_destroy (lvalue->var_location);
+ msg_location_destroy (lvalue->full_location);
+ for (size_t i = 0; i < lvalue->n_indexes; i++)
+ {
+ matrix_expr_destroy (lvalue->indexes[i]);
+ msg_location_destroy (lvalue->index_locations[i]);
+ }
+ free (lvalue);
+ }
+}
+
+/* Parses and returns an lvalue at the current position in S's lexer. Returns
+ null on parse failure. On success, the caller must eventually free the
+ lvalue. */
+static struct matrix_lvalue *
+matrix_lvalue_parse (struct matrix_state *s)
+{
+ if (!lex_force_id (s->lexer))
+ return NULL;
+
+ struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
+ int start_ofs = lex_ofs (s->lexer);
+ lvalue->var_location = lex_get_location (s->lexer, 0, 0);
+ lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (lex_next_token (s->lexer, 1) == T_LPAREN)
+ {
+ if (!lvalue->var)
+ {
+ lex_error (s->lexer, _("Undefined variable %s."),
+ lex_tokcstr (s->lexer));
+ goto error;
+ }
+
+ lex_get_n (s->lexer, 2);
+
+ if (!matrix_parse_index_expr (s, &lvalue->indexes[0],
+ &lvalue->index_locations[0]))
+ goto error;
+ lvalue->n_indexes++;
+
+ if (lex_match (s->lexer, T_COMMA))
+ {
+ if (!matrix_parse_index_expr (s, &lvalue->indexes[1],
+ &lvalue->index_locations[1]))
+ goto error;
+ lvalue->n_indexes++;
+ }
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ lvalue->full_location = lex_ofs_location (s->lexer, start_ofs,
+ lex_ofs (s->lexer) - 1);
+ }
+ else
+ {
+ if (!lvalue->var)
+ lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
+ lex_get (s->lexer);
+ }
+ return lvalue;
+
+error:
+ matrix_lvalue_destroy (lvalue);
+ return NULL;
+}
+
+static bool
+matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
+ enum index_type index_type, size_t other_size,
+ struct index_vector *iv)
+{
+ gsl_matrix *m;
+ if (e)
+ {
+ m = matrix_expr_evaluate (e);
+ if (!m)
+ return false;
+ }
+ else
+ m = NULL;
+
+ bool ok = matrix_normalize_index_vector (m, e, size, index_type,
+ other_size, iv);
+ gsl_matrix_free (m);
+ return ok;
+}
+
+/* Evaluates the indexes in LVALUE into IV0 and IV1, owned by the caller.
+ Returns true if successful, false if evaluating the expressions failed or if
+ LVALUE otherwise can't be used for an assignment.
+
+ On success, the caller retains ownership of the index vectors, which are
+ suitable for passing to matrix_lvalue_assign(). If not used for that
+ purpose then they need to eventually be freed (with
+ index_vector_uninit()). */
+static bool
+matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
+ struct index_vector *iv0,
+ struct index_vector *iv1)
+{
+ *iv0 = INDEX_VECTOR_INIT;
+ *iv1 = INDEX_VECTOR_INIT;
+ if (!lvalue->n_indexes)
+ return true;
+
+ /* Validate destination matrix exists and has the right shape. */
+ gsl_matrix *dm = lvalue->var->value;
+ if (!dm)
+ {
+ msg_at (SE, lvalue->var_location,
+ _("Undefined variable %s."), lvalue->var->name);
+ return false;
+ }
+ else if (dm->size1 == 0 || dm->size2 == 0)
+ {
+ msg_at (SE, lvalue->full_location, _("Cannot index %zu×%zu matrix %s."),
+ dm->size1, dm->size2, lvalue->var->name);
+ return false;
+ }
+ else if (lvalue->n_indexes == 1)
+ {
+ if (!is_vector (dm))
+ {
+ msg_at (SE, lvalue->full_location,
+ _("Can't use vector indexing on %zu×%zu matrix %s."),
+ dm->size1, dm->size2, lvalue->var->name);
+ return false;
+ }
+ return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
+ MAX (dm->size1, dm->size2),
+ IV_VECTOR, 0, iv0);
+ }
+ else
+ {
+ assert (lvalue->n_indexes == 2);
+ if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
+ IV_ROW, dm->size2, iv0))
+ return false;
+
+ if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
+ IV_COLUMN, dm->size1, iv1))
+ {
+ index_vector_uninit (iv0);
+ return false;
+ }
+ return true;
+ }
+}
+
+static bool
+matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
+ struct index_vector *iv,
+ gsl_matrix *sm, const struct msg_location *lsm)
+{
+ /* Convert source matrix 'sm' to source vector 'sv'. */
+ if (!is_vector (sm))
+ {
+ msg_at (SE, lvalue->full_location,
+ _("Only an %zu-element vector may be assigned to this "
+ "%zu-element subvector of %s."),
+ iv->n, iv->n, lvalue->var->name);
+ msg_at (SE, lsm,
+ _("The source is an %zu×%zu matrix."),
+ sm->size1, sm->size2);
+ return false;
+ }
+ gsl_vector sv = to_vector (sm);
+ if (iv->n != sv.size)
+ {
+ msg_at (SE, lvalue->full_location,
+ _("Only an %zu-element vector may be assigned to this "
+ "%zu-element subvector of %s."),
+ iv->n, iv->n, lvalue->var->name);
+ msg_at (SE, lsm, ngettext ("The source vector has %zu element.",
+ "The source vector has %zu elements.",
+ sv.size),
+ sv.size);
+ return false;
+ }
+
+ /* Assign elements. */
+ gsl_vector dv = to_vector (lvalue->var->value);
+ for (size_t x = 0; x < iv->n; x++)
+ gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
+ return true;
+}
+
+static bool
+matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
+ struct index_vector *iv0,
+ struct index_vector *iv1,
+ gsl_matrix *sm, const struct msg_location *lsm)
+{
+ gsl_matrix *dm = lvalue->var->value;
+
+ /* Convert source matrix 'sm' to source vector 'sv'. */
+ bool wrong_rows = iv0->n != sm->size1;
+ bool wrong_columns = iv1->n != sm->size2;
+ if (wrong_rows || wrong_columns)
+ {
+ if (wrong_rows && wrong_columns)
+ msg_at (SE, lvalue->full_location,
+ _("Numbers of indexes for assigning to %s differ from the "
+ "size of the source matrix."),
+ lvalue->var->name);
+ else if (wrong_rows)
+ msg_at (SE, lvalue->full_location,
+ _("Number of row indexes for assigning to %s differs from "
+ "number of rows in source matrix."),
+ lvalue->var->name);
+ else
+ msg_at (SE, lvalue->full_location,
+ _("Number of column indexes for assigning to %s differs from "
+ "number of columns in source matrix."),
+ lvalue->var->name);
+
+ if (wrong_rows)
+ {
+ if (lvalue->indexes[0])
+ msg_at (SN, lvalue->index_locations[0],
+ ngettext ("There is %zu row index.",
+ "There are %zu row indexes.",
+ iv0->n),
+ iv0->n);
+ else
+ msg_at (SN, lvalue->index_locations[0],
+ ngettext ("Destination matrix %s has %zu row.",
+ "Destination matrix %s has %zu rows.",
+ iv0->n),
+ lvalue->var->name, iv0->n);
+ }
+
+ if (wrong_columns)
+ {
+ if (lvalue->indexes[1])
+ msg_at (SN, lvalue->index_locations[1],
+ ngettext ("There is %zu column index.",
+ "There are %zu column indexes.",
+ iv1->n),
+ iv1->n);
+ else
+ msg_at (SN, lvalue->index_locations[1],
+ ngettext ("Destination matrix %s has %zu column.",
+ "Destination matrix %s has %zu columns.",
+ iv1->n),
+ lvalue->var->name, iv1->n);
+ }
+
+ msg_at (SN, lsm, _("The source matrix is %zu×%zu."),
+ sm->size1, sm->size2);
+ return false;
+ }
+
+ /* Assign elements. */
+ for (size_t y = 0; y < iv0->n; y++)
+ {
+ size_t f0 = iv0->indexes[y];
+ for (size_t x = 0; x < iv1->n; x++)
+ {
+ size_t f1 = iv1->indexes[x];
+ gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
+ }
+ }
+ return true;
+}
+
+/* Assigns rvalue SM to LVALUE, whose index vectors IV0 and IV1 were previously
+ obtained with matrix_lvalue_evaluate(). Returns true if successful, false
+ on error. Always takes ownership of M. LSM should be the source location
+ to use for errors related to SM. */
+static bool
+matrix_lvalue_assign (struct matrix_lvalue *lvalue,
+ struct index_vector *iv0, struct index_vector *iv1,
+ gsl_matrix *sm, const struct msg_location *lsm)
+{
+ if (!lvalue->n_indexes)
+ {
+ gsl_matrix_free (lvalue->var->value);
+ lvalue->var->value = sm;
+ return true;
+ }
+ else
+ {
+ bool ok = (lvalue->n_indexes == 1
+ ? matrix_lvalue_assign_vector (lvalue, iv0, sm, lsm)
+ : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm, lsm));
+ index_vector_uninit (iv0);
+ index_vector_uninit (iv1);
+ gsl_matrix_free (sm);
+ return ok;
+ }
+}
+
+/* Evaluates and then assigns SM to LVALUE. Always takes ownership of M. LSM
+ should be the source location to use for errors related to SM.*/
+static bool
+matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue,
+ gsl_matrix *sm,
+ const struct msg_location *lsm)
+{
+ struct index_vector iv0, iv1;
+ if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
+ {
+ gsl_matrix_free (sm);
+ return false;
+ }
+
+ return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm, lsm);
+}
+\f
+/* Matrix command data structure. */
+
+/* An array of matrix commands. */
+struct matrix_commands
+ {
+ struct matrix_command **commands;
+ size_t n;
+ };
+
+static bool matrix_commands_parse (struct matrix_state *,
+ struct matrix_commands *,
+ const char *command_name,
+ const char *stop1, const char *stop2);
+static void matrix_commands_uninit (struct matrix_commands *);
+
+/* A single matrix command. */
+struct matrix_command
+ {
+ /* The type of command. */
+ enum matrix_command_type
+ {
+ MCMD_COMPUTE,
+ MCMD_PRINT,
+ MCMD_DO_IF,
+ MCMD_LOOP,
+ MCMD_BREAK,
+ MCMD_DISPLAY,
+ MCMD_RELEASE,
+ MCMD_SAVE,
+ MCMD_READ,
+ MCMD_WRITE,
+ MCMD_GET,
+ MCMD_MSAVE,
+ MCMD_MGET,
+ MCMD_EIGEN,
+ MCMD_SETDIAG,
+ MCMD_SVD,
+ }
+ type;
+
+ /* Source lines for this command. */
+ struct msg_location *location;
+
+ union
+ {
+ struct matrix_compute
+ {
+ struct matrix_lvalue *lvalue;
+ struct matrix_expr *rvalue;
+ }
+ compute;
+
+ struct matrix_print
+ {
+ struct matrix_expr *expression;
+ bool use_default_format;
+ struct fmt_spec format;
+ char *title;
+ int space; /* -1 means NEWPAGE. */
+
+ struct print_labels
+ {
+ struct string_array literals; /* CLABELS/RLABELS. */
+ struct matrix_expr *expr; /* CNAMES/RNAMES. */
+ }
+ *rlabels, *clabels;
+ }
+ print;
+
+ struct matrix_do_if
+ {
+ struct do_if_clause
+ {
+ struct matrix_expr *condition;
+ struct matrix_commands commands;
+ }
+ *clauses;
+ size_t n_clauses;
+ }
+ do_if;
+
+ struct matrix_loop
+ {
+ /* Index. */
+ struct matrix_var *var;
+ struct matrix_expr *start;
+ struct matrix_expr *end;
+ struct matrix_expr *increment;
+
+ /* Loop conditions. */
+ struct matrix_expr *top_condition;
+ struct matrix_expr *bottom_condition;
+
+ /* Commands. */
+ struct matrix_commands commands;
+ }
+ loop;
+
+ struct matrix_display
+ {
+ struct matrix_state *state;
+ }
+ display;
+
+ struct matrix_release
+ {
+ struct matrix_var **vars;
+ size_t n_vars;
+ }
+ release;
+
+ struct matrix_save
+ {
+ struct matrix_expr *expression;
+ struct save_file *sf;
+ }
+ save;
+
+ struct matrix_read
+ {
+ struct read_file *rf;
+ struct matrix_lvalue *dst;
+ struct matrix_expr *size;
+ int c1, c2;
+ enum fmt_type format;
+ int w;
+ bool symmetric;
+ bool reread;
+ }
+ read;
+
+ struct matrix_write
+ {
+ struct write_file *wf;
+ struct matrix_expr *expression;
+ int c1, c2;
+
+ /* If this is nonnull, WRITE uses this format.
+
+ If this is NULL, WRITE uses free-field format with as many
+ digits of precision as needed. */
+ struct fmt_spec *format;
+
+ bool triangular;
+ bool hold;
+ }
+ write;
+
+ struct matrix_get
+ {
+ struct lexer *lexer;
+ struct matrix_lvalue *dst;
+ struct dataset *dataset;
+ struct file_handle *file;
+ char *encoding;
+ struct var_syntax *vars;
+ size_t n_vars;
+ struct matrix_var *names;
+
+ /* Treatment of missing values. */
+ struct
+ {
+ enum
+ {
+ MGET_ERROR, /* Flag error on command. */
+ MGET_ACCEPT, /* Accept without change, user-missing only. */
+ MGET_OMIT, /* Drop this case. */
+ MGET_RECODE /* Recode to 'substitute'. */
+ }
+ treatment;
+ double substitute; /* MGET_RECODE only. */
+ }
+ user, system;
+ }
+ get;
+
+ struct matrix_msave
+ {
+ struct msave_common *common;
+ struct matrix_expr *expr;
+ const char *rowtype;
+ const struct matrix_expr *factors;
+ const struct matrix_expr *splits;
+ }
+ msave;
+
+ struct matrix_mget
+ {
+ struct matrix_state *state;
+ struct file_handle *file;
+ char *encoding;
+ struct stringi_set rowtypes;
+ }
+ mget;
+
+ struct matrix_eigen
+ {
+ struct matrix_expr *expr;
+ struct matrix_var *evec;
+ struct matrix_var *eval;
+ }
+ eigen;
+
+ struct matrix_setdiag
+ {
+ struct matrix_var *dst;
+ struct matrix_expr *expr;
+ }
+ setdiag;
+
+ struct matrix_svd
+ {
+ struct matrix_expr *expr;
+ struct matrix_var *u;
+ struct matrix_var *s;
+ struct matrix_var *v;
+ }
+ svd;
+ };
+ };
+
+static struct matrix_command *matrix_command_parse (struct matrix_state *);
+static bool matrix_command_execute (struct matrix_command *);
+static void matrix_command_destroy (struct matrix_command *);
+\f
+/* COMPUTE. */
+
+static struct matrix_command *
+matrix_compute_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_COMPUTE,
+ .compute = { .lvalue = NULL },
+ };
+
+ struct matrix_compute *compute = &cmd->compute;
+ compute->lvalue = matrix_lvalue_parse (s);
+ if (!compute->lvalue)
+ goto error;
+
+ if (!lex_force_match (s->lexer, T_EQUALS))
+ goto error;
+
+ compute->rvalue = matrix_expr_parse (s);
+ if (!compute->rvalue)
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_compute_execute (struct matrix_command *cmd)
+{
+ struct matrix_compute *compute = &cmd->compute;
+ gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
+ if (!value)
+ return;
+
+ matrix_lvalue_evaluate_and_assign (compute->lvalue, value,
+ matrix_expr_location (compute->rvalue));
+}
+\f
+/* PRINT. */
+
+static void
+print_labels_destroy (struct print_labels *labels)
+{
+ if (labels)
+ {
+ string_array_destroy (&labels->literals);
+ matrix_expr_destroy (labels->expr);
+ free (labels);
+ }
+}
+
+static void
+parse_literal_print_labels (struct matrix_state *s,
+ struct print_labels **labelsp)
+{
+ lex_match (s->lexer, T_EQUALS);
+ print_labels_destroy (*labelsp);
+ *labelsp = xzalloc (sizeof **labelsp);
+ while (lex_token (s->lexer) != T_SLASH
+ && lex_token (s->lexer) != T_ENDCMD
+ && lex_token (s->lexer) != T_STOP)
+ {
+ struct string label = DS_EMPTY_INITIALIZER;
+ while (lex_token (s->lexer) != T_COMMA
+ && lex_token (s->lexer) != T_SLASH
+ && lex_token (s->lexer) != T_ENDCMD
+ && lex_token (s->lexer) != T_STOP)
+ {
+ if (!ds_is_empty (&label))
+ ds_put_byte (&label, ' ');
+
+ if (lex_token (s->lexer) == T_STRING)
+ ds_put_cstr (&label, lex_tokcstr (s->lexer));
+ else
+ {
+ char *rep = lex_next_representation (s->lexer, 0, 0);
+ ds_put_cstr (&label, rep);
+ free (rep);
+ }
+ lex_get (s->lexer);
+ }
+ string_array_append_nocopy (&(*labelsp)->literals,
+ ds_steal_cstr (&label));
+
+ if (!lex_match (s->lexer, T_COMMA))
+ break;
+ }
+}
+
+static bool
+parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
+{
+ lex_match (s->lexer, T_EQUALS);
+ struct matrix_expr *e = matrix_parse_exp (s);
+ if (!e)
+ return false;
+
+ print_labels_destroy (*labelsp);
+ *labelsp = xzalloc (sizeof **labelsp);
+ (*labelsp)->expr = e;
+ return true;
+}
+
+static struct matrix_command *
+matrix_print_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_PRINT,
+ .print = {
+ .use_default_format = true,
+ }
+ };
+
+ if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
+ {
+ int start_ofs = lex_ofs (s->lexer);
+ cmd->print.expression = matrix_parse_exp (s);
+ if (!cmd->print.expression)
+ goto error;
+ cmd->print.title = lex_ofs_representation (s->lexer, start_ofs,
+ lex_ofs (s->lexer) - 1);
+ }
+
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "FORMAT"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!parse_format_specifier (s->lexer, &cmd->print.format))
+ goto error;
+ cmd->print.use_default_format = false;
+ }
+ else if (lex_match_id (s->lexer, "TITLE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_string (s->lexer))
+ goto error;
+ free (cmd->print.title);
+ cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "SPACE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (lex_match_id (s->lexer, "NEWPAGE"))
+ cmd->print.space = -1;
+ else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
+ {
+ cmd->print.space = lex_integer (s->lexer);
+ lex_get (s->lexer);
+ }
+ else
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "RLABELS"))
+ parse_literal_print_labels (s, &cmd->print.rlabels);
+ else if (lex_match_id (s->lexer, "CLABELS"))
+ parse_literal_print_labels (s, &cmd->print.clabels);
+ else if (lex_match_id (s->lexer, "RNAMES"))
+ {
+ if (!parse_expr_print_labels (s, &cmd->print.rlabels))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "CNAMES"))
+ {
+ if (!parse_expr_print_labels (s, &cmd->print.clabels))
+ goto error;
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
+ "RLABELS", "CLABELS", "RNAMES", "CNAMES");
+ goto error;
+ }
+
+ }
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static bool
+matrix_is_integer (const gsl_matrix *m)
+{
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double d = gsl_matrix_get (m, y, x);
+ if (d != floor (d))
+ return false;
+ }
+ return true;
+}
+
+static double
+matrix_max_magnitude (const gsl_matrix *m)
+{
+ double max = 0.0;
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double d = fabs (gsl_matrix_get (m, y, x));
+ if (d > max)
+ max = d;
+ }
+ return max;
+}
+
+static bool
+format_fits (struct fmt_spec format, double x)
+{
+ char *s = data_out (&(union value) { .f = x }, NULL,
+ &format, settings_get_fmt_settings ());
+ bool fits = *s != '*' && !strchr (s, 'E');
+ free (s);
+ return fits;
+}
+
+static struct fmt_spec
+default_format (const gsl_matrix *m, int *log_scale)
+{
+ *log_scale = 0;
+
+ double max = matrix_max_magnitude (m);
+
+ if (matrix_is_integer (m))
+ for (int w = 1; w <= 12; w++)
+ {
+ struct fmt_spec format = { .type = FMT_F, .w = w };
+ if (format_fits (format, -max))
+ return format;
+ }
+
+ if (max >= 1e9 || max <= 1e-4)
+ {
+ char s[64];
+ snprintf (s, sizeof s, "%.1e", max);
+
+ const char *e = strchr (s, 'e');
+ if (e)
+ *log_scale = atoi (e + 1);
+ }
+
+ return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
+}
+
+static char *
+trimmed_string (double d)
+{
+ struct substring s = ss_buffer ((char *) &d, sizeof d);
+ ss_rtrim (&s, ss_cstr (" "));
+ return ss_xstrdup (s);
+}
+
+static struct string_array *
+print_labels_get (const struct print_labels *labels, size_t n,
+ const char *prefix, bool truncate)
+{
+ if (!labels)
+ return NULL;
+
+ struct string_array *out = xzalloc (sizeof *out);
+ if (labels->literals.n)
+ string_array_clone (out, &labels->literals);
+ else if (labels->expr)
+ {
+ gsl_matrix *m = matrix_expr_evaluate (labels->expr);
+ if (m && is_vector (m))
+ {
+ gsl_vector v = to_vector (m);
+ for (size_t i = 0; i < v.size; i++)
+ string_array_append_nocopy (out, trimmed_string (
+ gsl_vector_get (&v, i)));
+ }
+ gsl_matrix_free (m);
+ }
+
+ while (out->n < n)
+ {
+ if (labels->expr)
+ string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
+ out->n + 1));
+ else
+ string_array_append (out, "");
+ }
+ while (out->n > n)
+ string_array_delete (out, out->n - 1);
+
+ if (truncate)
+ for (size_t i = 0; i < out->n; i++)
+ {
+ char *s = out->strings[i];
+ s[strnlen (s, 8)] = '\0';
+ }
+
+ return out;
+}
+
+static void
+matrix_print_space (int space)
+{
+ if (space < 0)
+ output_item_submit (page_break_item_create ());
+ for (int i = 0; i < space; i++)
+ output_log ("%s", "");
+}
+
+static void
+matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
+ struct fmt_spec format, int log_scale)
+{
+ matrix_print_space (print->space);
+ if (print->title)
+ output_log ("%s", print->title);
+ if (log_scale != 0)
+ output_log (" 10 ** %d X", log_scale);
+
+ struct string_array *clabels = print_labels_get (print->clabels,
+ m->size2, "col", true);
+ if (clabels && format.w < 8)
+ format.w = 8;
+ struct string_array *rlabels = print_labels_get (print->rlabels,
+ m->size1, "row", true);
+
+ if (clabels)
+ {
+ struct string line = DS_EMPTY_INITIALIZER;
+ if (rlabels)
+ ds_put_byte_multiple (&line, ' ', 8);
+ for (size_t x = 0; x < m->size2; x++)
+ ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
+ output_log_nocopy (ds_steal_cstr (&line));
+ }
+
+ double scale = pow (10.0, log_scale);
+ bool numeric = fmt_is_numeric (format.type);
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ struct string line = DS_EMPTY_INITIALIZER;
+ if (rlabels)
+ ds_put_format (&line, "%-8s", rlabels->strings[y]);
+
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double f = gsl_matrix_get (m, y, x);
+ char *s = (numeric
+ ? data_out (&(union value) { .f = f / scale}, NULL,
+ &format, settings_get_fmt_settings ())
+ : trimmed_string (f));
+ ds_put_format (&line, " %s", s);
+ free (s);
+ }
+ output_log_nocopy (ds_steal_cstr (&line));
+ }
+
+ string_array_destroy (rlabels);
+ free (rlabels);
+ string_array_destroy (clabels);
+ free (clabels);
+}
+
+static void
+create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
+ const struct print_labels *print_labels, size_t n,
+ const char *name, const char *prefix)
+{
+ struct string_array *labels = print_labels_get (print_labels, n, prefix,
+ false);
+ struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
+ for (size_t i = 0; i < n; i++)
+ pivot_category_create_leaf (
+ d->root, (labels
+ ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
+ : pivot_value_new_integer (i + 1)));
+ if (!labels)
+ d->hide_all_labels = true;
+ string_array_destroy (labels);
+ free (labels);
+}
+
+static void
+matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
+ struct fmt_spec format, int log_scale)
+{
+ struct pivot_table *table = pivot_table_create__ (
+ pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
+ "Matrix Print");
+
+ create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
+ N_("Rows"), "row");
+ create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
+ N_("Columns"), "col");
+
+ struct pivot_footnote *footnote = NULL;
+ if (log_scale != 0)
+ {
+ char *s = xasprintf ("× 10**%d", log_scale);
+ footnote = pivot_table_create_footnote (
+ table, pivot_value_new_user_text_nocopy (s));
+ }
+
+ double scale = pow (10.0, log_scale);
+ bool numeric = fmt_is_numeric (format.type);
+ for (size_t y = 0; y < m->size1; y++)
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double f = gsl_matrix_get (m, y, x);
+ struct pivot_value *v;
+ if (numeric)
+ {
+ v = pivot_value_new_number (f / scale);
+ v->numeric.format = format;
+ }
+ else
+ v = pivot_value_new_user_text_nocopy (trimmed_string (f));
+ if (footnote)
+ pivot_value_add_footnote (v, footnote);
+ pivot_table_put2 (table, y, x, v);
+ }
+
+ pivot_table_submit (table);
+}
+
+static void
+matrix_print_execute (const struct matrix_print *print)
+{
+ if (print->expression)
+ {
+ gsl_matrix *m = matrix_expr_evaluate (print->expression);
+ if (!m)
+ return;
+
+ int log_scale = 0;
+ struct fmt_spec format = (print->use_default_format
+ ? default_format (m, &log_scale)
+ : print->format);
+
+ if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
+ matrix_print_text (print, m, format, log_scale);
+ else
+ matrix_print_tables (print, m, format, log_scale);
+
+ gsl_matrix_free (m);
+ }
+ else
+ {
+ matrix_print_space (print->space);
+ if (print->title)
+ output_log ("%s", print->title);
+ }
+}
+\f
+/* DO IF. */
+
+static bool
+matrix_do_if_clause_parse (struct matrix_state *s,
+ struct matrix_do_if *ifc,
+ bool parse_condition,
+ size_t *allocated_clauses)
+{
+ if (ifc->n_clauses >= *allocated_clauses)
+ ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
+ sizeof *ifc->clauses);
+ struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
+ *c = (struct do_if_clause) { .condition = NULL };
+
+ if (parse_condition)
+ {
+ c->condition = matrix_expr_parse (s);
+ if (!c->condition)
+ return false;
+ }
+
+ return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
+}
+
+static struct matrix_command *
+matrix_do_if_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_DO_IF,
+ .do_if = { .n_clauses = 0 }
+ };
+
+ size_t allocated_clauses = 0;
+ do
+ {
+ if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
+ goto error;
+ }
+ while (lex_match_phrase (s->lexer, "ELSE IF"));
+
+ if (lex_match_id (s->lexer, "ELSE")
+ && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
+ goto error;
+
+ if (!lex_match_phrase (s->lexer, "END IF"))
+ NOT_REACHED ();
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static bool
+matrix_do_if_execute (struct matrix_do_if *cmd)
+{
+ for (size_t i = 0; i < cmd->n_clauses; i++)
+ {
+ struct do_if_clause *c = &cmd->clauses[i];
+ if (c->condition)
+ {
+ double d;
+ if (!matrix_expr_evaluate_scalar (c->condition,
+ i ? "ELSE IF" : "DO IF",
+ &d) || d <= 0)
+ continue;
+ }
+
+ for (size_t j = 0; j < c->commands.n; j++)
+ if (!matrix_command_execute (c->commands.commands[j]))
+ return false;
+ return true;
+ }
+ return true;
+}
+\f
+/* LOOP. */
+
+static struct matrix_command *
+matrix_loop_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
+
+ struct matrix_loop *loop = &cmd->loop;
+ if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
+ {
+ struct substring name = lex_tokss (s->lexer);
+ loop->var = matrix_var_lookup (s, name);
+ if (!loop->var)
+ loop->var = matrix_var_create (s, name);
+
+ lex_get (s->lexer);
+ lex_get (s->lexer);
+
+ loop->start = matrix_expr_parse (s);
+ if (!loop->start || !lex_force_match (s->lexer, T_TO))
+ goto error;
+
+ loop->end = matrix_expr_parse (s);
+ if (!loop->end)
+ goto error;
+
+ if (lex_match (s->lexer, T_BY))
+ {
+ loop->increment = matrix_expr_parse (s);
+ if (!loop->increment)
+ goto error;
+ }
+ }
+
+ if (lex_match_id (s->lexer, "IF"))
+ {
+ loop->top_condition = matrix_expr_parse (s);
+ if (!loop->top_condition)
+ goto error;
+ }
+
+ bool was_in_loop = s->in_loop;
+ s->in_loop = true;
+ bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
+ "END LOOP", NULL);
+ s->in_loop = was_in_loop;
+ if (!ok)
+ goto error;
+
+ if (!lex_match_phrase (s->lexer, "END LOOP"))
+ NOT_REACHED ();
+
+ if (lex_match_id (s->lexer, "IF"))
+ {
+ loop->bottom_condition = matrix_expr_parse (s);
+ if (!loop->bottom_condition)
+ goto error;
+ }
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_loop_execute (struct matrix_loop *cmd)
+{
+ long int value = 0;
+ long int end = 0;
+ long int increment = 1;
+ if (cmd->var)
+ {
+ if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
+ || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
+ || (cmd->increment
+ && !matrix_expr_evaluate_integer (cmd->increment, "BY",
+ &increment)))
+ return;
+
+ if (increment > 0 ? value > end
+ : increment < 0 ? value < end
+ : true)
+ return;
+ }
+
+ int mxloops = settings_get_mxloops ();
+ for (int i = 0; i < mxloops; i++)
+ {
+ if (cmd->var)
+ {
+ if (cmd->var->value
+ && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
+ {
+ gsl_matrix_free (cmd->var->value);
+ cmd->var->value = NULL;
+ }
+ if (!cmd->var->value)
+ cmd->var->value = gsl_matrix_alloc (1, 1);
+
+ gsl_matrix_set (cmd->var->value, 0, 0, value);
+ }
+
+ if (cmd->top_condition)
+ {
+ double d;
+ if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
+ &d) || d <= 0)
+ return;
+ }
+
+ for (size_t j = 0; j < cmd->commands.n; j++)
+ if (!matrix_command_execute (cmd->commands.commands[j]))
+ return;
+
+ if (cmd->bottom_condition)
+ {
+ double d;
+ if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
+ "END LOOP IF", &d) || d > 0)
+ return;
+ }
+
+ if (cmd->var)
+ {
+ value += increment;
+ if (increment > 0 ? value > end : value < end)
+ return;
+ }
+ }
+}
+\f
+/* BREAK. */
+
+static struct matrix_command *
+matrix_break_parse (struct matrix_state *s)
+{
+ if (!s->in_loop)
+ {
+ lex_next_error (s->lexer, -1, -1, _("BREAK not inside LOOP."));
+ return NULL;
+ }
+
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) { .type = MCMD_BREAK };
+ return cmd;
+}
+\f
+/* DISPLAY. */
+
+static struct matrix_command *
+matrix_display_parse (struct matrix_state *s)
+{
+ while (lex_token (s->lexer) != T_ENDCMD)
+ {
+ if (!lex_match_id (s->lexer, "DICTIONARY")
+ && !lex_match_id (s->lexer, "STATUS"))
+ {
+ lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
+ return NULL;
+ }
+ }
+
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
+ return cmd;
+}
+
+static int
+compare_matrix_var_pointers (const void *a_, const void *b_)
+{
+ const struct matrix_var *const *ap = a_;
+ const struct matrix_var *const *bp = b_;
+ const struct matrix_var *a = *ap;
+ const struct matrix_var *b = *bp;
+ return strcmp (a->name, b->name);
+}
+
+static void
+matrix_display_execute (struct matrix_display *cmd)
+{
+ const struct matrix_state *s = cmd->state;
+
+ struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
+ struct pivot_dimension *attr_dimension
+ = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
+ pivot_category_create_group (attr_dimension->root, N_("Dimension"),
+ N_("Rows"), N_("Columns"));
+ pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
+
+ const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
+ size_t n_vars = 0;
+
+ const struct matrix_var *var;
+ HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
+ if (var->value)
+ vars[n_vars++] = var;
+ qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
+
+ struct pivot_dimension *rows = pivot_dimension_create (
+ table, PIVOT_AXIS_ROW, N_("Variable"));
+ for (size_t i = 0; i < n_vars; i++)
+ {
+ const struct matrix_var *var = vars[i];
+ pivot_category_create_leaf (
+ rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
+
+ size_t r = var->value->size1;
+ size_t c = var->value->size2;
+ double values[] = { r, c, r * c * sizeof (double) / 1024 };
+ for (size_t j = 0; j < sizeof values / sizeof *values; j++)
+ pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
+ }
+ free (vars);
+ pivot_table_submit (table);
+}
+\f
+/* RELEASE. */
+
+static struct matrix_command *
+matrix_release_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
+
+ size_t allocated_vars = 0;
+ while (lex_token (s->lexer) == T_ID)
+ {
+ struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (var)
+ {
+ if (cmd->release.n_vars >= allocated_vars)
+ cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
+ sizeof *cmd->release.vars);
+ cmd->release.vars[cmd->release.n_vars++] = var;
+ }
+ else
+ lex_error (s->lexer, _("Syntax error expecting variable name."));
+ lex_get (s->lexer);
+
+ if (!lex_match (s->lexer, T_COMMA))
+ break;
+ }
+
+ return cmd;
+}
+
+static void
+matrix_release_execute (struct matrix_release *cmd)
+{
+ for (size_t i = 0; i < cmd->n_vars; i++)
+ {
+ struct matrix_var *var = cmd->vars[i];
+ gsl_matrix_free (var->value);
+ var->value = NULL;
+ }
+}
+\f
+/* SAVE. */
+
+static struct save_file *
+save_file_create (struct matrix_state *s, struct file_handle *fh,
+ struct string_array *variables,
+ struct matrix_expr *names,
+ struct stringi_set *strings)
+{
+ for (size_t i = 0; i < s->n_save_files; i++)
+ {
+ struct save_file *sf = s->save_files[i];
+ if (fh_equal (sf->file, fh))
+ {
+ fh_unref (fh);
+
+ string_array_destroy (variables);
+ matrix_expr_destroy (names);
+ stringi_set_destroy (strings);
+
+ return sf;
+ }
+ }
+
+ struct save_file *sf = xmalloc (sizeof *sf);
+ *sf = (struct save_file) {
+ .file = fh,
+ .dataset = s->dataset,
+ .variables = *variables,
+ .names = names,
+ .strings = STRINGI_SET_INITIALIZER (sf->strings),
+ };
+
+ stringi_set_swap (&sf->strings, strings);
+ stringi_set_destroy (strings);
+
+ s->save_files = xrealloc (s->save_files,
+ (s->n_save_files + 1) * sizeof *s->save_files);
+ s->save_files[s->n_save_files++] = sf;
+ return sf;
+}
+
+static struct casewriter *
+save_file_open (struct save_file *sf, gsl_matrix *m,
+ const struct msg_location *save_location)
+{
+ if (sf->writer || sf->error)
+ {
+ if (sf->writer)
+ {
+ size_t n_variables = caseproto_get_n_widths (
+ casewriter_get_proto (sf->writer));
+ if (m->size2 != n_variables)
+ {
+ const char *file_name = (sf->file == fh_inline_file () ? "*"
+ : fh_get_name (sf->file));
+ msg_at (SE, save_location,
+ _("Cannot save %zu×%zu matrix to %s because the "
+ "first SAVE to %s in this matrix program wrote a "
+ "%zu-column matrix."),
+ m->size1, m->size2, file_name, file_name, n_variables);
+ msg_at (SE, sf->location,
+ _("This is the location of the first SAVE to %s."),
+ file_name);
+ return NULL;
+ }
+ }
+ return sf->writer;
+ }
+
+ bool ok = true;
+ struct dictionary *dict = dict_create (get_default_encoding ());
+
+ /* Fill 'names' with user-specified names if there were any, either from
+ sf->variables or sf->names. */
+ struct string_array names = { .n = 0 };
+ if (sf->variables.n)
+ string_array_clone (&names, &sf->variables);
+ else if (sf->names)
+ {
+ gsl_matrix *nm = matrix_expr_evaluate (sf->names);
+ if (nm && is_vector (nm))
+ {
+ gsl_vector nv = to_vector (nm);
+ for (size_t i = 0; i < nv.size; i++)
+ {
+ char *name = trimmed_string (gsl_vector_get (&nv, i));
+ char *error = dict_id_is_valid__ (dict, name);
+ if (!error)
+ string_array_append_nocopy (&names, name);
+ else
+ {
+ msg_at (SE, save_location, "%s", error);
+ free (error);
+ ok = false;
+ }
+ }
+ }
+ gsl_matrix_free (nm);
+ }
+
+ struct stringi_set strings;
+ stringi_set_clone (&strings, &sf->strings);
+
+ for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
+ {
+ char tmp_name[64];
+ const char *name;
+ if (i < names.n)
+ name = names.strings[i];
+ else
+ {
+ snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
+ name = tmp_name;
+ }
+
+ int width = stringi_set_delete (&strings, name) ? 8 : 0;
+ struct variable *var = dict_create_var (dict, name, width);
+ if (!var)
+ {
+ msg_at (ME, save_location,
+ _("Duplicate variable name %s in SAVE statement."), name);
+ ok = false;
+ }
+ }
+
+ if (!stringi_set_is_empty (&strings))
+ {
+ size_t count = stringi_set_count (&strings);
+ const char *example = stringi_set_node_get_string (
+ stringi_set_first (&strings));
+
+ if (count == 1)
+ msg_at (ME, save_location,
+ _("The SAVE command STRINGS subcommand specifies an unknown "
+ "variable %s."), example);
+ else
+ msg_at (ME, save_location,
+ ngettext ("The SAVE command STRINGS subcommand specifies %zu "
+ "unknown variable: %s.",
+ "The SAVE command STRINGS subcommand specifies %zu "
+ "unknown variables, including %s.",
+ count),
+ count, example);
+ ok = false;
+ }
+ stringi_set_destroy (&strings);
+ string_array_destroy (&names);
+
+ if (!ok)
+ {
+ dict_unref (dict);
+ sf->error = true;
+ return NULL;
+ }
+
+ if (sf->file == fh_inline_file ())
+ sf->writer = autopaging_writer_create (dict_get_proto (dict));
+ else
+ sf->writer = any_writer_open (sf->file, dict);
+ if (sf->writer)
+ {
+ sf->dict = dict;
+ sf->location = msg_location_dup (save_location);
+ }
+ else
+ {
+ dict_unref (dict);
+ sf->error = true;
+ }
+ return sf->writer;
+}
+
+static void
+save_file_destroy (struct save_file *sf)
+{
+ if (sf)
+ {
+ if (sf->file == fh_inline_file () && sf->writer && sf->dict)
+ {
+ dataset_set_dict (sf->dataset, sf->dict);
+ dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
+ }
+ else
+ {
+ casewriter_destroy (sf->writer);
+ dict_unref (sf->dict);
+ }
+ fh_unref (sf->file);
+ string_array_destroy (&sf->variables);
+ matrix_expr_destroy (sf->names);
+ stringi_set_destroy (&sf->strings);
+ msg_location_destroy (sf->location);
+ free (sf);
+ }
+}
+
+static struct matrix_command *
+matrix_save_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_SAVE,
+ .save = { .expression = NULL },
+ };
+
+ struct file_handle *fh = NULL;
+ struct matrix_save *save = &cmd->save;
+
+ struct string_array variables = STRING_ARRAY_INITIALIZER;
+ struct matrix_expr *names = NULL;
+ struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
+
+ save->expression = matrix_parse_exp (s);
+ if (!save->expression)
+ goto error;
+
+ int names_start = 0;
+ int names_end = 0;
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "OUTFILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (fh);
+ fh = (lex_match (s->lexer, T_ASTERISK)
+ ? fh_inline_file ()
+ : fh_parse (s->lexer, FH_REF_FILE, s->session));
+ if (!fh)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "VARIABLES"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ char **names;
+ size_t n;
+ struct dictionary *d = dict_create (get_default_encoding ());
+ bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
+ PV_NO_SCRATCH | PV_NO_DUPLICATE);
+ dict_unref (d);
+ if (!ok)
+ goto error;
+
+ string_array_clear (&variables);
+ variables = (struct string_array) {
+ .strings = names,
+ .n = n,
+ .allocated = n,
+ };
+ }
+ else if (lex_match_id (s->lexer, "NAMES"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ matrix_expr_destroy (names);
+ names_start = lex_ofs (s->lexer);
+ names = matrix_parse_exp (s);
+ names_end = lex_ofs (s->lexer) - 1;
+ if (!names)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "STRINGS"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ while (lex_token (s->lexer) == T_ID)
+ {
+ stringi_set_insert (&strings, lex_tokcstr (s->lexer));
+ lex_get (s->lexer);
+ if (!lex_match (s->lexer, T_COMMA))
+ break;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
+ "STRINGS");
+ goto error;
+ }
+ }
+
+ if (!fh)
+ {
+ if (s->prev_save_file)
+ fh = fh_ref (s->prev_save_file);
+ else
+ {
+ lex_sbc_missing (s->lexer, "OUTFILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_save_file);
+ s->prev_save_file = fh_ref (fh);
+
+ if (variables.n && names)
+ {
+ lex_ofs_msg (s->lexer, SW, names_start, names_end,
+ _("Ignoring NAMES because VARIABLES was also specified."));
+ matrix_expr_destroy (names);
+ names = NULL;
+ }
+
+ save->sf = save_file_create (s, fh, &variables, names, &strings);
+ return cmd;
+
+error:
+ string_array_destroy (&variables);
+ matrix_expr_destroy (names);
+ stringi_set_destroy (&strings);
+ fh_unref (fh);
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_save_execute (const struct matrix_command *cmd)
+{
+ const struct matrix_save *save = &cmd->save;
+
+ gsl_matrix *m = matrix_expr_evaluate (save->expression);
+ if (!m)
+ return;
+
+ struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
+ if (!writer)
+ {
+ gsl_matrix_free (m);
+ return;
+ }
+
+ const struct caseproto *proto = casewriter_get_proto (writer);
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ struct ccase *c = case_create (proto);
+ for (size_t x = 0; x < m->size2; x++)
+ {
+ double d = gsl_matrix_get (m, y, x);
+ int width = caseproto_get_width (proto, x);
+ union value *value = case_data_rw_idx (c, x);
+ if (width == 0)
+ value->f = d;
+ else
+ memcpy (value->s, &d, width);
+ }
+ casewriter_write (writer, c);
+ }
+ gsl_matrix_free (m);
+}
+\f
+/* READ. */
+
+static struct read_file *
+read_file_create (struct matrix_state *s, struct file_handle *fh)
+{
+ for (size_t i = 0; i < s->n_read_files; i++)
+ {
+ struct read_file *rf = s->read_files[i];
+ if (rf->file == fh)
+ {
+ fh_unref (fh);
+ return rf;
+ }
+ }
+
+ struct read_file *rf = xmalloc (sizeof *rf);
+ *rf = (struct read_file) { .file = fh };
+
+ s->read_files = xrealloc (s->read_files,
+ (s->n_read_files + 1) * sizeof *s->read_files);
+ s->read_files[s->n_read_files++] = rf;
+ return rf;
+}
+
+static struct dfm_reader *
+read_file_open (struct read_file *rf)
+{
+ if (!rf->reader)
+ rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
+ return rf->reader;
+}
+
+static void
+read_file_destroy (struct read_file *rf)
+{
+ if (rf)
+ {
+ fh_unref (rf->file);
+ dfm_close_reader (rf->reader);
+ free (rf->encoding);
+ free (rf);
+ }
+}
+
+static struct matrix_command *
+matrix_read_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_READ,
+ .read = { .format = FMT_F },
+ };
+
+ struct file_handle *fh = NULL;
+ char *encoding = NULL;
+ struct matrix_read *read = &cmd->read;
+ read->dst = matrix_lvalue_parse (s);
+ if (!read->dst)
+ goto error;
+
+ int by_ofs = 0;
+ int format_ofs = 0;
+ int record_width_start = 0, record_width_end = 0;
+
+ int by = 0;
+ int repetitions = 0;
+ int record_width = 0;
+ bool seen_format = false;
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "FILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (fh);
+ fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
+ if (!fh)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "ENCODING"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_string (s->lexer))
+ goto error;
+
+ free (encoding);
+ encoding = ss_xstrdup (lex_tokss (s->lexer));
+
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "FIELD"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ record_width_start = lex_ofs (s->lexer);
+ if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
+ goto error;
+ read->c1 = lex_integer (s->lexer);
+ lex_get (s->lexer);
+ if (!lex_force_match (s->lexer, T_TO)
+ || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
+ goto error;
+ read->c2 = lex_integer (s->lexer) + 1;
+ record_width_end = lex_ofs (s->lexer);
+ lex_get (s->lexer);
+
+ record_width = read->c2 - read->c1;
+ if (lex_match (s->lexer, T_BY))
+ {
+ if (!lex_force_int_range (s->lexer, "BY", 1,
+ read->c2 - read->c1))
+ goto error;
+ by = lex_integer (s->lexer);
+ by_ofs = lex_ofs (s->lexer);
+ int field_end = lex_ofs (s->lexer);
+ lex_get (s->lexer);
+
+ if (record_width % by)
+ {
+ lex_ofs_error (
+ s->lexer, record_width_start, field_end,
+ _("Field width %d does not evenly divide record width %d."),
+ by, record_width);
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates the record width."));
+ lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
+ _("This syntax specifies the field width."));
+ goto error;
+ }
+ }
+ else
+ by = 0;
+ }
+ else if (lex_match_id (s->lexer, "SIZE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ matrix_expr_destroy (read->size);
+ read->size = matrix_parse_exp (s);
+ if (!read->size)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "MODE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (lex_match_id (s->lexer, "RECTANGULAR"))
+ read->symmetric = false;
+ else if (lex_match_id (s->lexer, "SYMMETRIC"))
+ read->symmetric = true;
+ else
+ {
+ lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
+ goto error;
+ }
+ }
+ else if (lex_match_id (s->lexer, "REREAD"))
+ read->reread = true;
+ else if (lex_match_id (s->lexer, "FORMAT"))
+ {
+ if (seen_format)
+ {
+ lex_sbc_only_once (s->lexer, "FORMAT");
+ goto error;
+ }
+ seen_format = true;
+
+ lex_match (s->lexer, T_EQUALS);
+
+ if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
+ goto error;
+
+ format_ofs = lex_ofs (s->lexer);
+ const char *p = lex_tokcstr (s->lexer);
+ if (c_isdigit (p[0]))
+ {
+ repetitions = atoi (p);
+ p += strspn (p, "0123456789");
+ if (!fmt_from_name (p, &read->format))
+ {
+ lex_error (s->lexer, _("Unknown format %s."), p);
+ goto error;
+ }
+ lex_get (s->lexer);
+ }
+ else if (fmt_from_name (p, &read->format))
+ lex_get (s->lexer);
+ else
+ {
+ struct fmt_spec format;
+ if (!parse_format_specifier (s->lexer, &format))
+ goto error;
+ read->format = format.type;
+ read->w = format.w;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
+ "REREAD", "FORMAT");
+ goto error;
+ }
+ }
+
+ if (!read->c1)
+ {
+ lex_sbc_missing (s->lexer, "FIELD");
+ goto error;
+ }
+
+ if (!read->dst->n_indexes && !read->size)
+ {
+ msg (SE, _("SIZE is required for reading data into a full matrix "
+ "(as opposed to a submatrix)."));
+ msg_at (SN, read->dst->var_location,
+ _("This expression designates a full matrix."));
+ goto error;
+ }
+
+ if (!fh)
+ {
+ if (s->prev_read_file)
+ fh = fh_ref (s->prev_read_file);
+ else
+ {
+ lex_sbc_missing (s->lexer, "FILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_read_file);
+ s->prev_read_file = fh_ref (fh);
+
+ read->rf = read_file_create (s, fh);
+ fh = NULL;
+ if (encoding)
+ {
+ free (read->rf->encoding);
+ read->rf->encoding = encoding;
+ encoding = NULL;
+ }
+
+ /* Field width may be specified in multiple ways:
+
+ 1. BY on FIELD.
+ 2. The format on FORMAT.
+ 3. The repetition factor on FORMAT.
+
+ (2) and (3) are mutually exclusive.
+
+ If more than one of these is present, they must agree. If none of them is
+ present, then free-field format is used.
+ */
+ if (repetitions > record_width)
+ {
+ msg (SE, _("%d repetitions cannot fit in record width %d."),
+ repetitions, record_width);
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ _("This syntax designates the number of repetitions."));
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates the record width."));
+ goto error;
+ }
+ int w = (repetitions ? record_width / repetitions
+ : read->w ? read->w
+ : by);
+ if (by && w != by)
+ {
+ msg (SE, _("This command specifies two different field widths."));
+ if (repetitions)
+ {
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ ngettext ("This syntax specifies %d repetition.",
+ "This syntax specifies %d repetitions.",
+ repetitions),
+ repetitions);
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates record width %d, "
+ "which divided by %d repetitions implies "
+ "field width %d."),
+ record_width, repetitions, w);
+ }
+ else
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ _("This syntax specifies field width %d."), w);
+
+ lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
+ _("This syntax specifies field width %d."), by);
+ goto error;
+ }
+ read->w = w;
+ return cmd;
+
+error:
+ fh_unref (fh);
+ matrix_command_destroy (cmd);
+ free (encoding);
+ return NULL;
+}
+
+static void
+parse_error (const struct dfm_reader *reader, enum fmt_type format,
+ struct substring data, size_t y, size_t x,
+ int first_column, int last_column, char *error)
+{
+ int line_number = dfm_get_line_number (reader);
+ struct msg_location location = {
+ .file_name = intern_new (dfm_get_file_name (reader)),
+ .start = { .line = line_number, .column = first_column },
+ .end = { .line = line_number, .column = last_column },
+ };
+ msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
+ "for matrix row %zu, column %zu: %s"),
+ (int) data.length, data.string, fmt_name (format),
+ y + 1, x + 1, error);
+ msg_location_uninit (&location);
+ free (error);
+}
+
+static void
+matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
+ gsl_matrix *m, struct substring p, size_t y, size_t x,
+ const char *line_start)
+{
+ const char *input_encoding = dfm_reader_get_encoding (reader);
+ char *error;
+ double f;
+ if (fmt_is_numeric (read->format))
+ {
+ union value v;
+ error = data_in (p, input_encoding, read->format,
+ settings_get_fmt_settings (), &v, 0, NULL);
+ if (!error && v.f == SYSMIS)
+ error = xstrdup (_("Matrix data may not contain missing value."));
+ f = v.f;
+ }
+ else
+ {
+ uint8_t s[sizeof (double)];
+ union value v = { .s = s };
+ error = data_in (p, input_encoding, read->format,
+ settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
+ memcpy (&f, s, sizeof f);
+ }
+
+ if (error)
+ {
+ int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
+ int nc = ss_utf8_count_columns (p);
+ int c2 = c1 + MAX (1, nc) - 1;
+ parse_error (reader, read->format, p, y, x, c1, c2, error);
+ }
+ else
+ {
+ gsl_matrix_set (m, y, x, f);
+ if (read->symmetric && x != y)
+ gsl_matrix_set (m, x, y, f);
+ }
+}
+
+static bool
+matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
+ struct substring *line, const char **startp)
+{
+ struct matrix_read *read = &cmd->read;
+ if (dfm_eof (reader))
+ {
+ msg_at (SE, cmd->location,
+ _("Unexpected end of file reading matrix data."));
+ return false;
+ }
+ dfm_expand_tabs (reader);
+ struct substring record = dfm_get_record (reader);
+ /* XXX need to recode record into UTF-8 */
+ *startp = record.string;
+ *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
+ return true;
+}
+
+static void
+matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
+ gsl_matrix *m)
+{
+ struct matrix_read *read = &cmd->read;
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ size_t nx = read->symmetric ? y + 1 : m->size2;
+
+ struct substring line = ss_empty ();
+ const char *line_start = line.string;
+ for (size_t x = 0; x < nx; x++)
+ {
+ struct substring p;
+ if (!read->w)
+ {
+ for (;;)
+ {
+ ss_ltrim (&line, ss_cstr (" ,"));
+ if (!ss_is_empty (line))
+ break;
+ if (!matrix_read_line (cmd, reader, &line, &line_start))
+ return;
+ dfm_forward_record (reader);
+ }
+
+ ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
+ }
+ else
+ {
+ if (!matrix_read_line (cmd, reader, &line, &line_start))
+ return;
+ size_t fields_per_line = (read->c2 - read->c1) / read->w;
+ int f = x % fields_per_line;
+ if (f == fields_per_line - 1)
+ dfm_forward_record (reader);
+
+ p = ss_substr (line, read->w * f, read->w);
+ }
+
+ matrix_read_set_field (read, reader, m, p, y, x, line_start);
+ }
+
+ if (read->w)
+ dfm_forward_record (reader);
+ else
+ {
+ ss_ltrim (&line, ss_cstr (" ,"));
+ if (!ss_is_empty (line))
+ {
+ int line_number = dfm_get_line_number (reader);
+ int c1 = utf8_count_columns (line_start,
+ line.string - line_start) + 1;
+ int c2 = c1 + ss_utf8_count_columns (line) - 1;
+ struct msg_location location = {
+ .file_name = intern_new (dfm_get_file_name (reader)),
+ .start = { .line = line_number, .column = c1 },
+ .end = { .line = line_number, .column = c2 },
+ };
+ msg_at (DW, &location,
+ _("Trailing garbage following data for matrix row %zu."),
+ y + 1);
+ msg_location_uninit (&location);
+ }
+ }
+ }
+}
+
+static void
+matrix_read_execute (struct matrix_command *cmd)
+{
+ struct matrix_read *read = &cmd->read;
+ struct index_vector iv0, iv1;
+ if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
+ return;
+
+ size_t size[2] = { SIZE_MAX, SIZE_MAX };
+ if (read->size)
+ {
+ gsl_matrix *m = matrix_expr_evaluate (read->size);
+ if (!m)
+ return;
+
+ if (!is_vector (m))
+ {
+ msg_at (SE, matrix_expr_location (read->size),
+ _("SIZE must evaluate to a scalar or a 2-element vector, "
+ "not a %zu×%zu matrix."), m->size1, m->size2);
+ gsl_matrix_free (m);
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return;
+ }
+
+ gsl_vector v = to_vector (m);
+ double d[2];
+ if (v.size == 1)
+ {
+ d[0] = gsl_vector_get (&v, 0);
+ d[1] = 1;
+ }
+ else if (v.size == 2)
+ {
+ d[0] = gsl_vector_get (&v, 0);
+ d[1] = gsl_vector_get (&v, 1);
+ }
+ else
+ {
+ msg_at (SE, matrix_expr_location (read->size),
+ _("SIZE must evaluate to a scalar or a 2-element vector, "
+ "not a %zu×%zu matrix."),
+ m->size1, m->size2),
+ gsl_matrix_free (m);
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return;
+ }
+ gsl_matrix_free (m);
+
+ if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
+ {
+ msg_at (SE, matrix_expr_location (read->size),
+ _("Matrix dimensions %g×%g specified on SIZE "
+ "are outside valid range."),
+ d[0], d[1]);
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return;
+ }
+
+ size[0] = d[0];
+ size[1] = d[1];
+ }
+
+ if (read->dst->n_indexes)
+ {
+ size_t submatrix_size[2];
+ if (read->dst->n_indexes == 2)
+ {
+ submatrix_size[0] = iv0.n;
+ submatrix_size[1] = iv1.n;
+ }
+ else if (read->dst->var->value->size1 == 1)
+ {
+ submatrix_size[0] = 1;
+ submatrix_size[1] = iv0.n;
+ }
+ else
+ {
+ submatrix_size[0] = iv0.n;
+ submatrix_size[1] = 1;
+ }
+
+ if (read->size)
+ {
+ if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
+ {
+ msg_at (SE, cmd->location,
+ _("Dimensions specified on SIZE differ from dimensions "
+ "of destination submatrix."));
+ msg_at (SN, matrix_expr_location (read->size),
+ _("SIZE specifies dimensions %zu×%zu."),
+ size[0], size[1]);
+ msg_at (SN, read->dst->full_location,
+ _("Destination submatrix has dimensions %zu×%zu."),
+ submatrix_size[0], submatrix_size[1]);
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return;
+ }
+ }
+ else
+ {
+ size[0] = submatrix_size[0];
+ size[1] = submatrix_size[1];
+ }
+ }
+
+ struct dfm_reader *reader = read_file_open (read->rf);
+ if (read->reread)
+ dfm_reread_record (reader, 1);
+
+ if (read->symmetric && size[0] != size[1])
+ {
+ msg_at (SE, cmd->location,
+ _("Cannot read non-square %zu×%zu matrix "
+ "using READ with MODE=SYMMETRIC."),
+ size[0], size[1]);
+ index_vector_uninit (&iv0);
+ index_vector_uninit (&iv1);
+ return;
+ }
+ gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
+ matrix_read (cmd, reader, tmp);
+ matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
+}
+\f
+/* WRITE. */
+
+static struct write_file *
+write_file_create (struct matrix_state *s, struct file_handle *fh)
+{
+ for (size_t i = 0; i < s->n_write_files; i++)
+ {
+ struct write_file *wf = s->write_files[i];
+ if (wf->file == fh)
+ {
+ fh_unref (fh);
+ return wf;
+ }
+ }
+
+ struct write_file *wf = xmalloc (sizeof *wf);
+ *wf = (struct write_file) { .file = fh };
+
+ s->write_files = xrealloc (s->write_files,
+ (s->n_write_files + 1) * sizeof *s->write_files);
+ s->write_files[s->n_write_files++] = wf;
+ return wf;
+}
+
+static struct dfm_writer *
+write_file_open (struct write_file *wf)
+{
+ if (!wf->writer)
+ wf->writer = dfm_open_writer (wf->file, wf->encoding);
+ return wf->writer;
+}
+
+static void
+write_file_destroy (struct write_file *wf)
+{
+ if (wf)
+ {
+ if (wf->held)
+ {
+ dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
+ wf->held->s.ss.length);
+ u8_line_destroy (wf->held);
+ free (wf->held);
+ }
+
+ fh_unref (wf->file);
+ dfm_close_writer (wf->writer);
+ free (wf->encoding);
+ free (wf);
+ }
+}
+
+static struct matrix_command *
+matrix_write_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_WRITE,
+ };
+
+ struct file_handle *fh = NULL;
+ char *encoding = NULL;
+ struct matrix_write *write = &cmd->write;
+ write->expression = matrix_parse_exp (s);
+ if (!write->expression)
+ goto error;
+
+ int by_ofs = 0;
+ int format_ofs = 0;
+ int record_width_start = 0, record_width_end = 0;
+
+ int by = 0;
+ int repetitions = 0;
+ int record_width = 0;
+ enum fmt_type format = FMT_F;
+ bool has_format = false;
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "OUTFILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (fh);
+ fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
+ if (!fh)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "ENCODING"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_string (s->lexer))
+ goto error;
+
+ free (encoding);
+ encoding = ss_xstrdup (lex_tokss (s->lexer));
+
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "FIELD"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ record_width_start = lex_ofs (s->lexer);
+
+ if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
+ goto error;
+ write->c1 = lex_integer (s->lexer);
+ lex_get (s->lexer);
+ if (!lex_force_match (s->lexer, T_TO)
+ || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
+ goto error;
+ write->c2 = lex_integer (s->lexer) + 1;
+ record_width_end = lex_ofs (s->lexer);
+ lex_get (s->lexer);
+
+ record_width = write->c2 - write->c1;
+ if (lex_match (s->lexer, T_BY))
+ {
+ if (!lex_force_int_range (s->lexer, "BY", 1,
+ write->c2 - write->c1))
+ goto error;
+ by_ofs = lex_ofs (s->lexer);
+ int field_end = lex_ofs (s->lexer);
+ by = lex_integer (s->lexer);
+ lex_get (s->lexer);
+
+ if (record_width % by)
+ {
+ lex_ofs_error (
+ s->lexer, record_width_start, field_end,
+ _("Field width %d does not evenly divide record width %d."),
+ by, record_width);
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates the record width."));
+ lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
+ _("This syntax specifies the field width."));
+ goto error;
+ }
+ }
+ else
+ by = 0;
+ }
+ else if (lex_match_id (s->lexer, "MODE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (lex_match_id (s->lexer, "RECTANGULAR"))
+ write->triangular = false;
+ else if (lex_match_id (s->lexer, "TRIANGULAR"))
+ write->triangular = true;
+ else
+ {
+ lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
+ goto error;
+ }
+ }
+ else if (lex_match_id (s->lexer, "HOLD"))
+ write->hold = true;
+ else if (lex_match_id (s->lexer, "FORMAT"))
+ {
+ if (has_format || write->format)
+ {
+ lex_sbc_only_once (s->lexer, "FORMAT");
+ goto error;
+ }
+
+ lex_match (s->lexer, T_EQUALS);
+
+ if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
+ goto error;
+
+ format_ofs = lex_ofs (s->lexer);
+ const char *p = lex_tokcstr (s->lexer);
+ if (c_isdigit (p[0]))
+ {
+ repetitions = atoi (p);
+ p += strspn (p, "0123456789");
+ if (!fmt_from_name (p, &format))
+ {
+ lex_error (s->lexer, _("Unknown format %s."), p);
+ goto error;
+ }
+ has_format = true;
+ lex_get (s->lexer);
+ }
+ else if (fmt_from_name (p, &format))
+ {
+ has_format = true;
+ lex_get (s->lexer);
+ }
+ else
+ {
+ struct fmt_spec spec;
+ if (!parse_format_specifier (s->lexer, &spec))
+ goto error;
+ write->format = xmemdup (&spec, sizeof spec);
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
+ "HOLD", "FORMAT");
+ goto error;
+ }
+ }
+
+ if (!write->c1)
+ {
+ lex_sbc_missing (s->lexer, "FIELD");
+ goto error;
+ }
+
+ if (!fh)
+ {
+ if (s->prev_write_file)
+ fh = fh_ref (s->prev_write_file);
+ else
+ {
+ lex_sbc_missing (s->lexer, "OUTFILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_write_file);
+ s->prev_write_file = fh_ref (fh);
+
+ write->wf = write_file_create (s, fh);
+ fh = NULL;
+ if (encoding)
+ {
+ free (write->wf->encoding);
+ write->wf->encoding = encoding;
+ encoding = NULL;
+ }
+
+ /* Field width may be specified in multiple ways:
+
+ 1. BY on FIELD.
+ 2. The format on FORMAT.
+ 3. The repetition factor on FORMAT.
+
+ (2) and (3) are mutually exclusive.
+
+ If more than one of these is present, they must agree. If none of them is
+ present, then free-field format is used.
+ */
+ if (repetitions > record_width)
+ {
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ _("This syntax designates the number of repetitions."));
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates the record width."));
+ goto error;
+ }
+ int w = (repetitions ? record_width / repetitions
+ : write->format ? write->format->w
+ : by);
+ if (by && w != by)
+ {
+ msg (SE, _("This command specifies two different field widths."));
+ if (repetitions)
+ {
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ ngettext ("This syntax specifies %d repetition.",
+ "This syntax specifies %d repetitions.",
+ repetitions),
+ repetitions);
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates record width %d, "
+ "which divided by %d repetitions implies "
+ "field width %d."),
+ record_width, repetitions, w);
+ }
+ else
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ _("This syntax specifies field width %d."), w);
+
+ lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
+ _("This syntax specifies field width %d."), by);
+ goto error;
+ }
+ if (w && !write->format)
+ {
+ write->format = xmalloc (sizeof *write->format);
+ *write->format = (struct fmt_spec) { .type = format, .w = w };
+
+ char *error = fmt_check_output__ (write->format);
+ if (error)
+ {
+ msg (SE, "%s", error);
+ free (error);
+
+ if (has_format)
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ _("This syntax specifies format %s."),
+ fmt_name (format));
+
+ if (repetitions)
+ {
+ lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
+ ngettext ("This syntax specifies %d repetition.",
+ "This syntax specifies %d repetitions.",
+ repetitions),
+ repetitions);
+ lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
+ _("This syntax designates record width %d, "
+ "which divided by %d repetitions implies "
+ "field width %d."),
+ record_width, repetitions, w);
+ }
+
+ if (by)
+ lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
+ _("This syntax specifies field width %d."), by);
+
+ goto error;
+ }
+ }
+
+ if (write->format && fmt_var_width (write->format) > sizeof (double))
+ {
+ char fs[FMT_STRING_LEN_MAX + 1];
+ fmt_to_string (write->format, fs);
+ lex_ofs_error (s->lexer, format_ofs, format_ofs,
+ _("Format %s is too wide for %zu-byte matrix elements."),
+ fs, sizeof (double));
+ goto error;
+ }
+
+ return cmd;
+
+error:
+ fh_unref (fh);
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_write_execute (struct matrix_write *write)
+{
+ gsl_matrix *m = matrix_expr_evaluate (write->expression);
+ if (!m)
+ return;
+
+ if (write->triangular && m->size1 != m->size2)
+ {
+ msg_at (SE, matrix_expr_location (write->expression),
+ _("WRITE with MODE=TRIANGULAR requires a square matrix but "
+ "the matrix to be written has dimensions %zu×%zu."),
+ m->size1, m->size2);
+ gsl_matrix_free (m);
+ return;
+ }
+
+ struct dfm_writer *writer = write_file_open (write->wf);
+ if (!writer || !m->size1)
+ {
+ gsl_matrix_free (m);
+ return;
+ }
+
+ const struct fmt_settings *settings = settings_get_fmt_settings ();
+ struct u8_line *line = write->wf->held;
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ if (!line)
+ {
+ line = xmalloc (sizeof *line);
+ u8_line_init (line);
+ }
+ size_t nx = write->triangular ? y + 1 : m->size2;
+ int x0 = write->c1;
+ for (size_t x = 0; x < nx; x++)
+ {
+ char *s;
+ double f = gsl_matrix_get (m, y, x);
+ if (write->format)
+ {
+ union value v;
+ if (fmt_is_numeric (write->format->type))
+ v.f = f;
+ else
+ v.s = (uint8_t *) &f;
+ s = data_out (&v, NULL, write->format, settings);
+ }
+ else
+ {
+ s = xmalloc (DBL_BUFSIZE_BOUND);
+ if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
+ >= DBL_BUFSIZE_BOUND)
+ abort ();
+ }
+ size_t len = strlen (s);
+ int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
+ if (width + x0 > write->c2)
+ {
+ dfm_put_record_utf8 (writer, line->s.ss.string,
+ line->s.ss.length);
+ u8_line_clear (line);
+ x0 = write->c1;
+ }
+ u8_line_put (line, x0, x0 + width, s, len);
+ free (s);
+
+ x0 += write->format ? write->format->w : width + 1;
+ }
+
+ if (y + 1 >= m->size1 && write->hold)
+ break;
+ dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
+ u8_line_clear (line);
+ }
+ if (!write->hold)
+ {
+ u8_line_destroy (line);
+ free (line);
+ line = NULL;
+ }
+ write->wf->held = line;
+
+ gsl_matrix_free (m);
+}
+\f
+/* GET. */
+
+static struct matrix_command *
+matrix_get_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_GET,
+ .get = {
+ .lexer = s->lexer,
+ .dataset = s->dataset,
+ .user = { .treatment = MGET_ERROR },
+ .system = { .treatment = MGET_ERROR },
+ }
+ };
+
+ struct matrix_get *get = &cmd->get;
+ get->dst = matrix_lvalue_parse (s);
+ if (!get->dst)
+ goto error;
+
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "FILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (get->file);
+ if (lex_match (s->lexer, T_ASTERISK))
+ get->file = NULL;
+ else
+ {
+ get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
+ if (!get->file)
+ goto error;
+ }
+ }
+ else if (lex_match_id (s->lexer, "ENCODING"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_string (s->lexer))
+ goto error;
+
+ free (get->encoding);
+ get->encoding = ss_xstrdup (lex_tokss (s->lexer));
+
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "VARIABLES"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ if (get->n_vars)
+ {
+ lex_sbc_only_once (s->lexer, "VARIABLES");
+ goto error;
+ }
+
+ if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "NAMES"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_id (s->lexer))
+ goto error;
+
+ struct substring name = lex_tokss (s->lexer);
+ get->names = matrix_var_lookup (s, name);
+ if (!get->names)
+ get->names = matrix_var_create (s, name);
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "MISSING"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (lex_match_id (s->lexer, "ACCEPT"))
+ get->user.treatment = MGET_ACCEPT;
+ else if (lex_match_id (s->lexer, "OMIT"))
+ get->user.treatment = MGET_OMIT;
+ else if (lex_is_number (s->lexer))
+ {
+ get->user.treatment = MGET_RECODE;
+ get->user.substitute = lex_number (s->lexer);
+ lex_get (s->lexer);
+ }
+ else
+ {
+ lex_error (s->lexer, _("Syntax error expecting ACCEPT or OMIT or "
+ "a number for MISSING."));
+ goto error;
+ }
+ }
+ else if (lex_match_id (s->lexer, "SYSMIS"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (lex_match_id (s->lexer, "OMIT"))
+ get->system.treatment = MGET_OMIT;
+ else if (lex_is_number (s->lexer))
+ {
+ get->system.treatment = MGET_RECODE;
+ get->system.substitute = lex_number (s->lexer);
+ lex_get (s->lexer);
+ }
+ else
+ {
+ lex_error (s->lexer, _("Syntax error expecting OMIT or a number "
+ "for SYSMIS."));
+ goto error;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
+ "MISSING", "SYSMIS");
+ goto error;
+ }
+ }
+
+ if (get->user.treatment != MGET_ACCEPT)
+ get->system.treatment = MGET_ERROR;
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
+ const struct dictionary *dict)
+{
+ struct matrix_get *get = &cmd->get;
+ struct variable **vars;
+ size_t n_vars = 0;
+
+ if (get->n_vars)
+ {
+ if (!var_syntax_evaluate (get->lexer, get->vars, get->n_vars, dict,
+ &vars, &n_vars, PV_NUMERIC))
+ return;
+ }
+ else
+ {
+ n_vars = dict_get_n_vars (dict);
+ vars = xnmalloc (n_vars, sizeof *vars);
+ for (size_t i = 0; i < n_vars; i++)
+ {
+ struct variable *var = dict_get_var (dict, i);
+ if (!var_is_numeric (var))
+ {
+ msg_at (SE, cmd->location, _("Variable %s is not numeric."),
+ var_get_name (var));
+ free (vars);
+ return;
+ }
+ vars[i] = var;
+ }
+ }
+
+ if (get->names)
+ {
+ gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
+ for (size_t i = 0; i < n_vars; i++)
+ {
+ char s[sizeof (double)];
+ double f;
+ buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
+ memcpy (&f, s, sizeof f);
+ gsl_matrix_set (names, i, 0, f);
+ }
+
+ gsl_matrix_free (get->names->value);
+ get->names->value = names;
+ }
+
+ size_t n_rows = 0;
+ gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
+ long long int casenum = 1;
+ bool error = false;
+ for (struct ccase *c = casereader_read (reader); c;
+ c = casereader_read (reader), casenum++)
+ {
+ if (n_rows >= m->size1)
+ {
+ gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
+ for (size_t y = 0; y < n_rows; y++)
+ for (size_t x = 0; x < n_vars; x++)
+ gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
+ gsl_matrix_free (m);
+ m = p;
+ }
+
+ bool keep = true;
+ for (size_t x = 0; x < n_vars; x++)
+ {
+ const struct variable *var = vars[x];
+ double d = case_num (c, var);
+ if (d == SYSMIS)
+ {
+ if (get->system.treatment == MGET_RECODE)
+ d = get->system.substitute;
+ else if (get->system.treatment == MGET_OMIT)
+ keep = false;
+ else
+ {
+ msg_at (SE, cmd->location, _("Variable %s in case %lld "
+ "is system-missing."),
+ var_get_name (var), casenum);
+ error = true;
+ }
+ }
+ else if (var_is_num_missing (var, d) == MV_USER)
+ {
+ if (get->user.treatment == MGET_RECODE)
+ d = get->user.substitute;
+ else if (get->user.treatment == MGET_OMIT)
+ keep = false;
+ else if (get->user.treatment != MGET_ACCEPT)
+ {
+ msg_at (SE, cmd->location,
+ _("Variable %s in case %lld has user-missing "
+ "value %g."),
+ var_get_name (var), casenum, d);
+ error = true;
+ }
+ }
+ gsl_matrix_set (m, n_rows, x, d);
+ }
+ case_unref (c);
+ if (error)
+ break;
+ if (keep)
+ n_rows++;
+ }
+ if (!error)
+ {
+ m->size1 = n_rows;
+ matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
+ }
+ else
+ gsl_matrix_free (m);
+ free (vars);
+}
+
+static bool
+matrix_open_casereader (const struct matrix_command *cmd,
+ const char *command_name, struct file_handle *file,
+ const char *encoding, struct dataset *dataset,
+ struct casereader **readerp, struct dictionary **dictp)
+{
+ if (file)
+ {
+ *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
+ return *readerp != NULL;
+ }
+ else
+ {
+ if (dict_get_n_vars (dataset_dict (dataset)) == 0)
+ {
+ msg_at (SE, cmd->location,
+ _("The %s command cannot read an empty active file."),
+ command_name);
+ return false;
+ }
+ *readerp = proc_open (dataset);
+ *dictp = dict_ref (dataset_dict (dataset));
+ return true;
+ }
+}
+
+static void
+matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
+ struct casereader *reader, struct dictionary *dict)
+{
+ dict_unref (dict);
+ casereader_destroy (reader);
+ if (!file)
+ proc_commit (dataset);
+}
+
+static void
+matrix_get_execute (struct matrix_command *cmd)
+{
+ struct matrix_get *get = &cmd->get;
+ struct casereader *r;
+ struct dictionary *d;
+ if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
+ get->dataset, &r, &d))
+ {
+ matrix_get_execute__ (cmd, r, d);
+ matrix_close_casereader (get->file, get->dataset, r, d);
+ }
+}
+\f
+/* MSAVE. */
+
+static bool
+variables_changed (const char *keyword,
+ const struct string_array *new_vars,
+ const struct msg_location *new_vars_location,
+ const struct msg_location *new_location,
+ const struct string_array *old_vars,
+ const struct msg_location *old_vars_location,
+ const struct msg_location *old_location)
+{
+ if (new_vars->n)
+ {
+ if (!old_vars->n)
+ {
+ msg_at (SE, new_location,
+ _("%s may only be specified on MSAVE if it was specified "
+ "on the first MSAVE within MATRIX."), keyword);
+ msg_at (SN, old_location,
+ _("The first MSAVE in MATRIX did not specify %s."),
+ keyword);
+ msg_at (SN, new_vars_location,
+ _("This is the specification of %s on a later MSAVE."),
+ keyword);
+ return true;
+ }
+ if (!string_array_equal_case (old_vars, new_vars))
+ {
+ msg_at (SE, new_location,
+ _("%s must specify the same variables on each MSAVE "
+ "within a given MATRIX."), keyword);
+ msg_at (SE, old_vars_location,
+ _("This is the specification of %s on the first MSAVE."),
+ keyword);
+ msg_at (SE, new_vars_location,
+ _("This is a different specification of %s on a later MSAVE."),
+ keyword);
+ return true;
+ }
+ }
+ return false;
+}
+
+static bool
+msave_common_changed (const struct msave_common *old,
+ const struct msave_common *new)
+{
+ if (new->outfile && !fh_equal (old->outfile, new->outfile))
+ {
+ msg (SE, _("OUTFILE must name the same file on each MSAVE "
+ "within a single MATRIX command."));
+ msg_at (SN, old->outfile_location,
+ _("This is the OUTFILE on the first MSAVE command."));
+ msg_at (SN, new->outfile_location,
+ _("This is the OUTFILE on a later MSAVE command."));
+ return false;
+ }
+
+ if (!variables_changed ("VARIABLES",
+ &new->variables, new->variables_location, new->location,
+ &old->variables, old->variables_location, old->location)
+ && !variables_changed ("FNAMES",
+ &new->fnames, new->fnames_location, new->location,
+ &old->fnames, old->fnames_location, old->location)
+ && !variables_changed ("SNAMES",
+ &new->snames, new->snames_location, new->location,
+ &old->snames, old->snames_location, old->location))
+ return false;
+
+ return true;
+}
+
+static void
+msave_common_destroy (struct msave_common *common)
+{
+ if (common)
+ {
+ msg_location_destroy (common->location);
+ fh_unref (common->outfile);
+ msg_location_destroy (common->outfile_location);
+ string_array_destroy (&common->variables);
+ msg_location_destroy (common->variables_location);
+ string_array_destroy (&common->fnames);
+ msg_location_destroy (common->fnames_location);
+ string_array_destroy (&common->snames);
+ msg_location_destroy (common->snames_location);
+
+ for (size_t i = 0; i < common->n_factors; i++)
+ matrix_expr_destroy (common->factors[i]);
+ free (common->factors);
+
+ for (size_t i = 0; i < common->n_splits; i++)
+ matrix_expr_destroy (common->splits[i]);
+ free (common->splits);
+
+ dict_unref (common->dict);
+ casewriter_destroy (common->writer);
+
+ free (common);
+ }
+}
+
+static const char *
+match_rowtype (struct lexer *lexer)
+{
+ static const char *rowtypes[] = {
+ "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
+ };
+ size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
+
+ for (size_t i = 0; i < n_rowtypes; i++)
+ if (lex_match_id (lexer, rowtypes[i]))
+ return rowtypes[i];
+
+ lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
+ return NULL;
+}
+
+static bool
+parse_var_names (struct lexer *lexer, struct string_array *sa,
+ struct msg_location **locationp)
+{
+ lex_match (lexer, T_EQUALS);
+
+ string_array_clear (sa);
+ msg_location_destroy (*locationp);
+ *locationp = NULL;
+
+ struct dictionary *dict = dict_create (get_default_encoding ());
+ char **names;
+ size_t n_names;
+ int start_ofs = lex_ofs (lexer);
+ bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
+ PV_NO_DUPLICATE | PV_NO_SCRATCH);
+ int end_ofs = lex_ofs (lexer) - 1;
+ dict_unref (dict);
+
+ if (ok)
+ {
+ for (size_t i = 0; i < n_names; i++)
+ if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
+ || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
+ {
+ lex_ofs_error (lexer, start_ofs, end_ofs,
+ _("Variable name %s is reserved."), names[i]);
+ for (size_t j = 0; j < n_names; j++)
+ free (names[i]);
+ free (names);
+ return false;
+ }
+
+ sa->strings = names;
+ sa->n = sa->allocated = n_names;
+ *locationp = lex_ofs_location (lexer, start_ofs, end_ofs);
+ }
+ return ok;
+}
+
+static struct matrix_command *
+matrix_msave_parse (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+
+ struct msave_common *common = xmalloc (sizeof *common);
+ *common = (struct msave_common) { .outfile = NULL };
+
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
+
+ struct matrix_expr *splits = NULL;
+ struct matrix_expr *factors = NULL;
+
+ struct matrix_msave *msave = &cmd->msave;
+ msave->expr = matrix_parse_exp (s);
+ if (!msave->expr)
+ goto error;
+
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "TYPE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ msave->rowtype = match_rowtype (s->lexer);
+ if (!msave->rowtype)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "OUTFILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (common->outfile);
+ int start_ofs = lex_ofs (s->lexer);
+ common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
+ if (!common->outfile)
+ goto error;
+ msg_location_destroy (common->outfile_location);
+ common->outfile_location = lex_ofs_location (s->lexer, start_ofs,
+ lex_ofs (s->lexer) - 1);
+ }
+ else if (lex_match_id (s->lexer, "VARIABLES"))
+ {
+ if (!parse_var_names (s->lexer, &common->variables,
+ &common->variables_location))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "FNAMES"))
+ {
+ if (!parse_var_names (s->lexer, &common->fnames,
+ &common->fnames_location))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "SNAMES"))
+ {
+ if (!parse_var_names (s->lexer, &common->snames,
+ &common->snames_location))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "SPLIT"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ matrix_expr_destroy (splits);
+ splits = matrix_parse_exp (s);
+ if (!splits)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "FACTOR"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ matrix_expr_destroy (factors);
+ factors = matrix_parse_exp (s);
+ if (!factors)
+ goto error;
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
+ "FNAMES", "SNAMES", "SPLIT", "FACTOR");
+ goto error;
+ }
+ }
+ if (!msave->rowtype)
+ {
+ lex_sbc_missing (s->lexer, "TYPE");
+ goto error;
+ }
+
+ if (!s->msave_common)
+ {
+ if (common->fnames.n && !factors)
+ {
+ msg_at (SE, common->fnames_location, _("FNAMES requires FACTOR."));
+ goto error;
+ }
+ if (common->snames.n && !splits)
+ {
+ msg_at (SE, common->snames_location, _("SNAMES requires SPLIT."));
+ goto error;
+ }
+ if (!common->outfile)
+ {
+ lex_sbc_missing (s->lexer, "OUTFILE");
+ goto error;
+ }
+ common->location = lex_ofs_location (s->lexer, start_ofs,
+ lex_ofs (s->lexer));
+ msg_location_remove_columns (common->location);
+ s->msave_common = common;
+ }
+ else
+ {
+ if (msave_common_changed (s->msave_common, common))
+ goto error;
+ msave_common_destroy (common);
+ }
+ msave->common = s->msave_common;
+
+ struct msave_common *c = s->msave_common;
+ if (factors)
+ {
+ if (c->n_factors >= c->allocated_factors)
+ c->factors = x2nrealloc (c->factors, &c->allocated_factors,
+ sizeof *c->factors);
+ c->factors[c->n_factors++] = factors;
+ }
+ if (c->n_factors > 0)
+ msave->factors = c->factors[c->n_factors - 1];
+
+ if (splits)
+ {
+ if (c->n_splits >= c->allocated_splits)
+ c->splits = x2nrealloc (c->splits, &c->allocated_splits,
+ sizeof *c->splits);
+ c->splits[c->n_splits++] = splits;
+ }
+ if (c->n_splits > 0)
+ msave->splits = c->splits[c->n_splits - 1];
+
+ return cmd;
+
+error:
+ matrix_expr_destroy (splits);
+ matrix_expr_destroy (factors);
+ msave_common_destroy (common);
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static gsl_vector *
+matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
+{
+ gsl_matrix *m = matrix_expr_evaluate (e);
+ if (!m)
+ return NULL;
+
+ if (!is_vector (m))
+ {
+ msg_at (SE, matrix_expr_location (e),
+ _("%s expression must evaluate to vector, "
+ "not a %zu×%zu matrix."),
+ name, m->size1, m->size2);
+ gsl_matrix_free (m);
+ return NULL;
+ }
+
+ return matrix_to_vector (m);
+}
+
+static const char *
+msave_add_vars (struct dictionary *d, const struct string_array *vars)
+{
+ for (size_t i = 0; i < vars->n; i++)
+ if (!dict_create_var (d, vars->strings[i], 0))
+ return vars->strings[i];
+ return NULL;
+}
+
+static struct dictionary *
+msave_create_dict (const struct msave_common *common)
+{
+ struct dictionary *dict = dict_create (get_default_encoding ());
+
+ const char *dup_split = msave_add_vars (dict, &common->snames);
+ if (dup_split)
+ {
+ /* Should not be possible because the parser ensures that the names are
+ unique. */
+ NOT_REACHED ();
+ }
+
+ dict_create_var_assert (dict, "ROWTYPE_", 8);
+
+ const char *dup_factor = msave_add_vars (dict, &common->fnames);
+ if (dup_factor)
+ {
+ msg_at (SE, common->fnames_location,
+ _("Duplicate or invalid FACTOR variable name %s."),
+ dup_factor);
+ goto error;
+ }
+
+ dict_create_var_assert (dict, "VARNAME_", 8);
+
+ const char *dup_var = msave_add_vars (dict, &common->variables);
+ if (dup_var)
+ {
+ msg_at (SE, common->variables_location,
+ _("Duplicate or invalid variable name %s."),
+ dup_var);
+ goto error;
+ }
+
+ return dict;
+
+error:
+ dict_unref (dict);
+ return NULL;
+}
+
+static void
+matrix_msave_execute (struct matrix_command *cmd)
+{
+ struct matrix_msave *msave = &cmd->msave;
+ struct msave_common *common = msave->common;
+ gsl_matrix *m = NULL;
+ gsl_vector *factors = NULL;
+ gsl_vector *splits = NULL;
+
+ m = matrix_expr_evaluate (msave->expr);
+ if (!m)
+ goto error;
+
+ if (!common->variables.n)
+ for (size_t i = 0; i < m->size2; i++)
+ string_array_append_nocopy (&common->variables,
+ xasprintf ("COL%zu", i + 1));
+ else if (m->size2 != common->variables.n)
+ {
+ msg_at (SE, matrix_expr_location (msave->expr),
+ _("Matrix on MSAVE has %zu columns but there are %zu variables."),
+ m->size2, common->variables.n);
+ goto error;
+ }
+
+ if (msave->factors)
+ {
+ factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
+ if (!factors)
+ goto error;
+
+ if (!common->fnames.n)
+ for (size_t i = 0; i < factors->size; i++)
+ string_array_append_nocopy (&common->fnames,
+ xasprintf ("FAC%zu", i + 1));
+ else if (factors->size != common->fnames.n)
+ {
+ msg_at (SE, matrix_expr_location (msave->factors),
+ _("There are %zu factor variables, "
+ "but %zu factor values were supplied."),
+ common->fnames.n, factors->size);
+ goto error;
+ }
+ }
+
+ if (msave->splits)
+ {
+ splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
+ if (!splits)
+ goto error;
+
+ if (!common->snames.n)
+ for (size_t i = 0; i < splits->size; i++)
+ string_array_append_nocopy (&common->snames,
+ xasprintf ("SPL%zu", i + 1));
+ else if (splits->size != common->snames.n)
+ {
+ msg_at (SE, matrix_expr_location (msave->splits),
+ _("There are %zu split variables, "
+ "but %zu split values were supplied."),
+ common->snames.n, splits->size);
+ goto error;
+ }
+ }
+
+ if (!common->writer)
+ {
+ struct dictionary *dict = msave_create_dict (common);
+ if (!dict)
+ goto error;
+
+ common->writer = any_writer_open (common->outfile, dict);
+ if (!common->writer)
+ {
+ dict_unref (dict);
+ goto error;
+ }
+
+ common->dict = dict;
+ }
+
+ bool matrix = (!strcmp (msave->rowtype, "COV")
+ || !strcmp (msave->rowtype, "CORR"));
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ struct ccase *c = case_create (dict_get_proto (common->dict));
+ size_t idx = 0;
+
+ /* Split variables */
+ if (splits)
+ for (size_t i = 0; i < splits->size; i++)
+ case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
+
+ /* ROWTYPE_. */
+ buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
+ msave->rowtype, ' ');
+
+ /* Factors. */
+ if (factors)
+ for (size_t i = 0; i < factors->size; i++)
+ *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
+
+ /* VARNAME_. */
+ const char *varname_ = (matrix && y < common->variables.n
+ ? common->variables.strings[y]
+ : "");
+ buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
+ varname_, ' ');
+
+ /* Continuous variables. */
+ for (size_t x = 0; x < m->size2; x++)
+ case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
+ casewriter_write (common->writer, c);
+ }
+
+error:
+ gsl_matrix_free (m);
+ gsl_vector_free (factors);
+ gsl_vector_free (splits);
+}
+\f
+/* MGET. */
+
+static struct matrix_command *
+matrix_mget_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_MGET,
+ .mget = {
+ .state = s,
+ .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
+ },
+ };
+
+ struct matrix_mget *mget = &cmd->mget;
+
+ lex_match (s->lexer, T_SLASH);
+ while (lex_token (s->lexer) != T_ENDCMD)
+ {
+ if (lex_match_id (s->lexer, "FILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (mget->file);
+ mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
+ if (!mget->file)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "ENCODING"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ if (!lex_force_string (s->lexer))
+ goto error;
+
+ free (mget->encoding);
+ mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
+
+ lex_get (s->lexer);
+ }
+ else if (lex_match_id (s->lexer, "TYPE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ while (lex_token (s->lexer) != T_SLASH
+ && lex_token (s->lexer) != T_ENDCMD)
+ {
+ const char *rowtype = match_rowtype (s->lexer);
+ if (!rowtype)
+ goto error;
+
+ stringi_set_insert (&mget->rowtypes, rowtype);
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FILE", "TYPE");
+ goto error;
+ }
+ lex_match (s->lexer, T_SLASH);
+ }
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static const struct variable *
+get_a8_var (const struct msg_location *loc,
+ const struct dictionary *d, const char *name)
+{
+ const struct variable *v = dict_lookup_var (d, name);
+ if (!v)
+ {
+ msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
+ return NULL;
+ }
+ if (var_get_width (v) != 8)
+ {
+ msg_at (SE, loc, _("%s variable in matrix data file must be "
+ "8-byte string, but it has width %d."),
+ name, var_get_width (v));
+ return NULL;
+ }
+ return v;
+}
+
+static bool
+var_changed (const struct ccase *ca, const struct ccase *cb,
+ const struct variable *var)
+{
+ return (ca && cb
+ ? !value_equal (case_data (ca, var), case_data (cb, var),
+ var_get_width (var))
+ : ca || cb);
+}
+
+static bool
+vars_changed (const struct ccase *ca, const struct ccase *cb,
+ const struct dictionary *d,
+ size_t first_var, size_t n_vars)
+{
+ for (size_t i = 0; i < n_vars; i++)
+ {
+ const struct variable *v = dict_get_var (d, first_var + i);
+ if (var_changed (ca, cb, v))
+ return true;
+ }
+ return false;
+}
+
+static bool
+vars_all_missing (const struct ccase *c, const struct dictionary *d,
+ size_t first_var, size_t n_vars)
+{
+ for (size_t i = 0; i < n_vars; i++)
+ if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
+ return false;
+ return true;
+}
+
+static void
+matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
+ const struct dictionary *d,
+ const struct variable *rowtype_var,
+ const struct stringi_set *accepted_rowtypes,
+ struct matrix_state *s,
+ size_t ss, size_t sn, size_t si,
+ size_t fs, size_t fn, size_t fi,
+ size_t cs, size_t cn,
+ struct pivot_table *pt,
+ struct pivot_dimension *var_dimension)
+{
+ if (!n_rows)
+ goto exit;
+
+ /* Is this a matrix for pooled data, either where there are no factor
+ variables or the factor variables are missing? */
+ bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
+
+ struct substring rowtype = case_ss (rows[0], rowtype_var);
+ ss_rtrim (&rowtype, ss_cstr (" "));
+ if (!stringi_set_is_empty (accepted_rowtypes)
+ && !stringi_set_contains_len (accepted_rowtypes,
+ rowtype.string, rowtype.length))
+ goto exit;
+
+ const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
+ : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
+ : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
+ : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
+ : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
+ : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
+ : NULL);
+ if (!prefix)
+ {
+ msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
+ (int) rowtype.length, rowtype.string);
+ goto exit;
+ }
+
+ struct string name = DS_EMPTY_INITIALIZER;
+ ds_put_cstr (&name, prefix);
+ if (!pooled)
+ ds_put_format (&name, "F%zu", fi);
+ if (si > 0)
+ ds_put_format (&name, "S%zu", si);
+
+ struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
+ if (!mv)
+ mv = matrix_var_create (s, ds_ss (&name));
+ else if (mv->value)
+ {
+ msg (SW, _("Matrix data file contains variable with existing name %s."),
+ ds_cstr (&name));
+ goto exit_free_name;
+ }
+
+ gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
+ size_t n_missing = 0;
+ for (size_t y = 0; y < n_rows; y++)
+ {
+ for (size_t x = 0; x < cn; x++)
+ {
+ struct variable *var = dict_get_var (d, cs + x);
+ double value = case_num (rows[y], var);
+ if (var_is_num_missing (var, value))
+ {
+ n_missing++;
+ value = 0.0;
+ }
+ gsl_matrix_set (m, y, x, value);
+ }
+ }
+
+ int var_index = pivot_category_create_leaf (
+ var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
+ double values[] = { n_rows, cn };
+ for (size_t j = 0; j < sn; j++)
+ {
+ struct variable *var = dict_get_var (d, ss + j);
+ const union value *value = case_data (rows[0], var);
+ pivot_table_put2 (pt, j, var_index,
+ pivot_value_new_var_value (var, value));
+ }
+ for (size_t j = 0; j < fn; j++)
+ {
+ struct variable *var = dict_get_var (d, fs + j);
+ const union value sysmis = { .f = SYSMIS };
+ const union value *value = pooled ? &sysmis : case_data (rows[0], var);
+ pivot_table_put2 (pt, j + sn, var_index,
+ pivot_value_new_var_value (var, value));
+ }
+ for (size_t j = 0; j < sizeof values / sizeof *values; j++)
+ pivot_table_put2 (pt, j + sn + fn, var_index,
+ pivot_value_new_integer (values[j]));
+
+ if (n_missing)
+ msg (SE, ngettext ("Matrix data file variable %s contains a missing "
+ "value, which was treated as zero.",
+ "Matrix data file variable %s contains %zu missing "
+ "values, which were treated as zero.", n_missing),
+ ds_cstr (&name), n_missing);
+ mv->value = m;
+
+exit_free_name:
+ ds_destroy (&name);
+
+exit:
+ for (size_t y = 0; y < n_rows; y++)
+ case_unref (rows[y]);
+}
+
+static void
+matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
+ const struct dictionary *d)
+{
+ struct matrix_mget *mget = &cmd->mget;
+ const struct msg_location *loc = cmd->location;
+ const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
+ const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
+ if (!rowtype_ || !varname_)
+ return;
+
+ if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
+ {
+ msg_at (SE, loc,
+ _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
+ return;
+ }
+ if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
+ {
+ msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
+ return;
+ }
+
+ for (size_t i = 0; i < dict_get_n_vars (d); i++)
+ {
+ const struct variable *v = dict_get_var (d, i);
+ if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
+ {
+ msg_at (SE, loc,
+ _("Matrix data file contains unexpected string variable %s."),
+ var_get_name (v));
+ return;
+ }
+ }
+
+ /* SPLIT variables. */
+ size_t ss = 0;
+ size_t sn = var_get_dict_index (rowtype_);
+ struct ccase *sc = NULL;
+ size_t si = 0;
+
+ /* FACTOR variables. */
+ size_t fs = var_get_dict_index (rowtype_) + 1;
+ size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
+ struct ccase *fc = NULL;
+ size_t fi = 0;
+
+ /* Continuous variables. */
+ size_t cs = var_get_dict_index (varname_) + 1;
+ size_t cn = dict_get_n_vars (d) - cs;
+ struct ccase *cc = NULL;
+
+ /* Pivot table. */
+ struct pivot_table *pt = pivot_table_create (
+ N_("Matrix Variables Created by MGET"));
+ struct pivot_dimension *attr_dimension = pivot_dimension_create (
+ pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
+ struct pivot_dimension *var_dimension = pivot_dimension_create (
+ pt, PIVOT_AXIS_ROW, N_("Variable"));
+ if (sn > 0)
+ {
+ struct pivot_category *splits = pivot_category_create_group (
+ attr_dimension->root, N_("Split Values"));
+ for (size_t i = 0; i < sn; i++)
+ pivot_category_create_leaf (splits, pivot_value_new_variable (
+ dict_get_var (d, ss + i)));
+ }
+ if (fn > 0)
+ {
+ struct pivot_category *factors = pivot_category_create_group (
+ attr_dimension->root, N_("Factors"));
+ for (size_t i = 0; i < fn; i++)
+ pivot_category_create_leaf (factors, pivot_value_new_variable (
+ dict_get_var (d, fs + i)));
+ }
+ pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
+ N_("Rows"), N_("Columns"));
+
+ /* Matrix. */
+ struct ccase **rows = NULL;
+ size_t allocated_rows = 0;
+ size_t n_rows = 0;
+
+ struct ccase *c;
+ while ((c = casereader_read (r)) != NULL)
+ {
+ bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
+
+ enum
+ {
+ SPLITS_CHANGED,
+ FACTORS_CHANGED,
+ ROWTYPE_CHANGED,
+ NOTHING_CHANGED
+ }
+ change
+ = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
+ : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
+ : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
+ : NOTHING_CHANGED);
+
+ if (change != NOTHING_CHANGED)
+ {
+ matrix_mget_commit_var (rows, n_rows, d,
+ rowtype_, &mget->rowtypes,
+ mget->state,
+ ss, sn, si,
+ fs, fn, fi,
+ cs, cn,
+ pt, var_dimension);
+ n_rows = 0;
+ case_unref (cc);
+ cc = case_ref (c);
+ }
+
+ if (n_rows >= allocated_rows)
+ rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
+ rows[n_rows++] = c;
+
+ if (change == SPLITS_CHANGED)
+ {
+ si++;
+ case_unref (sc);
+ sc = case_ref (c);
+
+ /* Reset the factor number, if there are factors. */
+ if (fn)
+ {
+ fi = 0;
+ if (row_has_factors)
+ fi++;
+ case_unref (fc);
+ fc = case_ref (c);
+ }
+ }
+ else if (change == FACTORS_CHANGED)
+ {
+ if (row_has_factors)
+ fi++;
+ case_unref (fc);
+ fc = case_ref (c);
+ }
+ }
+ matrix_mget_commit_var (rows, n_rows, d,
+ rowtype_, &mget->rowtypes,
+ mget->state,
+ ss, sn, si,
+ fs, fn, fi,
+ cs, cn,
+ pt, var_dimension);
+ free (rows);
+
+ case_unref (sc);
+ case_unref (fc);
+ case_unref (cc);
+
+ if (var_dimension->n_leaves)
+ pivot_table_submit (pt);
+ else
+ pivot_table_unref (pt);
+}
+
+static void
+matrix_mget_execute (struct matrix_command *cmd)
+{
+ struct matrix_mget *mget = &cmd->mget;
+ struct casereader *r;
+ struct dictionary *d;
+ if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
+ mget->state->dataset, &r, &d))
+ {
+ matrix_mget_execute__ (cmd, r, d);
+ matrix_close_casereader (mget->file, mget->state->dataset, r, d);
+ }
+}
+\f
+/* CALL EIGEN. */
+
+static bool
+matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
+{
+ if (!lex_force_id (s->lexer))
+ return false;
+
+ *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (!*varp)
+ *varp = matrix_var_create (s, lex_tokss (s->lexer));
+ lex_get (s->lexer);
+ return true;
+}
+
+static struct matrix_command *
+matrix_eigen_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_EIGEN,
+ .eigen = { .expr = NULL }
+ };
+
+ struct matrix_eigen *eigen = &cmd->eigen;
+ if (!lex_force_match (s->lexer, T_LPAREN))
+ goto error;
+ eigen->expr = matrix_expr_parse (s);
+ if (!eigen->expr
+ || !lex_force_match (s->lexer, T_COMMA)
+ || !matrix_parse_dst_var (s, &eigen->evec)
+ || !lex_force_match (s->lexer, T_COMMA)
+ || !matrix_parse_dst_var (s, &eigen->eval)
+ || !lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_eigen_execute (struct matrix_command *cmd)
+{
+ struct matrix_eigen *eigen = &cmd->eigen;
+ gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
+ if (!A)
+ return;
+ if (!is_symmetric (A))
+ {
+ msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
+ gsl_matrix_free (A);
+ return;
+ }
+
+ gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
+ gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
+ gsl_vector v_eval = to_vector (eval);
+ gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
+ gsl_eigen_symmv (A, &v_eval, evec, w);
+ gsl_eigen_symmv_free (w);
+
+ gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
+
+ gsl_matrix_free (A);
+
+ gsl_matrix_free (eigen->eval->value);
+ eigen->eval->value = eval;
+
+ gsl_matrix_free (eigen->evec->value);
+ eigen->evec->value = evec;
+}
+\f
+/* CALL SETDIAG. */
+
+static struct matrix_command *
+matrix_setdiag_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_SETDIAG,
+ .setdiag = { .dst = NULL }
+ };
+
+ struct matrix_setdiag *setdiag = &cmd->setdiag;
+ if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
+ goto error;
+
+ setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (!setdiag->dst)
+ {
+ lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
+ goto error;
+ }
+ lex_get (s->lexer);
+
+ if (!lex_force_match (s->lexer, T_COMMA))
+ goto error;
+
+ setdiag->expr = matrix_expr_parse (s);
+ if (!setdiag->expr)
+ goto error;
+
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_setdiag_execute (struct matrix_command *cmd)
+{
+ struct matrix_setdiag *setdiag = &cmd->setdiag;
+ gsl_matrix *dst = setdiag->dst->value;
+ if (!dst)
+ {
+ msg_at (SE, cmd->location,
+ _("SETDIAG destination matrix %s is uninitialized."),
+ setdiag->dst->name);
+ return;
+ }
+
+ gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
+ if (!src)
+ return;
+
+ size_t n = MIN (dst->size1, dst->size2);
+ if (is_scalar (src))
+ {
+ double d = to_scalar (src);
+ for (size_t i = 0; i < n; i++)
+ gsl_matrix_set (dst, i, i, d);
+ }
+ else if (is_vector (src))
+ {
+ gsl_vector v = to_vector (src);
+ for (size_t i = 0; i < n && i < v.size; i++)
+ gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
+ }
+ else
+ msg_at (SE, matrix_expr_location (setdiag->expr),
+ _("SETDIAG argument 2 must be a scalar or a vector, "
+ "not a %zu×%zu matrix."),
+ src->size1, src->size2);
+ gsl_matrix_free (src);
+}
+\f
+/* CALL SVD. */
+
+static struct matrix_command *
+matrix_svd_parse (struct matrix_state *s)
+{
+ struct matrix_command *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_command) {
+ .type = MCMD_SVD,
+ .svd = { .expr = NULL }
+ };
+
+ struct matrix_svd *svd = &cmd->svd;
+ if (!lex_force_match (s->lexer, T_LPAREN))
+ goto error;
+ svd->expr = matrix_expr_parse (s);
+ if (!svd->expr
+ || !lex_force_match (s->lexer, T_COMMA)
+ || !matrix_parse_dst_var (s, &svd->u)
+ || !lex_force_match (s->lexer, T_COMMA)
+ || !matrix_parse_dst_var (s, &svd->s)
+ || !lex_force_match (s->lexer, T_COMMA)
+ || !matrix_parse_dst_var (s, &svd->v)
+ || !lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_command_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_svd_execute (struct matrix_svd *svd)
+{
+ gsl_matrix *m = matrix_expr_evaluate (svd->expr);
+ if (!m)
+ return;
+
+ if (m->size1 >= m->size2)
+ {
+ gsl_matrix *A = m;
+ gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
+ gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
+ gsl_vector Sv = gsl_matrix_diagonal (S).vector;
+ gsl_vector *work = gsl_vector_alloc (A->size2);
+ gsl_linalg_SV_decomp (A, V, &Sv, work);
+ gsl_vector_free (work);
+
+ matrix_var_set (svd->u, A);
+ matrix_var_set (svd->s, S);
+ matrix_var_set (svd->v, V);
+ }
+ else
+ {
+ gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
+ gsl_matrix_transpose_memcpy (At, m);
+ gsl_matrix_free (m);
+
+ gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
+ gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
+ gsl_vector Stv = gsl_matrix_diagonal (St).vector;
+ gsl_vector *work = gsl_vector_alloc (At->size2);
+ gsl_linalg_SV_decomp (At, Vt, &Stv, work);
+ gsl_vector_free (work);
+
+ matrix_var_set (svd->v, At);
+ matrix_var_set (svd->s, St);
+ matrix_var_set (svd->u, Vt);
+ }
+}
+\f
+/* The main MATRIX command logic. */
+
+static bool
+matrix_command_execute (struct matrix_command *cmd)
+{
+ switch (cmd->type)
+ {
+ case MCMD_COMPUTE:
+ matrix_compute_execute (cmd);
+ break;
+
+ case MCMD_PRINT:
+ matrix_print_execute (&cmd->print);
+ break;
+
+ case MCMD_DO_IF:
+ return matrix_do_if_execute (&cmd->do_if);
+
+ case MCMD_LOOP:
+ matrix_loop_execute (&cmd->loop);
+ break;
+
+ case MCMD_BREAK:
+ return false;
+
+ case MCMD_DISPLAY:
+ matrix_display_execute (&cmd->display);
+ break;
+
+ case MCMD_RELEASE:
+ matrix_release_execute (&cmd->release);
+ break;
+
+ case MCMD_SAVE:
+ matrix_save_execute (cmd);
+ break;
+
+ case MCMD_READ:
+ matrix_read_execute (cmd);
+ break;
+
+ case MCMD_WRITE:
+ matrix_write_execute (&cmd->write);
+ break;
+
+ case MCMD_GET:
+ matrix_get_execute (cmd);
+ break;
+
+ case MCMD_MSAVE:
+ matrix_msave_execute (cmd);
+ break;
+
+ case MCMD_MGET:
+ matrix_mget_execute (cmd);
+ break;
+
+ case MCMD_EIGEN:
+ matrix_eigen_execute (cmd);
+ break;
+
+ case MCMD_SETDIAG:
+ matrix_setdiag_execute (cmd);
+ break;
+
+ case MCMD_SVD:
+ matrix_svd_execute (&cmd->svd);
+ break;
+ }
+
+ return true;
+}
+
+static void
+matrix_command_destroy (struct matrix_command *cmd)
+{
+ if (!cmd)
+ return;
+
+ msg_location_destroy (cmd->location);
+
+ switch (cmd->type)
+ {
+ case MCMD_COMPUTE:
+ matrix_lvalue_destroy (cmd->compute.lvalue);
+ matrix_expr_destroy (cmd->compute.rvalue);
+ break;
+
+ case MCMD_PRINT:
+ matrix_expr_destroy (cmd->print.expression);
+ free (cmd->print.title);
+ print_labels_destroy (cmd->print.rlabels);
+ print_labels_destroy (cmd->print.clabels);
+ break;
+
+ case MCMD_DO_IF:
+ for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
+ {
+ matrix_expr_destroy (cmd->do_if.clauses[i].condition);
+ matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
+ }
+ free (cmd->do_if.clauses);
+ break;
+
+ case MCMD_LOOP:
+ matrix_expr_destroy (cmd->loop.start);
+ matrix_expr_destroy (cmd->loop.end);
+ matrix_expr_destroy (cmd->loop.increment);
+ matrix_expr_destroy (cmd->loop.top_condition);
+ matrix_expr_destroy (cmd->loop.bottom_condition);
+ matrix_commands_uninit (&cmd->loop.commands);
+ break;
+
+ case MCMD_BREAK:
+ break;
+
+ case MCMD_DISPLAY:
+ break;
+
+ case MCMD_RELEASE:
+ free (cmd->release.vars);
+ break;
+
+ case MCMD_SAVE:
+ matrix_expr_destroy (cmd->save.expression);
+ break;
+
+ case MCMD_READ:
+ matrix_lvalue_destroy (cmd->read.dst);
+ matrix_expr_destroy (cmd->read.size);
+ break;
+
+ case MCMD_WRITE:
+ matrix_expr_destroy (cmd->write.expression);
+ free (cmd->write.format);
+ break;
+
+ case MCMD_GET:
+ matrix_lvalue_destroy (cmd->get.dst);
+ fh_unref (cmd->get.file);
+ free (cmd->get.encoding);
+ var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
+ break;
+
+ case MCMD_MSAVE:
+ matrix_expr_destroy (cmd->msave.expr);
+ break;
+
+ case MCMD_MGET:
+ fh_unref (cmd->mget.file);
+ stringi_set_destroy (&cmd->mget.rowtypes);
+ break;
+
+ case MCMD_EIGEN:
+ matrix_expr_destroy (cmd->eigen.expr);
+ break;
+
+ case MCMD_SETDIAG:
+ matrix_expr_destroy (cmd->setdiag.expr);
+ break;
+
+ case MCMD_SVD:
+ matrix_expr_destroy (cmd->svd.expr);
+ break;
+ }
+ free (cmd);
+}
+
+static bool
+matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
+ const char *command_name,
+ const char *stop1, const char *stop2)
+{
+ lex_end_of_command (s->lexer);
+ lex_discard_rest_of_command (s->lexer);
+
+ size_t allocated = 0;
+ for (;;)
+ {
+ while (lex_token (s->lexer) == T_ENDCMD)
+ lex_get (s->lexer);
+
+ if (lex_at_phrase (s->lexer, stop1)
+ || (stop2 && lex_at_phrase (s->lexer, stop2)))
+ return true;
+
+ if (lex_at_phrase (s->lexer, "END MATRIX"))
+ {
+ lex_next_error (s->lexer, 0, 1,
+ _("Premature END MATRIX within %s."), command_name);
+ return false;
+ }
+
+ struct matrix_command *cmd = matrix_command_parse (s);
+ if (!cmd)
+ return false;
+
+ if (c->n >= allocated)
+ c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
+ c->commands[c->n++] = cmd;
+ }
+}
+
+static void
+matrix_commands_uninit (struct matrix_commands *cmds)
+{
+ for (size_t i = 0; i < cmds->n; i++)
+ matrix_command_destroy (cmds->commands[i]);
+ free (cmds->commands);
+}
+
+struct matrix_command_name
+ {
+ const char *name;
+ struct matrix_command *(*parse) (struct matrix_state *);
+ };
+
+static const struct matrix_command_name *
+matrix_command_name_parse (struct lexer *lexer)
+{
+ static const struct matrix_command_name commands[] = {
+ { "COMPUTE", matrix_compute_parse },
+ { "CALL EIGEN", matrix_eigen_parse },
+ { "CALL SETDIAG", matrix_setdiag_parse },
+ { "CALL SVD", matrix_svd_parse },
+ { "PRINT", matrix_print_parse },
+ { "DO IF", matrix_do_if_parse },
+ { "LOOP", matrix_loop_parse },
+ { "BREAK", matrix_break_parse },
+ { "READ", matrix_read_parse },
+ { "WRITE", matrix_write_parse },
+ { "GET", matrix_get_parse },
+ { "SAVE", matrix_save_parse },
+ { "MGET", matrix_mget_parse },
+ { "MSAVE", matrix_msave_parse },
+ { "DISPLAY", matrix_display_parse },
+ { "RELEASE", matrix_release_parse },
+ };
+ static size_t n = sizeof commands / sizeof *commands;
+
+ for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
+ if (lex_match_phrase (lexer, c->name))
+ return c;
+ return NULL;
+}
+
+static struct matrix_command *
+matrix_command_parse (struct matrix_state *s)
+{
+ int start_ofs = lex_ofs (s->lexer);
+ size_t nesting_level = SIZE_MAX;
+
+ struct matrix_command *c = NULL;
+ const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
+ if (!cmd)
+ lex_error (s->lexer, _("Unknown matrix command."));
+ else if (!cmd->parse)
+ lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
+ cmd->name);
+ else
+ {
+ nesting_level = output_open_group (
+ group_item_create_nocopy (utf8_to_title (cmd->name),
+ utf8_to_title (cmd->name)));
+ c = cmd->parse (s);
+ }
+
+ if (c)
+ {
+ c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
+ msg_location_remove_columns (c->location);
+ lex_end_of_command (s->lexer);
+ }
+ lex_discard_rest_of_command (s->lexer);
+ if (nesting_level != SIZE_MAX)
+ output_close_groups (nesting_level);
+
+ return c;
+}
+
+int
+cmd_matrix (struct lexer *lexer, struct dataset *ds)
+{
+ if (!lex_force_match (lexer, T_ENDCMD))
+ return CMD_FAILURE;
+
+ struct matrix_state state = {
+ .dataset = ds,
+ .session = dataset_session (ds),
+ .lexer = lexer,
+ .vars = HMAP_INITIALIZER (state.vars),
+ };
+
+ for (;;)
+ {
+ while (lex_match (lexer, T_ENDCMD))
+ continue;
+ if (lex_token (lexer) == T_STOP)
+ {
+ msg (SE, _("Unexpected end of input expecting matrix command."));
+ break;
+ }
+
+ if (lex_match_phrase (lexer, "END MATRIX"))
+ break;
+
+ struct matrix_command *c = matrix_command_parse (&state);
+ if (c)
+ {
+ matrix_command_execute (c);
+ matrix_command_destroy (c);
+ }
+ }
+
+ struct matrix_var *var, *next;
+ HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
+ {
+ free (var->name);
+ gsl_matrix_free (var->value);
+ hmap_delete (&state.vars, &var->hmap_node);
+ free (var);
+ }
+ hmap_destroy (&state.vars);
+ msave_common_destroy (state.msave_common);
+ fh_unref (state.prev_read_file);
+ for (size_t i = 0; i < state.n_read_files; i++)
+ read_file_destroy (state.read_files[i]);
+ free (state.read_files);
+ fh_unref (state.prev_write_file);
+ for (size_t i = 0; i < state.n_write_files; i++)
+ write_file_destroy (state.write_files[i]);
+ free (state.write_files);
+ fh_unref (state.prev_save_file);
+ for (size_t i = 0; i < state.n_save_files; i++)
+ save_file_destroy (state.save_files[i]);
+ free (state.save_files);
+
+ return CMD_SUCCESS;
+}
projects / pspp / blobdiff