+/* 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_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/data-io/data-reader.h"
+#include "language/data-io/data-writer.h"
+#include "language/data-io/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/misc.h"
+#include "libpspp/str.h"
+#include "libpspp/string-array.h"
+#include "libpspp/stringi-set.h"
+#include "libpspp/u8-line.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/minmax.h"
+#include "gl/xsize.h"
+
+#include "gettext.h"
+#define _(msgid) gettext (msgid)
+#define N_(msgid) (msgid)
+
+struct matrix_var
+ {
+ struct hmap_node hmap_node;
+ const char *name;
+ gsl_matrix *value;
+ };
+
+struct msave_common
+ {
+ /* Configuration. */
+ struct file_handle *outfile;
+ struct string_array variables;
+ struct string_array fnames;
+ struct string_array snames;
+ bool has_factors;
+ bool has_splits;
+ size_t n_varnames;
+
+ /* Execution state. */
+ struct dictionary *dict;
+ struct casewriter *writer;
+ };
+
+struct read_file
+ {
+ struct file_handle *file;
+ struct dfm_reader *reader;
+ char *encoding;
+ };
+
+struct matrix_state
+ {
+ struct dataset *dataset;
+ struct session *session;
+ struct lexer *lexer;
+ struct hmap vars;
+ bool in_loop;
+ struct file_handle *prev_save_outfile;
+ struct file_handle *prev_write_outfile;
+ struct msave_common *common;
+
+ struct file_handle *prev_read_file;
+ struct read_file **read_files;
+ size_t n_read_files;
+ };
+
+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;
+}
+
+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;
+}
+
+static void
+matrix_var_set (struct matrix_var *var, gsl_matrix *value)
+{
+ gsl_matrix_free (var->value);
+ var->value = value;
+}
+\f
+#define MATRIX_FUNCTIONS \
+ F(ABS, m_m) \
+ F(ALL, d_m) \
+ F(ANY, d_m) \
+ F(ARSIN, m_m) \
+ F(ARTAN, m_m) \
+ F(BLOCK, m_any) \
+ F(CHOL, m_m) \
+ F(CMIN, m_m) \
+ F(CMAX, m_m) \
+ F(COS, m_m) \
+ F(CSSQ, m_m) \
+ F(CSUM, m_m) \
+ F(DESIGN, m_m) \
+ F(DET, d_m) \
+ F(DIAG, m_m) \
+ F(EVAL, m_m) \
+ F(EXP, m_m) \
+ F(GINV, m_m) \
+ F(GRADE, m_m) \
+ F(GSCH, m_m) \
+ F(IDENT, IDENT) \
+ F(INV, m_m) \
+ F(KRONEKER, m_mm) \
+ F(LG10, m_m) \
+ F(LN, m_m) \
+ F(MAGIC, m_d) \
+ F(MAKE, m_ddd) \
+ F(MDIAG, m_v) \
+ F(MMAX, d_m) \
+ F(MMIN, d_m) \
+ F(MOD, m_md) \
+ F(MSSQ, d_m) \
+ F(MSUM, d_m) \
+ F(NCOL, d_m) \
+ F(NROW, d_m) \
+ F(RANK, d_m) \
+ F(RESHAPE, m_mdd) \
+ F(RMAX, m_m) \
+ F(RMIN, m_m) \
+ F(RND, m_m) \
+ F(RNKORDER, m_m) \
+ F(RSSQ, m_m) \
+ F(RSUM, m_m) \
+ F(SIN, m_m) \
+ F(SOLVE, m_mm) \
+ F(SQRT, m_m) \
+ F(SSCP, m_m) \
+ F(SVAL, m_m) \
+ F(SWEEP, m_md) \
+ F(T, m_m) \
+ F(TRACE, d_m) \
+ F(TRANSPOS, m_m) \
+ F(TRUNC, m_m) \
+ F(UNIFORM, m_dd)
+
+enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
+enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
+enum { m_ddd_MIN_ARGS = 3, m_ddd_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_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
+enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
+enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
+enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
+enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
+enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
+
+typedef gsl_matrix *matrix_proto_m_d (double);
+typedef gsl_matrix *matrix_proto_m_dd (double, double);
+typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
+typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
+typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
+typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
+typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
+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(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
+MATRIX_FUNCTIONS
+#undef F
+\f
+/* Matrix expressions. */
+
+struct matrix_expr
+ {
+ enum matrix_op
+ {
+ /* Functions. */
+#define F(NAME, PROTOTYPE) MOP_F_##NAME,
+ 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
+ {
+ struct
+ {
+ struct matrix_expr **subs;
+ size_t n_subs;
+ };
+
+ double number;
+ struct matrix_var *variable;
+ struct read_file *eof;
+ };
+ };
+
+static void
+matrix_expr_destroy (struct matrix_expr *e)
+{
+ if (!e)
+ return;
+
+ switch (e->op)
+ {
+#define F(NAME, PROTOTYPE) case MOP_F_##NAME:
+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]);
+ break;
+
+ case MOP_NUMBER:
+ case MOP_VARIABLE:
+ case MOP_EOF:
+ break;
+ }
+ free (e);
+}
+
+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);
+}
+
+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_parse_or_xor (struct matrix_state *);
+
+static struct matrix_expr *
+matrix_parse_curly_comma (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_or_xor (s);
+ if (!lhs)
+ return NULL;
+
+ while (lex_match (s->lexer, T_COMMA))
+ {
+ struct matrix_expr *rhs = matrix_parse_or_xor (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
+ }
+ return lhs;
+}
+
+static struct matrix_expr *
+matrix_parse_curly_semi (struct matrix_state *s)
+{
+ if (lex_token (s->lexer) == T_RCURLY)
+ return matrix_expr_create_0 (MOP_EMPTY);
+
+ struct matrix_expr *lhs = matrix_parse_curly_comma (s);
+ if (!lhs)
+ return NULL;
+
+ while (lex_match (s->lexer, T_SEMICOLON))
+ {
+ struct matrix_expr *rhs = matrix_parse_curly_comma (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
+ }
+ return lhs;
+}
+
+#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 (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 gsl_matrix *
+matrix_eval_ABS (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = fabs (*d);
+ return m;
+}
+
+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 !.0;
+ return 0.0;
+}
+
+static gsl_matrix *
+matrix_eval_ARSIN (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = asin (*d);
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_ARTAN (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = atan (*d);
+ return m;
+}
+
+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)
+{
+ 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 (SE, _("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 gsl_matrix *
+matrix_eval_COS (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = cos (*d);
+ return m;
+}
+
+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)
+{
+ 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 (MW, _("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)
+{
+ if (!is_symmetric (m))
+ {
+ msg (SE, _("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 gsl_matrix *
+matrix_eval_EXP (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = exp (*d);
+ return m;
+}
+
+/* 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 */
+ tmp_mat = gsl_matrix_alloc (m, n);
+ gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat);
+
+ gsl_matrix *A_pinv;
+ if (swap)
+ {
+ A_pinv = gsl_matrix_alloc (n, m);
+ gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat, 0., A_pinv);
+ }
+ else
+ {
+ A_pinv = gsl_matrix_alloc (m, n);
+ gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat, U, 0., A_pinv);
+ }
+
+ gsl_matrix_free (tmp_mat);
+ 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)
+{
+ if (v->size2 < v->size1)
+ {
+ msg (SE, _("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 (SE, _("Argument to GSCH with dimensions (%zu,%zu) contains only "
+ "%zu linearly independent columns."),
+ v->size1, v->size2, ux);
+ return NULL;
+ }
+
+ u->size2 = v->size1;
+ return u;
+}
+
+static gsl_matrix *
+matrix_eval_IDENT (double s1, double s2)
+{
+ if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
+ {
+ msg (SE, _("Arguments to IDENT must be non-negative integers."));
+ return NULL;
+ }
+
+ gsl_matrix *m = gsl_matrix_alloc (s1, s2);
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = x == y;
+ return m;
+}
+
+static void
+invert_matrix (gsl_matrix *x)
+{
+ gsl_permutation *p = gsl_permutation_alloc (x->size1);
+ int signum;
+ gsl_linalg_LU_decomp (x, p, &signum);
+ gsl_linalg_LU_invx (x, p);
+ gsl_permutation_free (p);
+}
+
+static gsl_matrix *
+matrix_eval_INV (gsl_matrix *m)
+{
+ invert_matrix (m);
+ return m;
+}
+
+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 gsl_matrix *
+matrix_eval_LG10 (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = log10 (*d);
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_LN (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = log (*d);
+ return m;
+}
+
+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_)
+{
+ if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
+ {
+ msg (SE, _("MAGIC argument must be an integer 3 or greater."));
+ return NULL;
+ }
+ 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)
+{
+ if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
+ {
+ msg (SE, _("Size arguments to MAKE must be integers."));
+ return NULL;
+ }
+
+ 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_alloc (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)
+{
+ if (divisor == 0.0)
+ {
+ msg (SE, _("Divisor argument to MOD function must be nonzero."));
+ return NULL;
+ }
+
+ 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_)
+{
+ if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
+ {
+ msg (SE, _("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)
+ {
+ msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
+ "product of matrix dimensions."));
+ 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 gsl_matrix *
+matrix_eval_RND (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = rint (*d);
+ return m;
+}
+
+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 gsl_matrix *
+matrix_eval_SIN (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = sin (*d);
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
+{
+ if (m1->size1 != m2->size1)
+ {
+ msg (SE, _("SOLVE requires its arguments to have the same number of "
+ "rows, but the first argument has dimensions (%zu,%zu) and "
+ "the second (%zu,%zu)."),
+ m1->size1, m1->size2,
+ m2->size1, m2->size2);
+ 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 gsl_matrix *
+matrix_eval_SQRT (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ {
+ if (*d < 0)
+ {
+ msg (SE, _("Argument to SQRT must be nonnegative."));
+ return NULL;
+ }
+ *d = sqrt (*d);
+ }
+ return m;
+}
+
+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)
+{
+ if (d < 1 || d > SIZE_MAX)
+ {
+ msg (SE, _("Scalar argument to SWEEP must be integer."));
+ return NULL;
+ }
+ size_t k = d - 1;
+ if (k >= MIN (m->size1, m->size2))
+ {
+ msg (SE, _("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 gsl_matrix *
+matrix_eval_TRUNC (gsl_matrix *m)
+{
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
+ *d = trunc (*d);
+ return m;
+}
+
+static gsl_matrix *
+matrix_eval_UNIFORM (double r_, double c_)
+{
+ if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
+ {
+ msg (SE, _("Arguments to UNIFORM must be integers."));
+ return NULL;
+ }
+ size_t r = r_;
+ size_t c = c_;
+ if (size_overflow_p (xtimes (r, xmax (c, 1))))
+ {
+ msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
+ 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;
+}
+
+struct matrix_function
+ {
+ const char *name;
+ enum matrix_op op;
+ size_t min_args, max_args;
+ };
+
+static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
+
+static const struct matrix_function *
+matrix_parse_function_name (const char *token)
+{
+ static const struct matrix_function functions[] =
+ {
+#define F(NAME, PROTO) { #NAME, MOP_F_##NAME, 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 (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
+ return &functions[i];
+ }
+ return NULL;
+}
+
+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 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;
+
+ 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 };
+ *exprp = e;
+ return true;
+ }
+
+ const struct matrix_function *f = matrix_parse_function_name (token);
+ if (!f)
+ return false;
+
+ lex_get_n (s->lexer, 2);
+
+ struct matrix_expr *e = xmalloc (sizeof *e);
+ *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
+
+ size_t allocated_subs = 0;
+ do
+ {
+ struct matrix_expr *sub = matrix_parse_expr (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;
+
+ *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_parse_or_xor (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, NULL);
+ return NULL;
+}
+
+static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
+
+static bool
+matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
+{
+ if (lex_match (s->lexer, T_COLON))
+ {
+ *indexp = NULL;
+ return true;
+ }
+ else
+ {
+ *indexp = matrix_parse_expr (s);
+ 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))
+ {
+ 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)
+ || !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)
+{
+ if (lex_match (s->lexer, T_DASH))
+ {
+ struct matrix_expr *lhs = matrix_parse_unary (s);
+ return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
+ }
+ else if (lex_match (s->lexer, T_PLUS))
+ return matrix_parse_unary (s);
+ else
+ return matrix_parse_postfix (s);
+}
+
+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)
+{
+ struct matrix_expr *lhs = matrix_parse_seq (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (lex_match (s->lexer, T_EXP))
+ op = MOP_EXP_MAT;
+ else if (lex_match_phrase (s->lexer, "&**"))
+ op = MOP_EXP_ELEMS;
+ else
+ return lhs;
+
+ struct matrix_expr *rhs = matrix_parse_seq (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_mul_div (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_exp (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (lex_match (s->lexer, T_ASTERISK))
+ op = MOP_MUL_MAT;
+ else if (lex_match (s->lexer, T_SLASH))
+ op = MOP_DIV_ELEMS;
+ else if (lex_match_phrase (s->lexer, "&*"))
+ op = MOP_MUL_ELEMS;
+ else if (lex_match_phrase (s->lexer, "&/"))
+ op = MOP_DIV_ELEMS;
+ else
+ return lhs;
+
+ struct matrix_expr *rhs = matrix_parse_exp (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+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)
+{
+ struct matrix_expr *lhs = matrix_parse_add_sub (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (lex_match (s->lexer, T_GT))
+ op = MOP_GT;
+ else if (lex_match (s->lexer, T_GE))
+ op = MOP_GE;
+ else if (lex_match (s->lexer, T_LT))
+ op = MOP_LT;
+ else if (lex_match (s->lexer, T_LE))
+ op = MOP_LE;
+ else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
+ op = MOP_EQ;
+ else if (lex_match (s->lexer, T_NE))
+ op = MOP_NE;
+ else
+ return lhs;
+
+ struct matrix_expr *rhs = matrix_parse_add_sub (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_not (struct matrix_state *s)
+{
+ if (lex_match (s->lexer, T_NOT))
+ {
+ struct matrix_expr *lhs = matrix_parse_not (s);
+ return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
+ }
+ else
+ return matrix_parse_relational (s);
+}
+
+static struct matrix_expr *
+matrix_parse_and (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_not (s);
+ if (!lhs)
+ return NULL;
+
+ while (lex_match (s->lexer, T_AND))
+ {
+ struct matrix_expr *rhs = matrix_parse_not (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
+ }
+ return lhs;
+}
+
+static struct matrix_expr *
+matrix_parse_or_xor (struct matrix_state *s)
+{
+ struct matrix_expr *lhs = matrix_parse_and (s);
+ if (!lhs)
+ return NULL;
+
+ for (;;)
+ {
+ enum matrix_op op;
+ if (lex_match (s->lexer, T_OR))
+ op = MOP_OR;
+ else if (lex_match_id (s->lexer, "XOR"))
+ op = MOP_XOR;
+ else
+ return lhs;
+
+ struct matrix_expr *rhs = matrix_parse_and (s);
+ if (!rhs)
+ {
+ matrix_expr_destroy (lhs);
+ return NULL;
+ }
+ lhs = matrix_expr_create_2 (op, lhs, rhs);
+ }
+}
+
+static struct matrix_expr *
+matrix_parse_expr (struct matrix_state *s)
+{
+ return matrix_parse_or_xor (s);
+}
+\f
+/* Expression evaluation. */
+
+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(NAME, PROTOTYPE) case MOP_F_##NAME:
+ 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(NAME, PROTOTYPE) case MOP_F_##NAME:
+ 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 (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 (SE, _("Operands to %s must have the same dimensions or one "
+ "must be a scalar, not matrices with dimensions (%zu,%zu) "
+ "and (%zu,%zu)."),
+ matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
+ return NULL;
+ }
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
+{
+ if (is_scalar (a) || is_scalar (b))
+ return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
+
+ if (a->size2 != b->size1)
+ {
+ msg (SE, _("Matrices with dimensions (%zu,%zu) and (%zu,%zu) are "
+ "not conformable for multiplication."),
+ a->size1, a->size2, 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 (gsl_matrix *x_, gsl_matrix *b)
+{
+ gsl_matrix *x = x_;
+ if (x->size1 != x->size2)
+ {
+ msg (SE, _("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 (SE, _("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 (SE, _("Exponent %.1f in matrix multiplication is non-integer "
+ "or outside the valid range."), bf);
+ return NULL;
+ }
+ long int bl = bf;
+
+ gsl_matrix *tmp = gsl_matrix_alloc (x->size1, x->size2);
+ gsl_matrix *y = gsl_matrix_alloc (x->size1, x->size2);
+ gsl_matrix_set_identity (y);
+ if (bl == 0)
+ return y;
+
+ for (unsigned long int n = labs (bl); n > 1; n /= 2)
+ if (n & 1)
+ {
+ mul_matrix (&y, x, y, &tmp);
+ square_matrix (&x, &tmp);
+ }
+ else
+ square_matrix (&x, &tmp);
+
+ mul_matrix (&y, x, y, &tmp);
+ if (bf < 0)
+ invert_matrix (y);
+
+ if (tmp != x_)
+ gsl_matrix_free (tmp);
+ return y;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
+ gsl_matrix *by_)
+{
+ if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
+ {
+ msg (SE, _("All operands of : operator must be scalars."));
+ 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 (SE, _("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)
+{
+ 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 = !(*ae > 0);
+ }
+ return a;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_paste_horz (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 (SE, _("All columns in a matrix must have the same number of rows, "
+ "but this tries to paste matrices with %zu and %zu rows."),
+ a->size1, b->size1);
+ 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 (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 (SE, _("All rows in a matrix must have the same number of columns, "
+ "but this tries to stack matrices with %zu and %zu columns."),
+ a->size2, b->size2);
+ 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;
+ 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 bool
+matrix_normalize_index_vector (gsl_matrix *m, 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 (SE, _("Vector index must be scalar or vector, not a "
+ "matrix with dimensions (%zu,%zu)."),
+ m->size1, m->size2);
+ break;
+
+ case IV_ROW:
+ msg (SE, _("Matrix row index must be scalar or vector, not a "
+ "matrix with dimensions (%zu,%zu)."),
+ m->size1, m->size2);
+ break;
+
+ case IV_COLUMN:
+ msg (SE, _("Matrix column index must be scalar or vector, not a "
+ "matrix with dimensions (%zu,%zu)."),
+ m->size1, m->size2);
+ break;
+ }
+ return false;
+ }
+
+ gsl_vector v = to_vector (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 (SE, _("Index %g is out of range for vector "
+ "with %zu elements."), index, size);
+ break;
+
+ case IV_ROW:
+ msg (SE, _("%g is not a valid row index for a matrix "
+ "with dimensions (%zu,%zu)."),
+ index, size, other_size);
+ break;
+
+ case IV_COLUMN:
+ msg (SE, _("%g is not a valid column index for a matrix "
+ "with dimensions (%zu,%zu)."),
+ index, other_size, size);
+ break;
+ }
+
+ free (iv->indexes);
+ 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 (gsl_matrix *sm)
+{
+ if (!is_vector (sm))
+ {
+ msg (SE, _("Vector index operator must be applied to vector, "
+ "not a matrix with dimensions (%zu,%zu)."),
+ sm->size1, sm->size2);
+ return NULL;
+ }
+
+ return sm;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
+{
+ if (!matrix_expr_evaluate_vec_all (sm))
+ return NULL;
+
+ gsl_vector sv = to_vector (sm);
+ struct index_vector iv;
+ if (!matrix_normalize_index_vector (im, 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));
+ }
+ free (iv.indexes);
+
+ return dm;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
+ gsl_matrix *im1)
+{
+ struct index_vector iv0;
+ if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
+ return NULL;
+
+ struct index_vector iv1;
+ if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
+ &iv1))
+ {
+ free (iv0.indexes);
+ 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));
+ }
+ }
+ free (iv0.indexes);
+ free (iv1.indexes);
+ return dm;
+}
+
+#define F(NAME, PROTOTYPE) \
+ static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
+ const char *name, gsl_matrix *[], size_t, \
+ matrix_proto_##PROTOTYPE *);
+MATRIX_FUNCTIONS
+#undef F
+
+static bool
+check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
+{
+ if (!is_scalar (subs[index]))
+ {
+ msg (SE, _("Function %s argument %zu must be a scalar, "
+ "but it has dimensions (%zu,%zu)."),
+ name, index + 1, subs[index]->size1, subs[index]->size2);
+ return false;
+ }
+ return true;
+}
+
+static bool
+check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
+{
+ if (!is_vector (subs[index]))
+ {
+ msg (SE, _("Function %s argument %zu must be a vector, "
+ "but it has dimensions (%zu,%zu)."),
+ name, index + 1, subs[index]->size1, subs[index]->size2);
+ return false;
+ }
+ return true;
+}
+
+static bool
+to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
+{
+ for (size_t i = 0; i < n_subs; i++)
+ {
+ if (!check_scalar_arg (name, subs, i))
+ return false;
+ d[i] = to_scalar (subs[i]);
+ }
+ return true;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_d (const char *name,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_d *f)
+{
+ assert (n_subs == 1);
+
+ double d;
+ return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_dd (const char *name,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_dd *f)
+{
+ assert (n_subs == 2);
+
+ double d[2];
+ return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_ddd (const char *name,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_ddd *f)
+{
+ assert (n_subs == 3);
+
+ double d[3];
+ return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_m (const char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_m *f)
+{
+ assert (n_subs == 1);
+ return f (subs[0]);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_md (const char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_md *f)
+{
+ assert (n_subs == 2);
+ if (!check_scalar_arg (name, subs, 1))
+ return NULL;
+ return f (subs[0], to_scalar (subs[1]));
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mdd (const char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_mdd *f)
+{
+ assert (n_subs == 3);
+ if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
+ return NULL;
+ return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_mm (const char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_mm *f)
+{
+ assert (n_subs == 2);
+ return f (subs[0], subs[1]);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_m_v (const char *name,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_v *f)
+{
+ assert (n_subs == 1);
+ if (!check_vector_arg (name, subs, 0))
+ return NULL;
+ gsl_vector v = to_vector (subs[0]);
+ return f (&v);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_d_m (const char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_d_m *f)
+{
+ assert (n_subs == 1);
+
+ 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 char *name UNUSED,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_any *f)
+{
+ return f (subs, n_subs);
+}
+
+static gsl_matrix *
+matrix_expr_evaluate_IDENT (const char *name,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_IDENT *f)
+{
+ assert (n_subs <= 2);
+
+ double d[2];
+ if (!to_scalar_args (name, subs, n_subs, d))
+ return NULL;
+
+ return f (d[0], d[n_subs - 1]);
+}
+
+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 (SE, _("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(NAME, PROTOTYPE) \
+ case MOP_F_##NAME: \
+ result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
+ subs, e->n_subs, \
+ matrix_eval_##NAME); \
+ 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->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 (subs[0], subs[1], NULL);
+ break;
+
+ case MOP_SEQ_BY:
+ result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
+ break;
+
+ case MOP_MUL_MAT:
+ result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
+ break;
+
+ case MOP_EXP_MAT:
+ result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
+ break;
+
+ case MOP_PASTE_HORZ:
+ result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
+ break;
+
+ case MOP_PASTE_VERT:
+ result = matrix_expr_evaluate_paste_vert (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 (subs[0], subs[1]);
+ break;
+
+ case MOP_VEC_ALL:
+ result = matrix_expr_evaluate_vec_all (subs[0]);
+ break;
+
+ case MOP_MAT_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
+ break;
+
+ case MOP_ROW_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
+ break;
+
+ case MOP_COL_INDEX:
+ result = matrix_expr_evaluate_mat_index (subs[0], NULL, 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 (SE, _("Expression for %s must evaluate to scalar, "
+ "not a matrix with dimensions (%zu,%zu)."),
+ 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 (SE, _("Expression for %s is outside the integer range."), context);
+ return false;
+ }
+ *integer = d;
+ return true;
+}
+\f
+struct matrix_lvalue
+ {
+ struct matrix_var *var;
+ struct matrix_expr *indexes[2];
+ size_t n_indexes;
+ };
+
+static void
+matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
+{
+ if (lvalue)
+ for (size_t i = 0; i < lvalue->n_indexes; i++)
+ matrix_expr_destroy (lvalue->indexes[i]);
+}
+
+static struct matrix_lvalue *
+matrix_lvalue_parse (struct matrix_state *s)
+{
+ struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
+ if (!lex_force_id (s->lexer))
+ goto error;
+ lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
+ if (lex_next_token (s->lexer, 1) == T_LPAREN)
+ {
+ if (!lvalue->var)
+ {
+ msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
+ return NULL;
+ }
+
+ lex_get_n (s->lexer, 2);
+
+ if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
+ goto error;
+ lvalue->n_indexes++;
+
+ if (lex_match (s->lexer, T_COMMA))
+ {
+ if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
+ goto error;
+ lvalue->n_indexes++;
+ }
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+ }
+ 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;
+
+ return matrix_normalize_index_vector (m, size, index_type, other_size, iv);
+}
+
+static bool
+matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
+ struct index_vector *iv, gsl_matrix *sm)
+{
+ gsl_vector dv = to_vector (lvalue->var->value);
+
+ /* Convert source matrix 'sm' to source vector 'sv'. */
+ if (!is_vector (sm))
+ {
+ msg (SE, _("Can't assign matrix with dimensions (%zu,%zu) to subvector."),
+ sm->size1, sm->size2);
+ return false;
+ }
+ gsl_vector sv = to_vector (sm);
+ if (iv->n != sv.size)
+ {
+ msg (SE, _("Can't assign vector with %zu elements "
+ "to subvector with %zu."), sv.size, iv->n);
+ return false;
+ }
+
+ /* Assign elements. */
+ 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)
+{
+ gsl_matrix *dm = lvalue->var->value;
+
+ /* Convert source matrix 'sm' to source vector 'sv'. */
+ if (iv0->n != sm->size1)
+ {
+ msg (SE, _("Row index vector for assignment to %s has %zu elements "
+ "but source matrix has %zu rows."),
+ lvalue->var->name, iv0->n, sm->size1);
+ return false;
+ }
+ if (iv1->n != sm->size2)
+ {
+ msg (SE, _("Column index vector for assignment to %s has %zu elements "
+ "but source matrix has %zu columns."),
+ lvalue->var->name, iv1->n, 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;
+}
+
+/* Takes ownership of M. */
+static bool
+matrix_lvalue_assign (struct matrix_lvalue *lvalue,
+ struct index_vector *iv0, struct index_vector *iv1,
+ gsl_matrix *sm)
+{
+ 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)
+ : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
+ free (iv0->indexes);
+ free (iv1->indexes);
+ gsl_matrix_free (sm);
+ return ok;
+ }
+}
+
+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 (SE, _("Undefined variable %s."), lvalue->var->name);
+ return false;
+ }
+ else if (dm->size1 == 0 || dm->size2 == 0)
+ {
+ msg (SE, _("Cannot index matrix with dimensions (%zu,%zu)."),
+ dm->size1, dm->size2);
+ return false;
+ }
+ else if (lvalue->n_indexes == 1)
+ {
+ if (!is_vector (dm))
+ {
+ msg (SE, _("Can't use vector indexing on matrix %s with "
+ "dimensions (%zu,%zu)."),
+ lvalue->var->name, dm->size1, dm->size2);
+ 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))
+ {
+ free (iv0->indexes);
+ return false;
+ }
+ return true;
+ }
+}
+
+/* Takes ownership of M. */
+static bool
+matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
+{
+ struct index_vector iv0, iv1;
+ return (matrix_lvalue_evaluate (lvalue, &iv0, &iv1)
+ && matrix_lvalue_assign (lvalue, &iv0, &iv1, sm));
+}
+
+\f
+/* Matrix command. */
+
+struct matrix_cmds
+ {
+ struct matrix_cmd **commands;
+ size_t n;
+ };
+
+struct matrix_cmd
+ {
+ enum matrix_cmd_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;
+
+ union
+ {
+ struct compute_command
+ {
+ struct matrix_lvalue *lvalue;
+ struct matrix_expr *rvalue;
+ }
+ compute;
+
+ struct print_command
+ {
+ 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 do_if_command
+ {
+ struct do_if_clause
+ {
+ struct matrix_expr *condition;
+ struct matrix_cmds commands;
+ }
+ *clauses;
+ size_t n_clauses;
+ }
+ do_if;
+
+ struct loop_command
+ {
+ /* 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_cmds commands;
+ }
+ loop;
+
+ struct save_command
+ {
+ struct matrix_expr *expression;
+ struct file_handle *outfile;
+ struct string_array *variables;
+ struct matrix_expr *names;
+ struct stringi_set strings;
+ }
+ save;
+
+ struct read_command
+ {
+ struct read_file *rf;
+ struct matrix_lvalue *dst;
+ struct matrix_expr *size;
+ int c1, c2;
+ enum fmt_type format;
+ int w;
+ int d;
+ bool symmetric;
+ bool reread;
+ }
+ read;
+
+ struct write_command
+ {
+ struct matrix_expr *expression;
+ struct file_handle *outfile;
+ char *encoding;
+ int c1, c2;
+ enum fmt_type format;
+ int w;
+ int d;
+ bool triangular;
+ bool hold;
+ }
+ write;
+
+ struct display_command
+ {
+ struct matrix_state *state;
+ bool dictionary;
+ bool status;
+ }
+ display;
+
+ struct release_command
+ {
+ struct matrix_var **vars;
+ size_t n_vars;
+ }
+ release;
+
+ struct get_command
+ {
+ struct matrix_lvalue *dst;
+ struct file_handle *file;
+ char *encoding;
+ struct string_array variables;
+ 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 msave_command
+ {
+ struct msave_common *common;
+ char *varname_;
+ struct matrix_expr *expr;
+ const char *rowtype;
+ struct matrix_expr *factors;
+ struct matrix_expr *splits;
+ }
+ msave;
+
+ struct mget_command
+ {
+ struct matrix_state *state;
+ struct file_handle *file;
+ struct stringi_set rowtypes;
+ }
+ mget;
+
+ struct eigen_command
+ {
+ struct matrix_expr *expr;
+ struct matrix_var *evec;
+ struct matrix_var *eval;
+ }
+ eigen;
+
+ struct setdiag_command
+ {
+ struct matrix_var *dst;
+ struct matrix_expr *expr;
+ }
+ setdiag;
+
+ struct svd_command
+ {
+ struct matrix_expr *expr;
+ struct matrix_var *u;
+ struct matrix_var *s;
+ struct matrix_var *v;
+ }
+ svd;
+ };
+ };
+
+static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
+static bool matrix_cmd_execute (struct matrix_cmd *);
+
+static void
+matrix_cmd_destroy (struct matrix_cmd *cmd)
+{
+ /* XXX */
+ free (cmd);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_compute (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_COMPUTE,
+ .compute = { .lvalue = NULL },
+ };
+
+ struct compute_command *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_parse_expr (s);
+ if (!compute->rvalue)
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_compute (struct compute_command *compute)
+{
+ gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
+ if (!value)
+ return;
+
+ matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
+}
+\f
+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_cmd *
+matrix_parse_print (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_PRINT,
+ .print = {
+ .use_default_format = true,
+ }
+ };
+
+ if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
+ {
+ size_t depth = 0;
+ for (size_t i = 0; ; i++)
+ {
+ enum token_type t = lex_next_token (s->lexer, i);
+ if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
+ depth++;
+ else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
+ depth--;
+ else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
+ {
+ if (i > 0)
+ cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
+ break;
+ }
+ }
+
+ cmd->print.expression = matrix_parse_exp (s);
+ if (!cmd->print.expression)
+ goto error;
+ }
+
+ 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_cmd_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_cmd_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_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
+ struct fmt_spec format, int log_scale)
+{
+ matrix_cmd_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);
+ string_array_destroy (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);
+}
+
+static void
+matrix_cmd_print_tables (const struct print_command *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_cmd_execute_print (const struct print_command *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_cmd_print_text (print, m, format, log_scale);
+ else
+ matrix_cmd_print_tables (print, m, format, log_scale);
+
+ gsl_matrix_free (m);
+ }
+ else
+ {
+ matrix_cmd_print_space (print->space);
+ if (print->title)
+ output_log ("%s", print->title);
+ }
+}
+\f
+/* DO IF. */
+
+static bool
+matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *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"))
+ {
+ msg (SE, _("Premature END MATRIX within %s."), command_name);
+ return false;
+ }
+
+ struct matrix_cmd *cmd = matrix_parse_command (s);
+ if (!cmd)
+ return false;
+
+ if (c->n >= allocated)
+ c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
+ c->commands[c->n++] = cmd;
+ }
+}
+
+static bool
+matrix_parse_do_if_clause (struct matrix_state *s,
+ struct do_if_command *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_parse_expr (s);
+ if (!c->condition)
+ return false;
+ }
+
+ return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
+}
+
+static struct matrix_cmd *
+matrix_parse_do_if (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_DO_IF,
+ .do_if = { .n_clauses = 0 }
+ };
+
+ size_t allocated_clauses = 0;
+ do
+ {
+ if (!matrix_parse_do_if_clause (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_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
+ goto error;
+
+ if (!lex_match_phrase (s->lexer, "END IF"))
+ NOT_REACHED ();
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static bool
+matrix_cmd_execute_do_if (struct do_if_command *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_cmd_execute (c->commands.commands[j]))
+ return false;
+ return true;
+ }
+ return true;
+}
+\f
+static struct matrix_cmd *
+matrix_parse_loop (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
+
+ struct loop_command *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_parse_expr (s);
+ if (!loop->start || !lex_force_match (s->lexer, T_TO))
+ goto error;
+
+ loop->end = matrix_parse_expr (s);
+ if (!loop->end)
+ goto error;
+
+ if (lex_match (s->lexer, T_BY))
+ {
+ loop->increment = matrix_parse_expr (s);
+ if (!loop->increment)
+ goto error;
+ }
+ }
+
+ if (lex_match_id (s->lexer, "IF"))
+ {
+ loop->top_condition = matrix_parse_expr (s);
+ if (!loop->top_condition)
+ goto error;
+ }
+
+ bool was_in_loop = s->in_loop;
+ s->in_loop = true;
+ bool ok = matrix_parse_commands (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_parse_expr (s);
+ if (!loop->bottom_condition)
+ goto error;
+ }
+
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_loop (struct loop_command *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_cmd_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
+static struct matrix_cmd *
+matrix_parse_break (struct matrix_state *s)
+{
+ if (!s->in_loop)
+ {
+ msg (SE, _("BREAK not inside LOOP."));
+ return NULL;
+ }
+
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
+ return cmd;
+}
+\f
+static struct matrix_cmd *
+matrix_parse_display (struct matrix_state *s)
+{
+ bool dictionary = false;
+ bool status = false;
+ while (lex_token (s->lexer) == T_ID)
+ {
+ if (lex_match_id (s->lexer, "DICTIONARY"))
+ dictionary = true;
+ else if (lex_match_id (s->lexer, "STATUS"))
+ status = true;
+ else
+ {
+ lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
+ return NULL;
+ }
+ }
+ if (!dictionary && !status)
+ dictionary = status = true;
+
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_DISPLAY,
+ .display = {
+ .state = s,
+ .dictionary = dictionary,
+ .status = status,
+ }
+ };
+ 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_cmd_execute_display (struct display_command *cmd)
+{
+ const struct matrix_state *s = cmd->state;
+
+ struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
+ pivot_dimension_create (
+ table, PIVOT_AXIS_COLUMN, N_("Property"),
+ N_("Rows"), N_("Columns"), 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]));
+ }
+ pivot_table_submit (table);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_release (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) { .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, _("Variable name expected."));
+ lex_get (s->lexer);
+
+ if (!lex_match (s->lexer, T_COMMA))
+ break;
+ }
+
+ return cmd;
+}
+
+static void
+matrix_cmd_execute_release (struct release_command *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
+static struct matrix_cmd *
+matrix_parse_save (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_SAVE,
+ .save = {
+ .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
+ }
+ };
+
+ struct save_command *save = &cmd->save;
+ save->expression = matrix_parse_exp (s);
+ if (!save->expression)
+ goto error;
+
+ while (lex_match (s->lexer, T_SLASH))
+ {
+ if (lex_match_id (s->lexer, "OUTFILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (save->outfile);
+ save->outfile = (lex_match (s->lexer, T_ASTERISK)
+ ? fh_inline_file ()
+ : fh_parse (s->lexer, FH_REF_FILE, s->session));
+ if (!save->outfile)
+ 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_destroy (save->variables);
+ if (!save->variables)
+ save->variables = xmalloc (sizeof *save->variables);
+ *save->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 (save->names);
+ save->names = matrix_parse_exp (s);
+ if (!save->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 (&save->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 (!save->outfile)
+ {
+ if (s->prev_save_outfile)
+ save->outfile = fh_ref (s->prev_save_outfile);
+ else
+ {
+ lex_sbc_missing ("OUTFILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_save_outfile);
+ s->prev_save_outfile = fh_ref (save->outfile);
+
+ if (save->variables && save->names)
+ {
+ msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
+ matrix_expr_destroy (save->names);
+ save->names = NULL;
+ }
+
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_save (const struct save_command *save)
+{
+ assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
+
+ gsl_matrix *m = matrix_expr_evaluate (save->expression);
+ if (!m)
+ return;
+
+ bool ok = true;
+ struct dictionary *dict = dict_create (get_default_encoding ());
+ struct string_array names = { .n = 0 };
+ if (save->variables)
+ string_array_clone (&names, save->variables);
+ else if (save->names)
+ {
+ gsl_matrix *nm = matrix_expr_evaluate (save->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));
+ if (dict_id_is_valid (dict, name, true))
+ string_array_append_nocopy (&names, name);
+ else
+ ok = false;
+ }
+ }
+ }
+
+ struct stringi_set strings;
+ stringi_set_clone (&strings, &save->strings);
+
+ for (size_t i = 0; dict_get_var_cnt (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 (SE, _("Duplicate variable name %s in SAVE statement."),
+ name);
+ ok = false;
+ }
+ }
+
+ if (!stringi_set_is_empty (&strings))
+ {
+ const char *example = stringi_set_node_get_string (
+ stringi_set_first (&strings));
+ msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
+ "with that name was found on SAVE.",
+ "STRINGS specified %2$zu variables, including %1$s, "
+ "whose names were not found on SAVE.",
+ stringi_set_count (&strings)),
+ example, stringi_set_count (&strings));
+ ok = false;
+ }
+ stringi_set_destroy (&strings);
+ string_array_destroy (&names);
+
+ if (!ok)
+ {
+ gsl_matrix_free (m);
+ dict_unref (dict);
+ return;
+ }
+
+ struct casewriter *writer = any_writer_open (save->outfile, dict);
+ if (!writer)
+ {
+ gsl_matrix_free (m);
+ dict_unref (dict);
+ return;
+ }
+
+ const struct caseproto *proto = dict_get_proto (dict);
+ 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);
+ }
+ casewriter_destroy (writer);
+
+ gsl_matrix_free (m);
+ dict_unref (dict);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_read (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_READ,
+ .read = { .format = FMT_F },
+ };
+
+ struct file_handle *fh = NULL;
+ char *encoding = NULL;
+ struct read_command *read = &cmd->read;
+ read->dst = matrix_lvalue_parse (s);
+ if (!read->dst)
+ goto error;
+
+ 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);
+
+ 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;
+ 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);
+ lex_get (s->lexer);
+
+ if (record_width % by)
+ {
+ msg (SE, _("BY %d does not evenly divide record width %d."),
+ by, record_width);
+ goto error;
+ }
+ }
+ else
+ by = 0;
+ }
+ else if (lex_match_id (s->lexer, "SIZE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ 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 ("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;
+
+ 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))
+ {
+ struct fmt_spec format;
+ if (!parse_format_specifier (s->lexer, &format))
+ goto error;
+ read->format = format.type;
+ read->w = format.w;
+ read->d = format.d;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
+ "REREAD", "FORMAT");
+ goto error;
+ }
+ }
+
+ if (!read->c1)
+ {
+ lex_sbc_missing ("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)."));
+ goto error;
+ }
+
+ if (!fh)
+ {
+ if (s->prev_read_file)
+ fh = fh_ref (s->prev_read_file);
+ else
+ {
+ lex_sbc_missing ("FILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_read_file);
+ s->prev_read_file = fh_ref (fh);
+
+ read->rf = read_file_create (s, fh);
+ 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);
+ goto error;
+ }
+ int w = (repetitions ? record_width / repetitions
+ : read->w ? read->w
+ : by);
+ if (by && w != by)
+ {
+ if (repetitions)
+ msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
+ "which implies field width %d, "
+ "but BY specifies field width %d."),
+ repetitions, record_width, w, by);
+ else
+ msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
+ w, by);
+ goto error;
+ }
+ read->w = w;
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ free (encoding);
+ return NULL;
+}
+
+static void
+matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
+ gsl_matrix *m, struct substring p, size_t y, size_t x)
+{
+ const char *input_encoding = dfm_reader_get_encoding (reader);
+ union value v;
+ char *error = data_in (p, input_encoding, read->format,
+ settings_get_fmt_settings (), &v, 0, NULL);
+ /* XXX report error if value is missing */
+ if (error)
+ msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
+ else
+ {
+ gsl_matrix_set (m, y, x, v.f);
+ if (read->symmetric && x != y)
+ gsl_matrix_set (m, x, y, v.f);
+ }
+}
+
+static bool
+matrix_read_line (struct read_command *read, struct dfm_reader *reader,
+ struct substring *line)
+{
+ if (dfm_eof (reader))
+ {
+ msg (SE, _("Unexpected end of file reading matrix data."));
+ return false;
+ }
+ dfm_expand_tabs (reader);
+ *line = ss_substr (dfm_get_record (reader),
+ read->c1 - 1, read->c2 - read->c1);
+ return true;
+}
+
+static void
+matrix_read (struct read_command *read, struct dfm_reader *reader,
+ gsl_matrix *m)
+{
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ size_t nx = read->symmetric ? y + 1 : m->size2;
+
+ struct substring line = ss_empty ();
+ 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 (read, reader, &line))
+ return;
+ dfm_forward_record (reader);
+ }
+
+ ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
+ }
+ else
+ {
+ if (!matrix_read_line (read, reader, &line))
+ 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);
+ }
+
+ if (read->w)
+ dfm_forward_record (reader);
+ else
+ {
+ ss_ltrim (&line, ss_cstr (" ,"));
+ if (!ss_is_empty (line))
+ msg (SW, _("Trailing garbage on line \"%.*s\""),
+ (int) line.length, line.string);
+ }
+ }
+}
+
+static void
+matrix_cmd_execute_read (struct read_command *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 (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
+ gsl_matrix_free (m);
+ free (iv0.indexes);
+ free (iv1.indexes);
+ 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 (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
+ gsl_matrix_free (m);
+ free (iv0.indexes);
+ free (iv1.indexes);
+ return;
+ }
+ gsl_matrix_free (m);
+
+ if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
+ {
+ msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
+ free (iv0.indexes);
+ free (iv1.indexes);
+ 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 (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
+ "(%zu,%zu)."),
+ size[0], size[1],
+ submatrix_size[0], submatrix_size[1]);
+ free (iv0.indexes);
+ free (iv1.indexes);
+ 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 (SE, _("Cannot read matrix with non-square dimensions (%zu,%zu) "
+ "using READ with MODE=SYMMETRIC."),
+ size[0], size[1]);
+ free (iv0.indexes);
+ free (iv1.indexes);
+ return;
+ }
+ gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
+ matrix_read (read, reader, tmp);
+ matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_write (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_WRITE,
+ .write = { .format = FMT_F },
+ };
+
+ struct write_command *write = &cmd->write;
+ write->expression = matrix_parse_exp (s);
+ if (!write->expression)
+ goto error;
+
+ 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, "OUTFILE"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (write->outfile);
+ write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
+ if (!write->outfile)
+ 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 (write->encoding);
+ write->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);
+
+ 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;
+ 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 = lex_integer (s->lexer);
+ lex_get (s->lexer);
+
+ if (record_width % by)
+ {
+ msg (SE, _("BY %d does not evenly divide record width %d."),
+ by, record_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 (seen_format)
+ {
+ lex_sbc_only_once ("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;
+
+ const char *p = lex_tokcstr (s->lexer);
+ if (c_isdigit (p[0]))
+ {
+ repetitions = atoi (p);
+ p += strspn (p, "0123456789");
+ if (!fmt_from_name (p, &write->format))
+ {
+ lex_error (s->lexer, _("Unknown format %s."), p);
+ goto error;
+ }
+ lex_get (s->lexer);
+ }
+ else if (!fmt_from_name (p, &write->format))
+ {
+ struct fmt_spec format;
+ if (!parse_format_specifier (s->lexer, &format))
+ goto error;
+ write->format = format.type;
+ write->w = format.w;
+ write->d = format.d;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
+ "HOLD", "FORMAT");
+ goto error;
+ }
+ }
+
+ if (!write->c1)
+ {
+ lex_sbc_missing ("FIELD");
+ goto error;
+ }
+
+ if (!write->outfile)
+ {
+ if (s->prev_write_outfile)
+ write->outfile = fh_ref (s->prev_write_outfile);
+ else
+ {
+ lex_sbc_missing ("OUTFILE");
+ goto error;
+ }
+ }
+ fh_unref (s->prev_write_outfile);
+ s->prev_write_outfile = fh_ref (write->outfile);
+
+ /* 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);
+ goto error;
+ }
+ int w = (repetitions ? record_width / repetitions
+ : write->w ? write->w
+ : by);
+ if (by && w != by)
+ {
+ if (repetitions)
+ msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
+ "which implies field width %d, "
+ "but BY specifies field width %d."),
+ repetitions, record_width, w, by);
+ else
+ msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
+ w, by);
+ goto error;
+ }
+ write->w = w;
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_write (struct write_command *write)
+{
+ gsl_matrix *m = matrix_expr_evaluate (write->expression);
+ if (!m)
+ return;
+
+ if (write->triangular && m->size1 != m->size2)
+ {
+ msg (SE, _("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 = dfm_open_writer (write->outfile, write->encoding);
+ if (!writer)
+ {
+ gsl_matrix_free (m);
+ return;
+ }
+
+ const struct fmt_settings *settings = settings_get_fmt_settings ();
+ struct fmt_spec format = {
+ .type = write->format,
+ .w = write->w ? write->w : 40,
+ .d = write->d
+ };
+ struct u8_line line = U8_LINE_INITIALIZER;
+ for (size_t y = 0; y < m->size1; y++)
+ {
+ size_t nx = write->triangular ? y + 1 : m->size2;
+ int x0 = write->c1;
+ for (size_t x = 0; x < nx; x++)
+ {
+ /* XXX string values */
+ union value v = { .f = gsl_matrix_get (m, y, x) };
+ char *s = (write->w
+ ? data_out (&v, NULL, &format, settings)
+ : data_out_stretchy (&v, NULL, &format, settings, NULL));
+ 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->w ? write->w : width + 1;
+ }
+
+ dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
+ u8_line_clear (&line);
+ }
+ u8_line_destroy (&line);
+ dfm_close_writer (writer);
+
+ gsl_matrix_free (m);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_get (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_GET,
+ .get = {
+ .user = { .treatment = MGET_ERROR },
+ .system = { .treatment = MGET_ERROR },
+ }
+ };
+
+ struct get_command *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"))
+ {
+ if (get->variables.n)
+ {
+ lex_error (s->lexer, _("FILE must precede VARIABLES"));
+ goto error;
+ }
+ lex_match (s->lexer, T_EQUALS);
+
+ fh_unref (get->file);
+ get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
+ if (!get->file)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "ENCODING"))
+ {
+ if (get->variables.n)
+ {
+ lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
+ goto error;
+ }
+ 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);
+
+ struct dictionary *dict = NULL;
+ if (!get->file)
+ {
+ dict = dataset_dict (s->dataset);
+ if (dict_get_var_cnt (dict) == 0)
+ {
+ lex_error (s->lexer, _("GET cannot read empty active file."));
+ goto error;
+ }
+ }
+ else
+ {
+ struct casereader *reader = any_reader_open_and_decode (
+ get->file, get->encoding, &dict, NULL);
+ if (!reader)
+ goto error;
+ casereader_destroy (reader);
+ }
+
+ struct variable **vars;
+ size_t n_vars;
+ bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
+ PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
+ if (!ok)
+ {
+ dict_unref (dict);
+ goto error;
+ }
+
+ string_array_clear (&get->variables);
+ for (size_t i = 0; i < n_vars; i++)
+ string_array_append (&get->variables, var_get_name (vars[i]));
+ free (vars);
+ dict_unref (dict);
+ }
+ 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, NULL);
+ goto error;
+ }
+ }
+ else if (lex_match_id (s->lexer, "SYSMIS"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+ 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, NULL);
+ goto error;
+ }
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
+ "MISSING", "SYSMIS");
+ goto error;
+ }
+ }
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_get (struct get_command *get)
+{
+ assert (get->file); /* XXX */
+
+ struct dictionary *dict;
+ struct casereader *reader = any_reader_open_and_decode (
+ get->file, get->encoding, &dict, NULL);
+ if (!reader)
+ return;
+
+ const struct variable **vars = xnmalloc (
+ get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
+ sizeof *vars);
+ size_t n_vars = 0;
+
+ if (get->variables.n)
+ {
+ for (size_t i = 0; i < get->variables.n; i++)
+ {
+ const char *name = get->variables.strings[i];
+ const struct variable *var = dict_lookup_var (dict, name);
+ if (!var)
+ {
+ msg (SE, _("GET: Data file does not contain variable %s."),
+ name);
+ dict_unref (dict);
+ free (vars);
+ return;
+ }
+ if (!var_is_numeric (var))
+ {
+ msg (SE, _("GET: Variable %s is not numeric."), name);
+ dict_unref (dict);
+ free (vars);
+ return;
+ }
+ vars[n_vars++] = var;
+ }
+ }
+ else
+ {
+ for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
+ {
+ const struct variable *var = dict_get_var (dict, i);
+ if (!var_is_numeric (var))
+ {
+ msg (SE, _("GET: Variable %s is not numeric."),
+ var_get_name (var));
+ dict_unref (dict);
+ free (vars);
+ return;
+ }
+ vars[n_vars++] = var;
+ }
+ }
+
+ 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 (SE, _("GET: 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 (SE, _("GET: 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++;
+ }
+ casereader_destroy (reader);
+ if (!error)
+ {
+ m->size1 = n_rows;
+ matrix_lvalue_evaluate_and_assign (get->dst, m);
+ }
+ else
+ gsl_matrix_free (m);
+ dict_unref (dict);
+ free (vars);
+}
+\f
+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)
+{
+ lex_match (lexer, T_EQUALS);
+
+ string_array_clear (sa);
+
+ struct dictionary *dict = dict_create (get_default_encoding ());
+ dict_create_var_assert (dict, "ROWTYPE_", 8);
+ dict_create_var_assert (dict, "VARNAME_", 8);
+ char **names;
+ size_t n_names;
+ bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
+ PV_NO_DUPLICATE | PV_NO_SCRATCH);
+ dict_unref (dict);
+
+ if (ok)
+ {
+ string_array_clear (sa);
+ sa->strings = names;
+ sa->n = sa->allocated = n_names;
+ }
+ return ok;
+}
+
+static void
+msave_common_uninit (struct msave_common *common)
+{
+ if (common)
+ {
+ fh_unref (common->outfile);
+ string_array_destroy (&common->variables);
+ string_array_destroy (&common->fnames);
+ string_array_destroy (&common->snames);
+ }
+}
+
+static bool
+compare_variables (const char *keyword,
+ const struct string_array *new,
+ const struct string_array *old)
+{
+ if (new->n)
+ {
+ if (!old->n)
+ {
+ msg (SE, _("%s may only be specified on MSAVE if it was specified "
+ "on the first MSAVE within MATRIX."), keyword);
+ return false;
+ }
+ else if (!string_array_equal_case (old, new))
+ {
+ msg (SE, _("%s must specify the same variables each time within "
+ "a given MATRIX."), keyword);
+ return false;
+ }
+ }
+ return true;
+}
+
+static struct matrix_cmd *
+matrix_parse_msave (struct matrix_state *s)
+{
+ struct msave_common common = { .outfile = NULL };
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
+
+ struct msave_command *msave = &cmd->msave;
+ if (lex_token (s->lexer) == T_ID)
+ msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
+ msave->expr = matrix_parse_exp (s);
+ if (!msave->expr)
+ return NULL;
+
+ 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);
+ common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
+ if (!common.outfile)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "VARIABLES"))
+ {
+ if (!parse_var_names (s->lexer, &common.variables))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "FNAMES"))
+ {
+ if (!parse_var_names (s->lexer, &common.fnames))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "SNAMES"))
+ {
+ if (!parse_var_names (s->lexer, &common.snames))
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "SPLIT"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ matrix_expr_destroy (msave->splits);
+ msave->splits = matrix_parse_exp (s);
+ if (!msave->splits)
+ goto error;
+ }
+ else if (lex_match_id (s->lexer, "FACTOR"))
+ {
+ lex_match (s->lexer, T_EQUALS);
+
+ matrix_expr_destroy (msave->factors);
+ msave->factors = matrix_parse_exp (s);
+ if (!msave->factors)
+ goto error;
+ }
+ else
+ {
+ lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
+ "FNAMES", "SNAMES", "SPLIT", "FACTOR");
+ goto error;
+ }
+ }
+ if (!msave->rowtype)
+ {
+ lex_sbc_missing ("TYPE");
+ goto error;
+ }
+ common.has_splits = msave->splits || common.snames.n;
+ common.has_factors = msave->factors || common.fnames.n;
+
+ struct msave_common *c = s->common ? s->common : &common;
+ if (c->fnames.n && !msave->factors)
+ {
+ msg (SE, _("FNAMES requires FACTOR."));
+ goto error;
+ }
+ if (c->snames.n && !msave->splits)
+ {
+ msg (SE, _("SNAMES requires SPLIT."));
+ goto error;
+ }
+ if (c->has_factors && !common.has_factors)
+ {
+ msg (SE, _("%s is required because it was present on the first "
+ "MSAVE in this MATRIX command."), "FACTOR");
+ goto error;
+ }
+ if (c->has_splits && !common.has_splits)
+ {
+ msg (SE, _("%s is required because it was present on the first "
+ "MSAVE in this MATRIX command."), "SPLIT");
+ goto error;
+ }
+
+ if (!s->common)
+ {
+ if (!common.outfile)
+ {
+ lex_sbc_missing ("OUTFILE");
+ goto error;
+ }
+ s->common = xmemdup (&common, sizeof common);
+ }
+ else
+ {
+ if (common.outfile)
+ {
+ bool same = common.outfile == s->common->outfile;
+ fh_unref (common.outfile);
+ if (!same)
+ {
+ msg (SE, _("OUTFILE must name the same file on each MSAVE "
+ "within a single MATRIX command."));
+ goto error;
+ }
+ }
+ if (!compare_variables ("VARIABLES",
+ &common.variables, &s->common->variables)
+ || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
+ || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
+ goto error;
+ msave_common_uninit (&common);
+ }
+ msave->common = s->common;
+ if (!msave->varname_)
+ msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
+ return cmd;
+
+error:
+ msave_common_uninit (&common);
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static gsl_vector *
+matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
+{
+ gsl_matrix *m = matrix_expr_evaluate (e);
+ if (!m)
+ return NULL;
+
+ if (!is_vector (m))
+ {
+ msg (SE, _("%s expression must evaluate to vector, not a matrix with "
+ "dimensions (%zu,%zu)."),
+ 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)
+ {
+ msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
+ goto error;
+ }
+
+ dict_create_var_assert (dict, "ROWTYPE_", 8);
+
+ const char *dup_factor = msave_add_vars (dict, &common->fnames);
+ if (dup_factor)
+ {
+ msg (SE, _("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 (SE, _("Duplicate or invalid variable name %s."), dup_var);
+ goto error;
+ }
+
+ return dict;
+
+error:
+ dict_unref (dict);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_msave (struct msave_command *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));
+
+ if (m->size2 != common->variables.n)
+ {
+ msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
+ 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));
+
+ if (factors->size != common->fnames.n)
+ {
+ msg (SE, _("There are %zu factor variables, "
+ "but %zu split 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->fnames.n)
+ for (size_t i = 0; i < splits->size; i++)
+ string_array_append_nocopy (&common->fnames,
+ xasprintf ("SPL%zu", i + 1));
+
+ if (splits->size != common->fnames.n)
+ {
+ msg (SE, _("There are %zu split variables, "
+ "but %zu split values were supplied."),
+ common->fnames.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;
+ }
+
+ 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_. */
+ buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
+ msave->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
+static struct matrix_cmd *
+matrix_parse_mget (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
+
+ struct mget_command *mget = &cmd->mget;
+
+ for (;;)
+ {
+ 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, "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;
+ }
+ if (lex_token (s->lexer) == T_ENDCMD)
+ break;
+
+ if (!lex_force_match (s->lexer, T_SLASH))
+ goto error;
+ }
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static const struct variable *
+get_a8_var (const struct dictionary *d, const char *name)
+{
+ const struct variable *v = dict_lookup_var (d, name);
+ if (!v)
+ {
+ msg (SE, _("Matrix data file lacks %s variable."), name);
+ return NULL;
+ }
+ if (var_get_width (v) != 8)
+ {
+ msg (SE, _("%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
+is_rowtype (const union value *v, const char *rowtype)
+{
+ struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
+ ss_rtrim (&vs, ss_cstr (" "));
+ return ss_equals_case (vs, ss_cstr (rowtype));
+}
+
+static void
+matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
+ const struct dictionary *d,
+ const struct variable *rowtype_var,
+ struct matrix_state *s, size_t si, size_t fi,
+ size_t cs, size_t cn)
+{
+ if (!n_rows)
+ return;
+
+ const union value *rowtype_ = case_data (rows[0], rowtype_var);
+ const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
+ : is_rowtype (rowtype_, "CORR") ? "CR"
+ : is_rowtype (rowtype_, "MEAN") ? "MN"
+ : is_rowtype (rowtype_, "STDDEV") ? "SD"
+ : is_rowtype (rowtype_, "N") ? "NC"
+ : "CN");
+
+ struct string name = DS_EMPTY_INITIALIZER;
+ ds_put_cstr (&name, name_prefix);
+ if (fi > 0)
+ 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;
+ }
+
+ 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, MV_ANY))
+ {
+ n_missing++;
+ value = 0.0;
+ }
+ gsl_matrix_set (m, y, x, value);
+ }
+ }
+
+ 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:
+ ds_destroy (&name);
+ for (size_t y = 0; y < n_rows; y++)
+ case_unref (rows[y]);
+}
+
+static bool
+var_changed (const struct ccase *ca, const struct ccase *cb,
+ const struct variable *var)
+{
+ return !value_equal (case_data (ca, var), case_data (cb, var),
+ var_get_width (var));
+}
+
+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 void
+matrix_cmd_execute_mget (struct mget_command *mget)
+{
+ struct dictionary *d;
+ struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
+ &d, NULL);
+ if (!r)
+ return;
+
+ const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
+ const struct variable *varname_ = get_a8_var (d, "VARNAME_");
+ if (!rowtype_ || !varname_)
+ goto exit;
+
+ if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
+ {
+ msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
+ goto exit;
+ }
+ if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
+ {
+ msg (SE, _("Matrix data file contains no continuous variables."));
+ goto exit;
+ }
+
+ for (size_t i = 0; i < dict_get_var_cnt (d); i++)
+ {
+ const struct variable *v = dict_get_var (d, i);
+ if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
+ {
+ msg (SE,
+ _("Matrix data file contains unexpected string variable %s."),
+ var_get_name (v));
+ goto exit;
+ }
+ }
+
+ /* 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_var_cnt (d) - cs;
+ struct ccase *cc = NULL;
+
+ /* 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 sd = vars_changed (sc, c, d, ss, sn);
+ bool fd = sd || vars_changed (fc, c, d, fs, fn);
+ bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
+ if (sd)
+ {
+ si++;
+ case_unref (sc);
+ sc = case_ref (c);
+ }
+ if (fd)
+ {
+ fi++;
+ case_unref (fc);
+ fc = case_ref (c);
+ }
+ if (md)
+ {
+ matrix_mget_commit_var (rows, n_rows, d, rowtype_,
+ mget->state, si, fi, cs, cn);
+ 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;
+ }
+ matrix_mget_commit_var (rows, n_rows, d, rowtype_,
+ mget->state, si, fi, cs, cn);
+ free (rows);
+
+exit:
+ casereader_destroy (r);
+}
+\f
+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_cmd *
+matrix_parse_eigen (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_EIGEN,
+ .eigen = { .expr = NULL }
+ };
+
+ struct eigen_command *eigen = &cmd->eigen;
+ if (!lex_force_match (s->lexer, T_LPAREN))
+ goto error;
+ eigen->expr = matrix_parse_expr (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_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_eigen (struct eigen_command *eigen)
+{
+ gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
+ if (!A)
+ return;
+ if (!is_symmetric (A))
+ {
+ msg (SE, _("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
+static struct matrix_cmd *
+matrix_parse_setdiag (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_SETDIAG,
+ .setdiag = { .dst = NULL }
+ };
+
+ struct setdiag_command *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_parse_expr (s);
+ if (!setdiag->expr)
+ goto error;
+
+ if (!lex_force_match (s->lexer, T_RPAREN))
+ goto error;
+
+ return cmd;
+
+error:
+ matrix_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
+{
+ gsl_matrix *dst = setdiag->dst->value;
+ if (!dst)
+ {
+ msg (SE, _("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 (SE, _("SETDIAG argument 2 must be a scalar or a vector but it "
+ "has dimensions (%zu,%zu)."),
+ src->size1, src->size2);
+ gsl_matrix_free (src);
+}
+\f
+static struct matrix_cmd *
+matrix_parse_svd (struct matrix_state *s)
+{
+ struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
+ *cmd = (struct matrix_cmd) {
+ .type = MCMD_SVD,
+ .svd = { .expr = NULL }
+ };
+
+ struct svd_command *svd = &cmd->svd;
+ if (!lex_force_match (s->lexer, T_LPAREN))
+ goto error;
+ svd->expr = matrix_parse_expr (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_cmd_destroy (cmd);
+ return NULL;
+}
+
+static void
+matrix_cmd_execute_svd (struct svd_command *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 *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
+static bool
+matrix_cmd_execute (struct matrix_cmd *cmd)
+{
+ switch (cmd->type)
+ {
+ case MCMD_COMPUTE:
+ matrix_cmd_execute_compute (&cmd->compute);
+ break;
+
+ case MCMD_PRINT:
+ matrix_cmd_execute_print (&cmd->print);
+ break;
+
+ case MCMD_DO_IF:
+ return matrix_cmd_execute_do_if (&cmd->do_if);
+
+ case MCMD_LOOP:
+ matrix_cmd_execute_loop (&cmd->loop);
+ break;
+
+ case MCMD_BREAK:
+ return false;
+
+ case MCMD_DISPLAY:
+ matrix_cmd_execute_display (&cmd->display);
+ break;
+
+ case MCMD_RELEASE:
+ matrix_cmd_execute_release (&cmd->release);
+ break;
+
+ case MCMD_SAVE:
+ matrix_cmd_execute_save (&cmd->save);
+ break;
+
+ case MCMD_READ:
+ matrix_cmd_execute_read (&cmd->read);
+ break;
+
+ case MCMD_WRITE:
+ matrix_cmd_execute_write (&cmd->write);
+ break;
+
+ case MCMD_GET:
+ matrix_cmd_execute_get (&cmd->get);
+ break;
+
+ case MCMD_MSAVE:
+ matrix_cmd_execute_msave (&cmd->msave);
+ break;
+
+ case MCMD_MGET:
+ matrix_cmd_execute_mget (&cmd->mget);
+ break;
+
+ case MCMD_EIGEN:
+ matrix_cmd_execute_eigen (&cmd->eigen);
+ break;
+
+ case MCMD_SETDIAG:
+ matrix_cmd_execute_setdiag (&cmd->setdiag);
+ break;
+
+ case MCMD_SVD:
+ matrix_cmd_execute_svd (&cmd->svd);
+ break;
+ }
+
+ return true;
+}
+
+struct matrix_command_name
+ {
+ const char *name;
+ struct matrix_cmd *(*parse) (struct matrix_state *);
+ };
+
+static const struct matrix_command_name *
+matrix_parse_command_name (struct lexer *lexer)
+{
+ static const struct matrix_command_name commands[] = {
+ { "COMPUTE", matrix_parse_compute },
+ { "CALL EIGEN", matrix_parse_eigen },
+ { "CALL SETDIAG", matrix_parse_setdiag },
+ { "CALL SVD", matrix_parse_svd },
+ { "PRINT", matrix_parse_print },
+ { "DO IF", matrix_parse_do_if },
+ { "LOOP", matrix_parse_loop },
+ { "BREAK", matrix_parse_break },
+ { "READ", matrix_parse_read },
+ { "WRITE", matrix_parse_write },
+ { "GET", matrix_parse_get },
+ { "SAVE", matrix_parse_save },
+ { "MGET", matrix_parse_mget },
+ { "MSAVE", matrix_parse_msave },
+ { "DISPLAY", matrix_parse_display },
+ { "RELEASE", matrix_parse_release },
+ };
+ 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_cmd *
+matrix_parse_command (struct matrix_state *s)
+{
+ size_t nesting_level = SIZE_MAX;
+
+ struct matrix_cmd *c = NULL;
+ const struct matrix_command_name *cmd = matrix_parse_command_name (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)
+ 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 = {
+ .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_cmd *c = matrix_parse_command (&state);
+ if (c)
+ {
+ matrix_cmd_execute (c);
+ matrix_cmd_destroy (c);
+ }
+ }
+
+ if (state.common)
+ {
+ dict_unref (state.common->dict);
+ casewriter_destroy (state.common->writer);
+ free (state.common);
+ }
+
+ return CMD_SUCCESS;
+}
projects / pspp / commitdiff