1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_eigen.h>
21 #include <gsl/gsl_linalg.h>
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_permutation.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_vector.h>
30 #include "data/any-reader.h"
31 #include "data/any-writer.h"
32 #include "data/casereader.h"
33 #include "data/casewriter.h"
34 #include "data/data-in.h"
35 #include "data/data-out.h"
36 #include "data/dataset.h"
37 #include "data/dictionary.h"
38 #include "data/file-handle-def.h"
39 #include "language/command.h"
40 #include "language/data-io/data-reader.h"
41 #include "language/data-io/data-writer.h"
42 #include "language/data-io/file-handle.h"
43 #include "language/lexer/format-parser.h"
44 #include "language/lexer/lexer.h"
45 #include "language/lexer/variable-parser.h"
46 #include "libpspp/array.h"
47 #include "libpspp/assertion.h"
48 #include "libpspp/compiler.h"
49 #include "libpspp/hmap.h"
50 #include "libpspp/i18n.h"
51 #include "libpspp/misc.h"
52 #include "libpspp/str.h"
53 #include "libpspp/string-array.h"
54 #include "libpspp/stringi-set.h"
55 #include "libpspp/u8-line.h"
56 #include "math/random.h"
57 #include "output/driver.h"
58 #include "output/output-item.h"
59 #include "output/pivot-table.h"
61 #include "gl/c-ctype.h"
62 #include "gl/minmax.h"
66 #define _(msgid) gettext (msgid)
67 #define N_(msgid) (msgid)
71 struct hmap_node hmap_node;
79 struct file_handle *outfile;
80 struct string_array variables;
81 struct string_array fnames;
82 struct string_array snames;
87 /* Execution state. */
88 struct dictionary *dict;
89 struct casewriter *writer;
94 struct file_handle *file;
95 struct dfm_reader *reader;
101 struct file_handle *file;
102 struct dfm_writer *writer;
104 struct u8_line *held;
109 struct dataset *dataset;
110 struct session *session;
114 struct file_handle *prev_save_outfile;
115 struct msave_common *common;
117 struct file_handle *prev_read_file;
118 struct read_file **read_files;
121 struct file_handle *prev_write_file;
122 struct write_file **write_files;
123 size_t n_write_files;
126 static struct matrix_var *
127 matrix_var_lookup (struct matrix_state *s, struct substring name)
129 struct matrix_var *var;
131 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
132 utf8_hash_case_substring (name, 0), &s->vars)
133 if (!utf8_sscasecmp (ss_cstr (var->name), name))
139 static struct matrix_var *
140 matrix_var_create (struct matrix_state *s, struct substring name)
142 struct matrix_var *var = xmalloc (sizeof *var);
143 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
144 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
149 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
151 gsl_matrix_free (var->value);
155 #define MATRIX_FUNCTIONS \
211 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
212 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
213 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
214 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
215 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
216 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
217 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
218 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
219 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
220 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
221 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
223 typedef gsl_matrix *matrix_proto_m_d (double);
224 typedef gsl_matrix *matrix_proto_m_dd (double, double);
225 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
226 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
227 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
228 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
229 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
230 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
231 typedef double matrix_proto_d_m (gsl_matrix *);
232 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
233 typedef gsl_matrix *matrix_proto_IDENT (double, double);
235 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
239 /* Matrix expressions. */
246 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
250 /* Elementwise and scalar arithmetic. */
251 MOP_NEGATE, /* unary - */
252 MOP_ADD_ELEMS, /* + */
253 MOP_SUB_ELEMS, /* - */
254 MOP_MUL_ELEMS, /* &* */
255 MOP_DIV_ELEMS, /* / and &/ */
256 MOP_EXP_ELEMS, /* &** */
258 MOP_SEQ_BY, /* a:b:c */
260 /* Matrix arithmetic. */
262 MOP_EXP_MAT, /* ** */
279 MOP_PASTE_HORZ, /* a, b, c, ... */
280 MOP_PASTE_VERT, /* a; b; c; ... */
284 MOP_VEC_INDEX, /* x(y) */
285 MOP_VEC_ALL, /* x(:) */
286 MOP_MAT_INDEX, /* x(y,z) */
287 MOP_ROW_INDEX, /* x(y,:) */
288 MOP_COL_INDEX, /* x(:,z) */
295 MOP_EOF, /* EOF('file') */
303 struct matrix_expr **subs;
308 struct matrix_var *variable;
309 struct read_file *eof;
314 matrix_expr_destroy (struct matrix_expr *e)
321 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
352 for (size_t i = 0; i < e->n_subs; i++)
353 matrix_expr_destroy (e->subs[i]);
365 static struct matrix_expr *
366 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
369 struct matrix_expr *e = xmalloc (sizeof *e);
370 *e = (struct matrix_expr) {
372 .subs = xmemdup (subs, n_subs * sizeof *subs),
378 static struct matrix_expr *
379 matrix_expr_create_0 (enum matrix_op op)
381 struct matrix_expr *sub;
382 return matrix_expr_create_subs (op, &sub, 0);
385 static struct matrix_expr *
386 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
388 return matrix_expr_create_subs (op, &sub, 1);
391 static struct matrix_expr *
392 matrix_expr_create_2 (enum matrix_op op,
393 struct matrix_expr *sub0, struct matrix_expr *sub1)
395 struct matrix_expr *subs[] = { sub0, sub1 };
396 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
399 static struct matrix_expr *
400 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
401 struct matrix_expr *sub1, struct matrix_expr *sub2)
403 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
404 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
407 static struct matrix_expr *
408 matrix_expr_create_number (double number)
410 struct matrix_expr *e = xmalloc (sizeof *e);
411 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
415 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
417 static struct matrix_expr *
418 matrix_parse_curly_comma (struct matrix_state *s)
420 struct matrix_expr *lhs = matrix_parse_or_xor (s);
424 while (lex_match (s->lexer, T_COMMA))
426 struct matrix_expr *rhs = matrix_parse_or_xor (s);
429 matrix_expr_destroy (lhs);
432 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
437 static struct matrix_expr *
438 matrix_parse_curly_semi (struct matrix_state *s)
440 if (lex_token (s->lexer) == T_RCURLY)
441 return matrix_expr_create_0 (MOP_EMPTY);
443 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
447 while (lex_match (s->lexer, T_SEMICOLON))
449 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
452 matrix_expr_destroy (lhs);
455 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
460 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
461 for (size_t Y = 0; Y < (M)->size1; Y++) \
462 for (size_t X = 0; X < (M)->size2; X++) \
463 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
466 is_vector (const gsl_matrix *m)
468 return m->size1 <= 1 || m->size2 <= 1;
472 to_vector (gsl_matrix *m)
474 return (m->size1 == 1
475 ? gsl_matrix_row (m, 0).vector
476 : gsl_matrix_column (m, 0).vector);
481 matrix_eval_ABS (gsl_matrix *m)
483 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
489 matrix_eval_ALL (gsl_matrix *m)
491 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
498 matrix_eval_ANY (gsl_matrix *m)
500 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
507 matrix_eval_ARSIN (gsl_matrix *m)
509 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
515 matrix_eval_ARTAN (gsl_matrix *m)
517 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
523 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
527 for (size_t i = 0; i < n; i++)
532 gsl_matrix *block = gsl_matrix_calloc (r, c);
534 for (size_t i = 0; i < n; i++)
536 for (size_t y = 0; y < m[i]->size1; y++)
537 for (size_t x = 0; x < m[i]->size2; x++)
538 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
546 matrix_eval_CHOL (gsl_matrix *m)
548 if (!gsl_linalg_cholesky_decomp1 (m))
550 for (size_t y = 0; y < m->size1; y++)
551 for (size_t x = y + 1; x < m->size2; x++)
552 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
554 for (size_t y = 0; y < m->size1; y++)
555 for (size_t x = 0; x < y; x++)
556 gsl_matrix_set (m, y, x, 0);
561 msg (SE, _("Input to CHOL function is not positive-definite."));
567 matrix_eval_col_extremum (gsl_matrix *m, bool min)
572 return gsl_matrix_alloc (1, 0);
574 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
575 for (size_t x = 0; x < m->size2; x++)
577 double ext = gsl_matrix_get (m, 0, x);
578 for (size_t y = 1; y < m->size1; y++)
580 double value = gsl_matrix_get (m, y, x);
581 if (min ? value < ext : value > ext)
584 gsl_matrix_set (cext, 0, x, ext);
590 matrix_eval_CMAX (gsl_matrix *m)
592 return matrix_eval_col_extremum (m, false);
596 matrix_eval_CMIN (gsl_matrix *m)
598 return matrix_eval_col_extremum (m, true);
602 matrix_eval_COS (gsl_matrix *m)
604 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
610 matrix_eval_col_sum (gsl_matrix *m, bool square)
615 return gsl_matrix_alloc (1, 0);
617 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
618 for (size_t x = 0; x < m->size2; x++)
621 for (size_t y = 0; y < m->size1; y++)
623 double d = gsl_matrix_get (m, y, x);
624 sum += square ? pow2 (d) : d;
626 gsl_matrix_set (result, 0, x, sum);
632 matrix_eval_CSSQ (gsl_matrix *m)
634 return matrix_eval_col_sum (m, true);
638 matrix_eval_CSUM (gsl_matrix *m)
640 return matrix_eval_col_sum (m, false);
644 compare_double_3way (const void *a_, const void *b_)
646 const double *a = a_;
647 const double *b = b_;
648 return *a < *b ? -1 : *a > *b;
652 matrix_eval_DESIGN (gsl_matrix *m)
654 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
655 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
656 gsl_matrix_transpose_memcpy (&m2, m);
658 for (size_t y = 0; y < m2.size1; y++)
659 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
661 size_t *n = xcalloc (m2.size1, sizeof *n);
663 for (size_t i = 0; i < m2.size1; i++)
665 double *row = tmp + m2.size2 * i;
666 for (size_t j = 0; j < m2.size2; )
669 for (k = j + 1; k < m2.size2; k++)
670 if (row[j] != row[k])
672 row[n[i]++] = row[j];
677 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
682 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
684 for (size_t i = 0; i < m->size2; i++)
689 const double *unique = tmp + m2.size2 * i;
690 for (size_t j = 0; j < n[i]; j++, x++)
692 double value = unique[j];
693 for (size_t y = 0; y < m->size1; y++)
694 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
705 matrix_eval_DET (gsl_matrix *m)
707 gsl_permutation *p = gsl_permutation_alloc (m->size1);
709 gsl_linalg_LU_decomp (m, p, &signum);
710 gsl_permutation_free (p);
711 return gsl_linalg_LU_det (m, signum);
715 matrix_eval_DIAG (gsl_matrix *m)
717 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
718 for (size_t i = 0; i < diag->size1; i++)
719 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
724 is_symmetric (const gsl_matrix *m)
726 if (m->size1 != m->size2)
729 for (size_t y = 0; y < m->size1; y++)
730 for (size_t x = 0; x < y; x++)
731 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
738 compare_double_desc (const void *a_, const void *b_)
740 const double *a = a_;
741 const double *b = b_;
742 return *a > *b ? -1 : *a < *b;
746 matrix_eval_EVAL (gsl_matrix *m)
748 if (!is_symmetric (m))
750 msg (SE, _("Argument of EVAL must be symmetric."));
754 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
755 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
756 gsl_vector v_eval = to_vector (eval);
757 gsl_eigen_symm (m, &v_eval, w);
758 gsl_eigen_symm_free (w);
760 assert (v_eval.stride == 1);
761 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
767 matrix_eval_EXP (gsl_matrix *m)
769 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
774 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
777 Charl Linssen <charl@itfromb.it>
781 matrix_eval_GINV (gsl_matrix *A)
786 gsl_matrix *tmp_mat = NULL;
789 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
790 tmp_mat = gsl_matrix_alloc (m, n);
791 gsl_matrix_transpose_memcpy (tmp_mat, A);
799 gsl_matrix *V = gsl_matrix_alloc (m, m);
800 gsl_vector *u = gsl_vector_alloc (m);
802 gsl_vector *tmp_vec = gsl_vector_alloc (m);
803 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
804 gsl_vector_free (tmp_vec);
807 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
808 gsl_matrix_set_zero (Sigma_pinv);
809 double cutoff = 1e-15 * gsl_vector_max (u);
811 for (size_t i = 0; i < m; ++i)
813 double x = gsl_vector_get (u, i);
814 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
817 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
818 gsl_matrix *U = gsl_matrix_calloc (n, n);
819 for (size_t i = 0; i < n; i++)
820 for (size_t j = 0; j < m; j++)
821 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
823 /* two dot products to obtain pseudoinverse */
824 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
825 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
830 A_pinv = gsl_matrix_alloc (n, m);
831 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
835 A_pinv = gsl_matrix_alloc (m, n);
836 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
839 gsl_matrix_free (tmp_mat);
840 gsl_matrix_free (tmp_mat2);
842 gsl_matrix_free (Sigma_pinv);
856 grade_compare_3way (const void *a_, const void *b_)
858 const struct grade *a = a_;
859 const struct grade *b = b_;
861 return (a->value < b->value ? -1
862 : a->value > b->value ? 1
870 matrix_eval_GRADE (gsl_matrix *m)
872 size_t n = m->size1 * m->size2;
873 struct grade *grades = xmalloc (n * sizeof *grades);
876 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
877 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
878 qsort (grades, n, sizeof *grades, grade_compare_3way);
880 for (size_t i = 0; i < n; i++)
881 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
889 dot (gsl_vector *a, gsl_vector *b)
892 for (size_t i = 0; i < a->size; i++)
893 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
898 norm2 (gsl_vector *v)
901 for (size_t i = 0; i < v->size; i++)
902 result += pow2 (gsl_vector_get (v, i));
909 return sqrt (norm2 (v));
913 matrix_eval_GSCH (gsl_matrix *v)
915 if (v->size2 < v->size1)
917 msg (SE, _("GSCH requires its argument to have at least as many columns "
918 "as rows, but it has dimensions %zu×%zu."),
922 if (!v->size1 || !v->size2)
925 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
927 for (size_t vx = 0; vx < v->size2; vx++)
929 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
930 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
932 gsl_vector_memcpy (&u_i, &v_i);
933 for (size_t j = 0; j < ux; j++)
935 gsl_vector u_j = gsl_matrix_column (u, j).vector;
936 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
937 for (size_t k = 0; k < u_i.size; k++)
938 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
939 - scale * gsl_vector_get (&u_j, k)));
942 double len = norm (&u_i);
945 gsl_vector_scale (&u_i, 1.0 / len);
946 if (++ux >= v->size1)
953 msg (SE, _("%zu×%zu argument to GSCH contains only "
954 "%zu linearly independent columns."),
955 v->size1, v->size2, ux);
965 matrix_eval_IDENT (double s1, double s2)
967 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
969 msg (SE, _("Arguments to IDENT must be non-negative integers."));
973 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
974 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
980 invert_matrix (gsl_matrix *x)
982 gsl_permutation *p = gsl_permutation_alloc (x->size1);
984 gsl_linalg_LU_decomp (x, p, &signum);
985 gsl_linalg_LU_invx (x, p);
986 gsl_permutation_free (p);
990 matrix_eval_INV (gsl_matrix *m)
997 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
999 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1000 a->size2 * b->size2);
1002 for (size_t ar = 0; ar < a->size1; ar++)
1003 for (size_t br = 0; br < b->size1; br++, y++)
1006 for (size_t ac = 0; ac < a->size2; ac++)
1007 for (size_t bc = 0; bc < b->size2; bc++, x++)
1009 double av = gsl_matrix_get (a, ar, ac);
1010 double bv = gsl_matrix_get (b, br, bc);
1011 gsl_matrix_set (k, y, x, av * bv);
1018 matrix_eval_LG10 (gsl_matrix *m)
1020 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1026 matrix_eval_LN (gsl_matrix *m)
1028 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1034 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1036 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1039 for (size_t i = 1; i <= n * n; i++)
1041 gsl_matrix_set (m, y, x, i);
1043 size_t y1 = !y ? n - 1 : y - 1;
1044 size_t x1 = x + 1 >= n ? 0 : x + 1;
1045 if (gsl_matrix_get (m, y1, x1) == 0)
1051 y = y + 1 >= n ? 0 : y + 1;
1056 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1058 double a = gsl_matrix_get (m, y1, x1);
1059 double b = gsl_matrix_get (m, y2, x2);
1060 gsl_matrix_set (m, y1, x1, b);
1061 gsl_matrix_set (m, y2, x2, a);
1065 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1069 /* A. Umar, "On the Construction of Even Order Magic Squares",
1070 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1072 for (size_t i = 1; i <= n * n / 2; i++)
1074 gsl_matrix_set (m, y, x, i);
1084 for (size_t i = n * n; i > n * n / 2; i--)
1086 gsl_matrix_set (m, y, x, i);
1094 for (size_t y = 0; y < n; y++)
1095 for (size_t x = 0; x < n / 2; x++)
1097 unsigned int d = gsl_matrix_get (m, y, x);
1098 if (d % 2 != (y < n / 2))
1099 magic_exchange (m, y, x, y, n - x - 1);
1104 size_t x1 = n / 2 - 1;
1106 magic_exchange (m, y1, x1, y2, x1);
1107 magic_exchange (m, y1, x2, y2, x2);
1111 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1113 /* A. Umar, "On the Construction of Even Order Magic Squares",
1114 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1118 for (size_t i = 1; ; i++)
1120 gsl_matrix_set (m, y, x, i);
1121 if (++y == n / 2 - 1)
1133 for (size_t i = n * n; ; i--)
1135 gsl_matrix_set (m, y, x, i);
1136 if (++y == n / 2 - 1)
1145 for (size_t y = 0; y < n; y++)
1146 if (y != n / 2 - 1 && y != n / 2)
1147 for (size_t x = 0; x < n / 2; x++)
1149 unsigned int d = gsl_matrix_get (m, y, x);
1150 if (d % 2 != (y < n / 2))
1151 magic_exchange (m, y, x, y, n - x - 1);
1154 size_t a0 = (n * n - 2 * n) / 2 + 1;
1155 for (size_t i = 0; i < n / 2; i++)
1158 gsl_matrix_set (m, n / 2, i, a);
1159 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1161 for (size_t i = 0; i < n / 2; i++)
1163 size_t a = a0 + i + n / 2;
1164 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1165 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1167 for (size_t x = 1; x < n / 2; x += 2)
1168 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1169 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1170 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1171 size_t x1 = n / 2 - 2;
1172 size_t x2 = n / 2 + 1;
1173 size_t y1 = n / 2 - 2;
1174 size_t y2 = n / 2 + 1;
1175 magic_exchange (m, y1, x1, y2, x1);
1176 magic_exchange (m, y1, x2, y2, x2);
1180 matrix_eval_MAGIC (double n_)
1182 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1184 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1189 gsl_matrix *m = gsl_matrix_calloc (n, n);
1191 matrix_eval_MAGIC_odd (m, n);
1193 matrix_eval_MAGIC_singly_even (m, n);
1195 matrix_eval_MAGIC_doubly_even (m, n);
1200 matrix_eval_MAKE (double r, double c, double value)
1202 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1204 msg (SE, _("Size arguments to MAKE must be integers."));
1208 gsl_matrix *m = gsl_matrix_alloc (r, c);
1209 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1215 matrix_eval_MDIAG (gsl_vector *v)
1217 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1218 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1219 gsl_vector_memcpy (&diagonal, v);
1224 matrix_eval_MMAX (gsl_matrix *m)
1226 return gsl_matrix_max (m);
1230 matrix_eval_MMIN (gsl_matrix *m)
1232 return gsl_matrix_min (m);
1236 matrix_eval_MOD (gsl_matrix *m, double divisor)
1240 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1244 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1245 *d = fmod (*d, divisor);
1250 matrix_eval_MSSQ (gsl_matrix *m)
1253 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1259 matrix_eval_MSUM (gsl_matrix *m)
1262 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1268 matrix_eval_NCOL (gsl_matrix *m)
1274 matrix_eval_NROW (gsl_matrix *m)
1280 matrix_eval_RANK (gsl_matrix *m)
1282 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1283 gsl_linalg_QR_decomp (m, tau);
1284 gsl_vector_free (tau);
1286 return gsl_linalg_QRPT_rank (m, -1);
1290 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1292 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1294 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1299 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1301 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1302 "product of matrix dimensions."));
1306 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1309 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1311 gsl_matrix_set (dst, y1, x1, *d);
1322 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1327 return gsl_matrix_alloc (0, 1);
1329 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1330 for (size_t y = 0; y < m->size1; y++)
1332 double ext = gsl_matrix_get (m, y, 0);
1333 for (size_t x = 1; x < m->size2; x++)
1335 double value = gsl_matrix_get (m, y, x);
1336 if (min ? value < ext : value > ext)
1339 gsl_matrix_set (rext, y, 0, ext);
1345 matrix_eval_RMAX (gsl_matrix *m)
1347 return matrix_eval_row_extremum (m, false);
1351 matrix_eval_RMIN (gsl_matrix *m)
1353 return matrix_eval_row_extremum (m, true);
1357 matrix_eval_RND (gsl_matrix *m)
1359 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1371 rank_compare_3way (const void *a_, const void *b_)
1373 const struct rank *a = a_;
1374 const struct rank *b = b_;
1376 return a->value < b->value ? -1 : a->value > b->value;
1380 matrix_eval_RNKORDER (gsl_matrix *m)
1382 size_t n = m->size1 * m->size2;
1383 struct rank *ranks = xmalloc (n * sizeof *ranks);
1385 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1386 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1387 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1389 for (size_t i = 0; i < n; )
1392 for (j = i + 1; j < n; j++)
1393 if (ranks[i].value != ranks[j].value)
1396 double rank = (i + j + 1.0) / 2.0;
1397 for (size_t k = i; k < j; k++)
1398 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1409 matrix_eval_row_sum (gsl_matrix *m, bool square)
1414 return gsl_matrix_alloc (0, 1);
1416 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1417 for (size_t y = 0; y < m->size1; y++)
1420 for (size_t x = 0; x < m->size2; x++)
1422 double d = gsl_matrix_get (m, y, x);
1423 sum += square ? pow2 (d) : d;
1425 gsl_matrix_set (result, y, 0, sum);
1431 matrix_eval_RSSQ (gsl_matrix *m)
1433 return matrix_eval_row_sum (m, true);
1437 matrix_eval_RSUM (gsl_matrix *m)
1439 return matrix_eval_row_sum (m, false);
1443 matrix_eval_SIN (gsl_matrix *m)
1445 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1451 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1453 if (m1->size1 != m2->size1)
1455 msg (SE, _("SOLVE requires its arguments to have the same number of "
1456 "rows, but the first argument has dimensions %zu×%zu and "
1457 "the second %zu×%zu."),
1458 m1->size1, m1->size2,
1459 m2->size1, m2->size2);
1463 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1464 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1466 gsl_linalg_LU_decomp (m1, p, &signum);
1467 for (size_t i = 0; i < m2->size2; i++)
1469 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1470 gsl_vector xi = gsl_matrix_column (x, i).vector;
1471 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1473 gsl_permutation_free (p);
1478 matrix_eval_SQRT (gsl_matrix *m)
1480 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1484 msg (SE, _("Argument to SQRT must be nonnegative."));
1493 matrix_eval_SSCP (gsl_matrix *m)
1495 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1496 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1501 matrix_eval_SVAL (gsl_matrix *m)
1503 gsl_matrix *tmp_mat = NULL;
1504 if (m->size2 > m->size1)
1506 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1507 gsl_matrix_transpose_memcpy (tmp_mat, m);
1512 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1513 gsl_vector *S = gsl_vector_alloc (m->size2);
1514 gsl_vector *work = gsl_vector_alloc (m->size2);
1515 gsl_linalg_SV_decomp (m, V, S, work);
1517 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1518 for (size_t i = 0; i < m->size2; i++)
1519 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1521 gsl_matrix_free (V);
1522 gsl_vector_free (S);
1523 gsl_vector_free (work);
1524 gsl_matrix_free (tmp_mat);
1530 matrix_eval_SWEEP (gsl_matrix *m, double d)
1532 if (d < 1 || d > SIZE_MAX)
1534 msg (SE, _("Scalar argument to SWEEP must be integer."));
1538 if (k >= MIN (m->size1, m->size2))
1540 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1541 "equal to the smaller of the matrix argument's rows and "
1546 double m_kk = gsl_matrix_get (m, k, k);
1547 if (fabs (m_kk) > 1e-19)
1549 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1550 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1552 double m_ij = gsl_matrix_get (m, i, j);
1553 double m_ik = gsl_matrix_get (m, i, k);
1554 double m_kj = gsl_matrix_get (m, k, j);
1555 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1564 for (size_t i = 0; i < m->size1; i++)
1566 gsl_matrix_set (m, i, k, 0);
1567 gsl_matrix_set (m, k, i, 0);
1574 matrix_eval_TRACE (gsl_matrix *m)
1577 size_t n = MIN (m->size1, m->size2);
1578 for (size_t i = 0; i < n; i++)
1579 sum += gsl_matrix_get (m, i, i);
1584 matrix_eval_T (gsl_matrix *m)
1586 return matrix_eval_TRANSPOS (m);
1590 matrix_eval_TRANSPOS (gsl_matrix *m)
1592 if (m->size1 == m->size2)
1594 gsl_matrix_transpose (m);
1599 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1600 gsl_matrix_transpose_memcpy (t, m);
1606 matrix_eval_TRUNC (gsl_matrix *m)
1608 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1614 matrix_eval_UNIFORM (double r_, double c_)
1616 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1618 msg (SE, _("Arguments to UNIFORM must be integers."));
1623 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1625 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1629 gsl_matrix *m = gsl_matrix_alloc (r, c);
1630 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1631 *d = gsl_ran_flat (get_rng (), 0, 1);
1635 struct matrix_function
1639 size_t min_args, max_args;
1642 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1644 static const struct matrix_function *
1645 matrix_parse_function_name (const char *token)
1647 static const struct matrix_function functions[] =
1649 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1653 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1655 for (size_t i = 0; i < N_FUNCTIONS; i++)
1657 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1658 return &functions[i];
1663 static struct read_file *
1664 read_file_create (struct matrix_state *s, struct file_handle *fh)
1666 for (size_t i = 0; i < s->n_read_files; i++)
1668 struct read_file *rf = s->read_files[i];
1676 struct read_file *rf = xmalloc (sizeof *rf);
1677 *rf = (struct read_file) { .file = fh };
1679 s->read_files = xrealloc (s->read_files,
1680 (s->n_read_files + 1) * sizeof *s->read_files);
1681 s->read_files[s->n_read_files++] = rf;
1685 static struct dfm_reader *
1686 read_file_open (struct read_file *rf)
1689 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1694 read_file_destroy (struct read_file *rf)
1698 fh_unref (rf->file);
1699 dfm_close_reader (rf->reader);
1700 free (rf->encoding);
1705 static struct write_file *
1706 write_file_create (struct matrix_state *s, struct file_handle *fh)
1708 for (size_t i = 0; i < s->n_write_files; i++)
1710 struct write_file *wf = s->write_files[i];
1718 struct write_file *wf = xmalloc (sizeof *wf);
1719 *wf = (struct write_file) { .file = fh };
1721 s->write_files = xrealloc (s->write_files,
1722 (s->n_write_files + 1) * sizeof *s->write_files);
1723 s->write_files[s->n_write_files++] = wf;
1727 static struct dfm_writer *
1728 write_file_open (struct write_file *wf)
1731 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1736 write_file_destroy (struct write_file *wf)
1742 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1743 wf->held->s.ss.length);
1744 u8_line_destroy (wf->held);
1748 fh_unref (wf->file);
1749 dfm_close_writer (wf->writer);
1750 free (wf->encoding);
1756 matrix_parse_function (struct matrix_state *s, const char *token,
1757 struct matrix_expr **exprp)
1760 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1763 if (lex_match_id (s->lexer, "EOF"))
1766 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1770 if (!lex_force_match (s->lexer, T_RPAREN))
1776 struct read_file *rf = read_file_create (s, fh);
1778 struct matrix_expr *e = xmalloc (sizeof *e);
1779 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1784 const struct matrix_function *f = matrix_parse_function_name (token);
1788 lex_get_n (s->lexer, 2);
1790 struct matrix_expr *e = xmalloc (sizeof *e);
1791 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1793 size_t allocated_subs = 0;
1796 struct matrix_expr *sub = matrix_parse_expr (s);
1800 if (e->n_subs >= allocated_subs)
1801 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1802 e->subs[e->n_subs++] = sub;
1804 while (lex_match (s->lexer, T_COMMA));
1805 if (!lex_force_match (s->lexer, T_RPAREN))
1812 matrix_expr_destroy (e);
1816 static struct matrix_expr *
1817 matrix_parse_primary (struct matrix_state *s)
1819 if (lex_is_number (s->lexer))
1821 double number = lex_number (s->lexer);
1823 return matrix_expr_create_number (number);
1825 else if (lex_is_string (s->lexer))
1827 char string[sizeof (double)];
1828 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1832 memcpy (&number, string, sizeof number);
1833 return matrix_expr_create_number (number);
1835 else if (lex_match (s->lexer, T_LPAREN))
1837 struct matrix_expr *e = matrix_parse_or_xor (s);
1838 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1840 matrix_expr_destroy (e);
1845 else if (lex_match (s->lexer, T_LCURLY))
1847 struct matrix_expr *e = matrix_parse_curly_semi (s);
1848 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1850 matrix_expr_destroy (e);
1855 else if (lex_token (s->lexer) == T_ID)
1857 struct matrix_expr *retval;
1858 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1861 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1864 lex_error (s->lexer, _("Unknown variable %s."),
1865 lex_tokcstr (s->lexer));
1870 struct matrix_expr *e = xmalloc (sizeof *e);
1871 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1874 else if (lex_token (s->lexer) == T_ALL)
1876 struct matrix_expr *retval;
1877 if (matrix_parse_function (s, "ALL", &retval))
1881 lex_error (s->lexer, NULL);
1885 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1888 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1890 if (lex_match (s->lexer, T_COLON))
1897 *indexp = matrix_parse_expr (s);
1898 return *indexp != NULL;
1902 static struct matrix_expr *
1903 matrix_parse_postfix (struct matrix_state *s)
1905 struct matrix_expr *lhs = matrix_parse_primary (s);
1906 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1909 struct matrix_expr *i0;
1910 if (!matrix_parse_index_expr (s, &i0))
1912 matrix_expr_destroy (lhs);
1915 if (lex_match (s->lexer, T_RPAREN))
1917 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1918 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1919 else if (lex_match (s->lexer, T_COMMA))
1921 struct matrix_expr *i1;
1922 if (!matrix_parse_index_expr (s, &i1)
1923 || !lex_force_match (s->lexer, T_RPAREN))
1925 matrix_expr_destroy (lhs);
1926 matrix_expr_destroy (i0);
1927 matrix_expr_destroy (i1);
1930 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1931 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1932 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1937 lex_error_expecting (s->lexer, "`)'", "`,'");
1942 static struct matrix_expr *
1943 matrix_parse_unary (struct matrix_state *s)
1945 if (lex_match (s->lexer, T_DASH))
1947 struct matrix_expr *lhs = matrix_parse_unary (s);
1948 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1950 else if (lex_match (s->lexer, T_PLUS))
1951 return matrix_parse_unary (s);
1953 return matrix_parse_postfix (s);
1956 static struct matrix_expr *
1957 matrix_parse_seq (struct matrix_state *s)
1959 struct matrix_expr *start = matrix_parse_unary (s);
1960 if (!start || !lex_match (s->lexer, T_COLON))
1963 struct matrix_expr *end = matrix_parse_unary (s);
1966 matrix_expr_destroy (start);
1970 if (lex_match (s->lexer, T_COLON))
1972 struct matrix_expr *increment = matrix_parse_unary (s);
1975 matrix_expr_destroy (start);
1976 matrix_expr_destroy (end);
1979 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1982 return matrix_expr_create_2 (MOP_SEQ, start, end);
1985 static struct matrix_expr *
1986 matrix_parse_exp (struct matrix_state *s)
1988 struct matrix_expr *lhs = matrix_parse_seq (s);
1995 if (lex_match (s->lexer, T_EXP))
1997 else if (lex_match_phrase (s->lexer, "&**"))
2002 struct matrix_expr *rhs = matrix_parse_seq (s);
2005 matrix_expr_destroy (lhs);
2008 lhs = matrix_expr_create_2 (op, lhs, rhs);
2012 static struct matrix_expr *
2013 matrix_parse_mul_div (struct matrix_state *s)
2015 struct matrix_expr *lhs = matrix_parse_exp (s);
2022 if (lex_match (s->lexer, T_ASTERISK))
2024 else if (lex_match (s->lexer, T_SLASH))
2026 else if (lex_match_phrase (s->lexer, "&*"))
2028 else if (lex_match_phrase (s->lexer, "&/"))
2033 struct matrix_expr *rhs = matrix_parse_exp (s);
2036 matrix_expr_destroy (lhs);
2039 lhs = matrix_expr_create_2 (op, lhs, rhs);
2043 static struct matrix_expr *
2044 matrix_parse_add_sub (struct matrix_state *s)
2046 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2053 if (lex_match (s->lexer, T_PLUS))
2055 else if (lex_match (s->lexer, T_DASH))
2057 else if (lex_token (s->lexer) == T_NEG_NUM)
2062 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2065 matrix_expr_destroy (lhs);
2068 lhs = matrix_expr_create_2 (op, lhs, rhs);
2072 static struct matrix_expr *
2073 matrix_parse_relational (struct matrix_state *s)
2075 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2082 if (lex_match (s->lexer, T_GT))
2084 else if (lex_match (s->lexer, T_GE))
2086 else if (lex_match (s->lexer, T_LT))
2088 else if (lex_match (s->lexer, T_LE))
2090 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2092 else if (lex_match (s->lexer, T_NE))
2097 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2100 matrix_expr_destroy (lhs);
2103 lhs = matrix_expr_create_2 (op, lhs, rhs);
2107 static struct matrix_expr *
2108 matrix_parse_not (struct matrix_state *s)
2110 if (lex_match (s->lexer, T_NOT))
2112 struct matrix_expr *lhs = matrix_parse_not (s);
2113 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2116 return matrix_parse_relational (s);
2119 static struct matrix_expr *
2120 matrix_parse_and (struct matrix_state *s)
2122 struct matrix_expr *lhs = matrix_parse_not (s);
2126 while (lex_match (s->lexer, T_AND))
2128 struct matrix_expr *rhs = matrix_parse_not (s);
2131 matrix_expr_destroy (lhs);
2134 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2139 static struct matrix_expr *
2140 matrix_parse_or_xor (struct matrix_state *s)
2142 struct matrix_expr *lhs = matrix_parse_and (s);
2149 if (lex_match (s->lexer, T_OR))
2151 else if (lex_match_id (s->lexer, "XOR"))
2156 struct matrix_expr *rhs = matrix_parse_and (s);
2159 matrix_expr_destroy (lhs);
2162 lhs = matrix_expr_create_2 (op, lhs, rhs);
2166 static struct matrix_expr *
2167 matrix_parse_expr (struct matrix_state *s)
2169 return matrix_parse_or_xor (s);
2172 /* Expression evaluation. */
2175 matrix_op_eval (enum matrix_op op, double a, double b)
2179 case MOP_ADD_ELEMS: return a + b;
2180 case MOP_SUB_ELEMS: return a - b;
2181 case MOP_MUL_ELEMS: return a * b;
2182 case MOP_DIV_ELEMS: return a / b;
2183 case MOP_EXP_ELEMS: return pow (a, b);
2184 case MOP_GT: return a > b;
2185 case MOP_GE: return a >= b;
2186 case MOP_LT: return a < b;
2187 case MOP_LE: return a <= b;
2188 case MOP_EQ: return a == b;
2189 case MOP_NE: return a != b;
2190 case MOP_AND: return (a > 0) && (b > 0);
2191 case MOP_OR: return (a > 0) || (b > 0);
2192 case MOP_XOR: return (a > 0) != (b > 0);
2194 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2203 case MOP_PASTE_HORZ:
2204 case MOP_PASTE_VERT:
2220 matrix_op_name (enum matrix_op op)
2224 case MOP_ADD_ELEMS: return "+";
2225 case MOP_SUB_ELEMS: return "-";
2226 case MOP_MUL_ELEMS: return "&*";
2227 case MOP_DIV_ELEMS: return "&/";
2228 case MOP_EXP_ELEMS: return "&**";
2229 case MOP_GT: return ">";
2230 case MOP_GE: return ">=";
2231 case MOP_LT: return "<";
2232 case MOP_LE: return "<=";
2233 case MOP_EQ: return "=";
2234 case MOP_NE: return "<>";
2235 case MOP_AND: return "AND";
2236 case MOP_OR: return "OR";
2237 case MOP_XOR: return "XOR";
2239 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2248 case MOP_PASTE_HORZ:
2249 case MOP_PASTE_VERT:
2265 is_scalar (const gsl_matrix *m)
2267 return m->size1 == 1 && m->size2 == 1;
2271 to_scalar (const gsl_matrix *m)
2273 assert (is_scalar (m));
2274 return gsl_matrix_get (m, 0, 0);
2278 matrix_expr_evaluate_elementwise (enum matrix_op op,
2279 gsl_matrix *a, gsl_matrix *b)
2283 double be = to_scalar (b);
2284 for (size_t r = 0; r < a->size1; r++)
2285 for (size_t c = 0; c < a->size2; c++)
2287 double *ae = gsl_matrix_ptr (a, r, c);
2288 *ae = matrix_op_eval (op, *ae, be);
2292 else if (is_scalar (a))
2294 double ae = to_scalar (a);
2295 for (size_t r = 0; r < b->size1; r++)
2296 for (size_t c = 0; c < b->size2; c++)
2298 double *be = gsl_matrix_ptr (b, r, c);
2299 *be = matrix_op_eval (op, ae, *be);
2303 else if (a->size1 == b->size1 && a->size2 == b->size2)
2305 for (size_t r = 0; r < a->size1; r++)
2306 for (size_t c = 0; c < a->size2; c++)
2308 double *ae = gsl_matrix_ptr (a, r, c);
2309 double be = gsl_matrix_get (b, r, c);
2310 *ae = matrix_op_eval (op, *ae, be);
2316 msg (SE, _("Operands to %s must have the same dimensions or one "
2317 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2318 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2324 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2326 if (is_scalar (a) || is_scalar (b))
2327 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2329 if (a->size2 != b->size1)
2331 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2332 "not conformable for multiplication."),
2333 a->size1, a->size2, b->size1, b->size2);
2337 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2338 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2343 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2345 gsl_matrix *tmp = *a;
2351 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2354 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2355 swap_matrix (z, tmp);
2359 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2361 mul_matrix (x, *x, *x, tmp);
2365 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2368 if (x->size1 != x->size2)
2370 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2371 "the left-hand size, not one with dimensions %zu×%zu."),
2372 x->size1, x->size2);
2377 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2378 "right-hand side, not a matrix with dimensions %zu×%zu."),
2379 b->size1, b->size2);
2382 double bf = to_scalar (b);
2383 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2385 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2386 "or outside the valid range."), bf);
2391 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2393 gsl_matrix_set_identity (y);
2397 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2399 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2402 mul_matrix (&y, x, y, &t);
2403 square_matrix (&x, &t);
2406 square_matrix (&x, &t);
2408 mul_matrix (&y, x, y, &t);
2412 /* Garbage collection.
2414 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2415 a permutation of them. We are returning one of them; that one must not be
2416 destroyed. We must not destroy 'x_' because the caller owns it. */
2418 gsl_matrix_free (y_);
2420 gsl_matrix_free (t_);
2426 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2429 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2431 msg (SE, _("All operands of : operator must be scalars."));
2435 long int start = to_scalar (start_);
2436 long int end = to_scalar (end_);
2437 long int by = by_ ? to_scalar (by_) : 1;
2441 msg (SE, _("The increment operand to : must be nonzero."));
2445 long int n = (end >= start && by > 0 ? (end - start + by) / by
2446 : end <= start && by < 0 ? (start - end - by) / -by
2448 gsl_matrix *m = gsl_matrix_alloc (1, n);
2449 for (long int i = 0; i < n; i++)
2450 gsl_matrix_set (m, 0, i, start + i * by);
2455 matrix_expr_evaluate_not (gsl_matrix *a)
2457 for (size_t r = 0; r < a->size1; r++)
2458 for (size_t c = 0; c < a->size2; c++)
2460 double *ae = gsl_matrix_ptr (a, r, c);
2467 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2469 if (a->size1 != b->size1)
2471 if (!a->size1 || !a->size2)
2473 else if (!b->size1 || !b->size2)
2476 msg (SE, _("All columns in a matrix must have the same number of rows, "
2477 "but this tries to paste matrices with %zu and %zu rows."),
2478 a->size1, b->size1);
2482 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2483 for (size_t y = 0; y < a->size1; y++)
2485 for (size_t x = 0; x < a->size2; x++)
2486 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2487 for (size_t x = 0; x < b->size2; x++)
2488 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2494 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2496 if (a->size2 != b->size2)
2498 if (!a->size1 || !a->size2)
2500 else if (!b->size1 || !b->size2)
2503 msg (SE, _("All rows in a matrix must have the same number of columns, "
2504 "but this tries to stack matrices with %zu and %zu columns."),
2505 a->size2, b->size2);
2509 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2510 for (size_t x = 0; x < a->size2; x++)
2512 for (size_t y = 0; y < a->size1; y++)
2513 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2514 for (size_t y = 0; y < b->size1; y++)
2515 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2521 matrix_to_vector (gsl_matrix *m)
2524 gsl_vector v = to_vector (m);
2525 assert (v.block == m->block || !v.block);
2529 return xmemdup (&v, sizeof v);
2543 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2546 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2547 enum index_type index_type, size_t other_size,
2548 struct index_vector *iv)
2557 msg (SE, _("Vector index must be scalar or vector, not a "
2559 m->size1, m->size2);
2563 msg (SE, _("Matrix row index must be scalar or vector, not a "
2565 m->size1, m->size2);
2569 msg (SE, _("Matrix column index must be scalar or vector, not a "
2571 m->size1, m->size2);
2577 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2578 *iv = (struct index_vector) {
2579 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2582 for (size_t i = 0; i < v.size; i++)
2584 double index = gsl_vector_get (&v, i);
2585 if (index < 1 || index >= size + 1)
2590 msg (SE, _("Index %g is out of range for vector "
2591 "with %zu elements."), index, size);
2595 msg (SE, _("%g is not a valid row index for "
2596 "a %zu×%zu matrix."),
2597 index, size, other_size);
2601 msg (SE, _("%g is not a valid column index for "
2602 "a %zu×%zu matrix."),
2603 index, other_size, size);
2610 iv->indexes[i] = index - 1;
2616 *iv = (struct index_vector) {
2617 .indexes = xnmalloc (size, sizeof *iv->indexes),
2620 for (size_t i = 0; i < size; i++)
2627 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2629 if (!is_vector (sm))
2631 msg (SE, _("Vector index operator may not be applied to "
2632 "a %zu×%zu matrix."),
2633 sm->size1, sm->size2);
2641 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2643 if (!matrix_expr_evaluate_vec_all (sm))
2646 gsl_vector sv = to_vector (sm);
2647 struct index_vector iv;
2648 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2651 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2652 sm->size1 == 1 ? iv.n : 1);
2653 gsl_vector dv = to_vector (dm);
2654 for (size_t dx = 0; dx < iv.n; dx++)
2656 size_t sx = iv.indexes[dx];
2657 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2665 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2668 struct index_vector iv0;
2669 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2672 struct index_vector iv1;
2673 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2680 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2681 for (size_t dy = 0; dy < iv0.n; dy++)
2683 size_t sy = iv0.indexes[dy];
2685 for (size_t dx = 0; dx < iv1.n; dx++)
2687 size_t sx = iv1.indexes[dx];
2688 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2696 #define F(NAME, PROTOTYPE) \
2697 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2698 const char *name, gsl_matrix *[], size_t, \
2699 matrix_proto_##PROTOTYPE *);
2704 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2706 if (!is_scalar (subs[index]))
2708 msg (SE, _("Function %s argument %zu must be a scalar, "
2709 "not a %zu×%zu matrix."),
2710 name, index + 1, subs[index]->size1, subs[index]->size2);
2717 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2719 if (!is_vector (subs[index]))
2721 msg (SE, _("Function %s argument %zu must be a vector, "
2722 "not a %zu×%zu matrix."),
2723 name, index + 1, subs[index]->size1, subs[index]->size2);
2730 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2732 for (size_t i = 0; i < n_subs; i++)
2734 if (!check_scalar_arg (name, subs, i))
2736 d[i] = to_scalar (subs[i]);
2742 matrix_expr_evaluate_m_d (const char *name,
2743 gsl_matrix *subs[], size_t n_subs,
2744 matrix_proto_m_d *f)
2746 assert (n_subs == 1);
2749 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2753 matrix_expr_evaluate_m_dd (const char *name,
2754 gsl_matrix *subs[], size_t n_subs,
2755 matrix_proto_m_dd *f)
2757 assert (n_subs == 2);
2760 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2764 matrix_expr_evaluate_m_ddd (const char *name,
2765 gsl_matrix *subs[], size_t n_subs,
2766 matrix_proto_m_ddd *f)
2768 assert (n_subs == 3);
2771 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2775 matrix_expr_evaluate_m_m (const char *name UNUSED,
2776 gsl_matrix *subs[], size_t n_subs,
2777 matrix_proto_m_m *f)
2779 assert (n_subs == 1);
2784 matrix_expr_evaluate_m_md (const char *name UNUSED,
2785 gsl_matrix *subs[], size_t n_subs,
2786 matrix_proto_m_md *f)
2788 assert (n_subs == 2);
2789 if (!check_scalar_arg (name, subs, 1))
2791 return f (subs[0], to_scalar (subs[1]));
2795 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2796 gsl_matrix *subs[], size_t n_subs,
2797 matrix_proto_m_mdd *f)
2799 assert (n_subs == 3);
2800 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2802 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2806 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2807 gsl_matrix *subs[], size_t n_subs,
2808 matrix_proto_m_mm *f)
2810 assert (n_subs == 2);
2811 return f (subs[0], subs[1]);
2815 matrix_expr_evaluate_m_v (const char *name,
2816 gsl_matrix *subs[], size_t n_subs,
2817 matrix_proto_m_v *f)
2819 assert (n_subs == 1);
2820 if (!check_vector_arg (name, subs, 0))
2822 gsl_vector v = to_vector (subs[0]);
2827 matrix_expr_evaluate_d_m (const char *name UNUSED,
2828 gsl_matrix *subs[], size_t n_subs,
2829 matrix_proto_d_m *f)
2831 assert (n_subs == 1);
2833 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2834 gsl_matrix_set (m, 0, 0, f (subs[0]));
2839 matrix_expr_evaluate_m_any (const char *name UNUSED,
2840 gsl_matrix *subs[], size_t n_subs,
2841 matrix_proto_m_any *f)
2843 return f (subs, n_subs);
2847 matrix_expr_evaluate_IDENT (const char *name,
2848 gsl_matrix *subs[], size_t n_subs,
2849 matrix_proto_IDENT *f)
2851 assert (n_subs <= 2);
2854 if (!to_scalar_args (name, subs, n_subs, d))
2857 return f (d[0], d[n_subs - 1]);
2861 matrix_expr_evaluate (const struct matrix_expr *e)
2863 if (e->op == MOP_NUMBER)
2865 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2866 gsl_matrix_set (m, 0, 0, e->number);
2869 else if (e->op == MOP_VARIABLE)
2871 const gsl_matrix *src = e->variable->value;
2874 msg (SE, _("Uninitialized variable %s used in expression."),
2879 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2880 gsl_matrix_memcpy (dst, src);
2883 else if (e->op == MOP_EOF)
2885 struct dfm_reader *reader = read_file_open (e->eof);
2886 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2887 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2891 enum { N_LOCAL = 3 };
2892 gsl_matrix *local_subs[N_LOCAL];
2893 gsl_matrix **subs = (e->n_subs < N_LOCAL
2895 : xmalloc (e->n_subs * sizeof *subs));
2897 for (size_t i = 0; i < e->n_subs; i++)
2899 subs[i] = matrix_expr_evaluate (e->subs[i]);
2902 for (size_t j = 0; j < i; j++)
2903 gsl_matrix_free (subs[j]);
2904 if (subs != local_subs)
2910 gsl_matrix *result = NULL;
2913 #define F(NAME, PROTOTYPE) \
2914 case MOP_F_##NAME: \
2915 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2917 matrix_eval_##NAME); \
2923 gsl_matrix_scale (subs[0], -1.0);
2941 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2945 result = matrix_expr_evaluate_not (subs[0]);
2949 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2953 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2957 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2961 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2964 case MOP_PASTE_HORZ:
2965 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2968 case MOP_PASTE_VERT:
2969 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2973 result = gsl_matrix_alloc (0, 0);
2977 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2981 result = matrix_expr_evaluate_vec_all (subs[0]);
2985 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2989 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2993 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3002 for (size_t i = 0; i < e->n_subs; i++)
3003 if (subs[i] != result)
3004 gsl_matrix_free (subs[i]);
3005 if (subs != local_subs)
3011 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3014 gsl_matrix *m = matrix_expr_evaluate (e);
3020 msg (SE, _("Expression for %s must evaluate to scalar, "
3021 "not a %zu×%zu matrix."),
3022 context, m->size1, m->size2);
3023 gsl_matrix_free (m);
3028 gsl_matrix_free (m);
3033 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3037 if (!matrix_expr_evaluate_scalar (e, context, &d))
3041 if (d < LONG_MIN || d > LONG_MAX)
3043 msg (SE, _("Expression for %s is outside the integer range."), context);
3050 struct matrix_lvalue
3052 struct matrix_var *var;
3053 struct matrix_expr *indexes[2];
3058 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3062 for (size_t i = 0; i < lvalue->n_indexes; i++)
3063 matrix_expr_destroy (lvalue->indexes[i]);
3068 static struct matrix_lvalue *
3069 matrix_lvalue_parse (struct matrix_state *s)
3071 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3072 if (!lex_force_id (s->lexer))
3074 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3075 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3079 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3084 lex_get_n (s->lexer, 2);
3086 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3088 lvalue->n_indexes++;
3090 if (lex_match (s->lexer, T_COMMA))
3092 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3094 lvalue->n_indexes++;
3096 if (!lex_force_match (s->lexer, T_RPAREN))
3102 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3108 matrix_lvalue_destroy (lvalue);
3113 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3114 enum index_type index_type, size_t other_size,
3115 struct index_vector *iv)
3120 m = matrix_expr_evaluate (e);
3127 bool ok = matrix_normalize_index_vector (m, size, index_type,
3129 gsl_matrix_free (m);
3134 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3135 struct index_vector *iv, gsl_matrix *sm)
3137 gsl_vector dv = to_vector (lvalue->var->value);
3139 /* Convert source matrix 'sm' to source vector 'sv'. */
3140 if (!is_vector (sm))
3142 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3143 sm->size1, sm->size2);
3146 gsl_vector sv = to_vector (sm);
3147 if (iv->n != sv.size)
3149 msg (SE, _("Can't assign %zu-element vector "
3150 "to %zu-element subvector."), sv.size, iv->n);
3154 /* Assign elements. */
3155 for (size_t x = 0; x < iv->n; x++)
3156 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3161 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3162 struct index_vector *iv0,
3163 struct index_vector *iv1, gsl_matrix *sm)
3165 gsl_matrix *dm = lvalue->var->value;
3167 /* Convert source matrix 'sm' to source vector 'sv'. */
3168 if (iv0->n != sm->size1)
3170 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3171 "but source matrix has %zu rows."),
3172 lvalue->var->name, iv0->n, sm->size1);
3175 if (iv1->n != sm->size2)
3177 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3178 "but source matrix has %zu columns."),
3179 lvalue->var->name, iv1->n, sm->size2);
3183 /* Assign elements. */
3184 for (size_t y = 0; y < iv0->n; y++)
3186 size_t f0 = iv0->indexes[y];
3187 for (size_t x = 0; x < iv1->n; x++)
3189 size_t f1 = iv1->indexes[x];
3190 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3196 /* Takes ownership of M. */
3198 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3199 struct index_vector *iv0, struct index_vector *iv1,
3202 if (!lvalue->n_indexes)
3204 gsl_matrix_free (lvalue->var->value);
3205 lvalue->var->value = sm;
3210 bool ok = (lvalue->n_indexes == 1
3211 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3212 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3213 free (iv0->indexes);
3214 free (iv1->indexes);
3215 gsl_matrix_free (sm);
3221 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3222 struct index_vector *iv0,
3223 struct index_vector *iv1)
3225 *iv0 = INDEX_VECTOR_INIT;
3226 *iv1 = INDEX_VECTOR_INIT;
3227 if (!lvalue->n_indexes)
3230 /* Validate destination matrix exists and has the right shape. */
3231 gsl_matrix *dm = lvalue->var->value;
3234 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3237 else if (dm->size1 == 0 || dm->size2 == 0)
3239 msg (SE, _("Cannot index %zu×%zu matrix."),
3240 dm->size1, dm->size2);
3243 else if (lvalue->n_indexes == 1)
3245 if (!is_vector (dm))
3247 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3248 dm->size1, dm->size2, lvalue->var->name);
3251 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3252 MAX (dm->size1, dm->size2),
3257 assert (lvalue->n_indexes == 2);
3258 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3259 IV_ROW, dm->size2, iv0))
3262 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3263 IV_COLUMN, dm->size1, iv1))
3265 free (iv0->indexes);
3272 /* Takes ownership of M. */
3274 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3276 struct index_vector iv0, iv1;
3277 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3279 gsl_matrix_free (sm);
3283 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3287 /* Matrix command. */
3291 struct matrix_cmd **commands;
3295 static void matrix_cmds_uninit (struct matrix_cmds *);
3299 enum matrix_cmd_type
3322 struct compute_command
3324 struct matrix_lvalue *lvalue;
3325 struct matrix_expr *rvalue;
3329 struct print_command
3331 struct matrix_expr *expression;
3332 bool use_default_format;
3333 struct fmt_spec format;
3335 int space; /* -1 means NEWPAGE. */
3339 struct string_array literals; /* CLABELS/RLABELS. */
3340 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3346 struct do_if_command
3350 struct matrix_expr *condition;
3351 struct matrix_cmds commands;
3361 struct matrix_var *var;
3362 struct matrix_expr *start;
3363 struct matrix_expr *end;
3364 struct matrix_expr *increment;
3366 /* Loop conditions. */
3367 struct matrix_expr *top_condition;
3368 struct matrix_expr *bottom_condition;
3371 struct matrix_cmds commands;
3375 struct display_command
3377 struct matrix_state *state;
3383 struct release_command
3385 struct matrix_var **vars;
3392 struct matrix_expr *expression;
3393 struct file_handle *outfile;
3394 struct string_array *variables;
3395 struct matrix_expr *names;
3396 struct stringi_set strings;
3402 struct read_file *rf;
3403 struct matrix_lvalue *dst;
3404 struct matrix_expr *size;
3406 enum fmt_type format;
3414 struct write_command
3416 struct write_file *wf;
3417 struct matrix_expr *expression;
3419 enum fmt_type format;
3423 bool hold; /* XXX */
3429 struct matrix_lvalue *dst;
3430 struct file_handle *file;
3432 struct string_array variables;
3433 struct matrix_var *names;
3435 /* Treatment of missing values. */
3440 MGET_ERROR, /* Flag error on command. */
3441 MGET_ACCEPT, /* Accept without change, user-missing only. */
3442 MGET_OMIT, /* Drop this case. */
3443 MGET_RECODE /* Recode to 'substitute'. */
3446 double substitute; /* MGET_RECODE only. */
3452 struct msave_command
3454 struct msave_common *common;
3456 struct matrix_expr *expr;
3457 const char *rowtype;
3458 struct matrix_expr *factors;
3459 struct matrix_expr *splits;
3465 struct matrix_state *state;
3466 struct file_handle *file;
3467 struct stringi_set rowtypes;
3471 struct eigen_command
3473 struct matrix_expr *expr;
3474 struct matrix_var *evec;
3475 struct matrix_var *eval;
3479 struct setdiag_command
3481 struct matrix_var *dst;
3482 struct matrix_expr *expr;
3488 struct matrix_expr *expr;
3489 struct matrix_var *u;
3490 struct matrix_var *s;
3491 struct matrix_var *v;
3497 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3498 static bool matrix_cmd_execute (struct matrix_cmd *);
3499 static void matrix_cmd_destroy (struct matrix_cmd *);
3502 static struct matrix_cmd *
3503 matrix_parse_compute (struct matrix_state *s)
3505 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3506 *cmd = (struct matrix_cmd) {
3507 .type = MCMD_COMPUTE,
3508 .compute = { .lvalue = NULL },
3511 struct compute_command *compute = &cmd->compute;
3512 compute->lvalue = matrix_lvalue_parse (s);
3513 if (!compute->lvalue)
3516 if (!lex_force_match (s->lexer, T_EQUALS))
3519 compute->rvalue = matrix_parse_expr (s);
3520 if (!compute->rvalue)
3526 matrix_cmd_destroy (cmd);
3531 matrix_cmd_execute_compute (struct compute_command *compute)
3533 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3537 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3541 print_labels_destroy (struct print_labels *labels)
3545 string_array_destroy (&labels->literals);
3546 matrix_expr_destroy (labels->expr);
3552 parse_literal_print_labels (struct matrix_state *s,
3553 struct print_labels **labelsp)
3555 lex_match (s->lexer, T_EQUALS);
3556 print_labels_destroy (*labelsp);
3557 *labelsp = xzalloc (sizeof **labelsp);
3558 while (lex_token (s->lexer) != T_SLASH
3559 && lex_token (s->lexer) != T_ENDCMD
3560 && lex_token (s->lexer) != T_STOP)
3562 struct string label = DS_EMPTY_INITIALIZER;
3563 while (lex_token (s->lexer) != T_COMMA
3564 && lex_token (s->lexer) != T_SLASH
3565 && lex_token (s->lexer) != T_ENDCMD
3566 && lex_token (s->lexer) != T_STOP)
3568 if (!ds_is_empty (&label))
3569 ds_put_byte (&label, ' ');
3571 if (lex_token (s->lexer) == T_STRING)
3572 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3575 char *rep = lex_next_representation (s->lexer, 0, 0);
3576 ds_put_cstr (&label, rep);
3581 string_array_append_nocopy (&(*labelsp)->literals,
3582 ds_steal_cstr (&label));
3584 if (!lex_match (s->lexer, T_COMMA))
3590 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3592 lex_match (s->lexer, T_EQUALS);
3593 struct matrix_expr *e = matrix_parse_exp (s);
3597 print_labels_destroy (*labelsp);
3598 *labelsp = xzalloc (sizeof **labelsp);
3599 (*labelsp)->expr = e;
3603 static struct matrix_cmd *
3604 matrix_parse_print (struct matrix_state *s)
3606 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3607 *cmd = (struct matrix_cmd) {
3610 .use_default_format = true,
3614 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3617 for (size_t i = 0; ; i++)
3619 enum token_type t = lex_next_token (s->lexer, i);
3620 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3622 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3624 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3627 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3632 cmd->print.expression = matrix_parse_exp (s);
3633 if (!cmd->print.expression)
3637 while (lex_match (s->lexer, T_SLASH))
3639 if (lex_match_id (s->lexer, "FORMAT"))
3641 lex_match (s->lexer, T_EQUALS);
3642 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3644 cmd->print.use_default_format = false;
3646 else if (lex_match_id (s->lexer, "TITLE"))
3648 lex_match (s->lexer, T_EQUALS);
3649 if (!lex_force_string (s->lexer))
3651 free (cmd->print.title);
3652 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3655 else if (lex_match_id (s->lexer, "SPACE"))
3657 lex_match (s->lexer, T_EQUALS);
3658 if (lex_match_id (s->lexer, "NEWPAGE"))
3659 cmd->print.space = -1;
3660 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3662 cmd->print.space = lex_integer (s->lexer);
3668 else if (lex_match_id (s->lexer, "RLABELS"))
3669 parse_literal_print_labels (s, &cmd->print.rlabels);
3670 else if (lex_match_id (s->lexer, "CLABELS"))
3671 parse_literal_print_labels (s, &cmd->print.clabels);
3672 else if (lex_match_id (s->lexer, "RNAMES"))
3674 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3677 else if (lex_match_id (s->lexer, "CNAMES"))
3679 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3684 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3685 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3693 matrix_cmd_destroy (cmd);
3698 matrix_is_integer (const gsl_matrix *m)
3700 for (size_t y = 0; y < m->size1; y++)
3701 for (size_t x = 0; x < m->size2; x++)
3703 double d = gsl_matrix_get (m, y, x);
3711 matrix_max_magnitude (const gsl_matrix *m)
3714 for (size_t y = 0; y < m->size1; y++)
3715 for (size_t x = 0; x < m->size2; x++)
3717 double d = fabs (gsl_matrix_get (m, y, x));
3725 format_fits (struct fmt_spec format, double x)
3727 char *s = data_out (&(union value) { .f = x }, NULL,
3728 &format, settings_get_fmt_settings ());
3729 bool fits = *s != '*' && !strchr (s, 'E');
3734 static struct fmt_spec
3735 default_format (const gsl_matrix *m, int *log_scale)
3739 double max = matrix_max_magnitude (m);
3741 if (matrix_is_integer (m))
3742 for (int w = 1; w <= 12; w++)
3744 struct fmt_spec format = { .type = FMT_F, .w = w };
3745 if (format_fits (format, -max))
3749 if (max >= 1e9 || max <= 1e-4)
3752 snprintf (s, sizeof s, "%.1e", max);
3754 const char *e = strchr (s, 'e');
3756 *log_scale = atoi (e + 1);
3759 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3763 trimmed_string (double d)
3765 struct substring s = ss_buffer ((char *) &d, sizeof d);
3766 ss_rtrim (&s, ss_cstr (" "));
3767 return ss_xstrdup (s);
3770 static struct string_array *
3771 print_labels_get (const struct print_labels *labels, size_t n,
3772 const char *prefix, bool truncate)
3777 struct string_array *out = xzalloc (sizeof *out);
3778 if (labels->literals.n)
3779 string_array_clone (out, &labels->literals);
3780 else if (labels->expr)
3782 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3783 if (m && is_vector (m))
3785 gsl_vector v = to_vector (m);
3786 for (size_t i = 0; i < v.size; i++)
3787 string_array_append_nocopy (out, trimmed_string (
3788 gsl_vector_get (&v, i)));
3790 gsl_matrix_free (m);
3796 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3799 string_array_append (out, "");
3802 string_array_delete (out, out->n - 1);
3805 for (size_t i = 0; i < out->n; i++)
3807 char *s = out->strings[i];
3808 s[strnlen (s, 8)] = '\0';
3815 matrix_cmd_print_space (int space)
3818 output_item_submit (page_break_item_create ());
3819 for (int i = 0; i < space; i++)
3820 output_log ("%s", "");
3824 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3825 struct fmt_spec format, int log_scale)
3827 matrix_cmd_print_space (print->space);
3829 output_log ("%s", print->title);
3831 output_log (" 10 ** %d X", log_scale);
3833 struct string_array *clabels = print_labels_get (print->clabels,
3834 m->size2, "col", true);
3835 if (clabels && format.w < 8)
3837 struct string_array *rlabels = print_labels_get (print->rlabels,
3838 m->size1, "row", true);
3842 struct string line = DS_EMPTY_INITIALIZER;
3844 ds_put_byte_multiple (&line, ' ', 8);
3845 for (size_t x = 0; x < m->size2; x++)
3846 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3847 output_log_nocopy (ds_steal_cstr (&line));
3850 double scale = pow (10.0, log_scale);
3851 bool numeric = fmt_is_numeric (format.type);
3852 for (size_t y = 0; y < m->size1; y++)
3854 struct string line = DS_EMPTY_INITIALIZER;
3856 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3858 for (size_t x = 0; x < m->size2; x++)
3860 double f = gsl_matrix_get (m, y, x);
3862 ? data_out (&(union value) { .f = f / scale}, NULL,
3863 &format, settings_get_fmt_settings ())
3864 : trimmed_string (f));
3865 ds_put_format (&line, " %s", s);
3868 output_log_nocopy (ds_steal_cstr (&line));
3871 string_array_destroy (rlabels);
3873 string_array_destroy (clabels);
3878 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3879 const struct print_labels *print_labels, size_t n,
3880 const char *name, const char *prefix)
3882 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3884 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3885 for (size_t i = 0; i < n; i++)
3886 pivot_category_create_leaf (
3888 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3889 : pivot_value_new_integer (i + 1)));
3891 d->hide_all_labels = true;
3892 string_array_destroy (labels);
3897 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3898 struct fmt_spec format, int log_scale)
3900 struct pivot_table *table = pivot_table_create__ (
3901 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3904 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3906 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3907 N_("Columns"), "col");
3909 struct pivot_footnote *footnote = NULL;
3912 char *s = xasprintf ("× 10**%d", log_scale);
3913 footnote = pivot_table_create_footnote (
3914 table, pivot_value_new_user_text_nocopy (s));
3917 double scale = pow (10.0, log_scale);
3918 bool numeric = fmt_is_numeric (format.type);
3919 for (size_t y = 0; y < m->size1; y++)
3920 for (size_t x = 0; x < m->size2; x++)
3922 double f = gsl_matrix_get (m, y, x);
3923 struct pivot_value *v;
3926 v = pivot_value_new_number (f / scale);
3927 v->numeric.format = format;
3930 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3932 pivot_value_add_footnote (v, footnote);
3933 pivot_table_put2 (table, y, x, v);
3936 pivot_table_submit (table);
3940 matrix_cmd_execute_print (const struct print_command *print)
3942 if (print->expression)
3944 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3949 struct fmt_spec format = (print->use_default_format
3950 ? default_format (m, &log_scale)
3953 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3954 matrix_cmd_print_text (print, m, format, log_scale);
3956 matrix_cmd_print_tables (print, m, format, log_scale);
3958 gsl_matrix_free (m);
3962 matrix_cmd_print_space (print->space);
3964 output_log ("%s", print->title);
3971 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3972 const char *command_name,
3973 const char *stop1, const char *stop2)
3975 lex_end_of_command (s->lexer);
3976 lex_discard_rest_of_command (s->lexer);
3978 size_t allocated = 0;
3981 while (lex_token (s->lexer) == T_ENDCMD)
3984 if (lex_at_phrase (s->lexer, stop1)
3985 || (stop2 && lex_at_phrase (s->lexer, stop2)))
3988 if (lex_at_phrase (s->lexer, "END MATRIX"))
3990 msg (SE, _("Premature END MATRIX within %s."), command_name);
3994 struct matrix_cmd *cmd = matrix_parse_command (s);
3998 if (c->n >= allocated)
3999 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4000 c->commands[c->n++] = cmd;
4005 matrix_parse_do_if_clause (struct matrix_state *s,
4006 struct do_if_command *ifc,
4007 bool parse_condition,
4008 size_t *allocated_clauses)
4010 if (ifc->n_clauses >= *allocated_clauses)
4011 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4012 sizeof *ifc->clauses);
4013 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4014 *c = (struct do_if_clause) { .condition = NULL };
4016 if (parse_condition)
4018 c->condition = matrix_parse_expr (s);
4023 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4026 static struct matrix_cmd *
4027 matrix_parse_do_if (struct matrix_state *s)
4029 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4030 *cmd = (struct matrix_cmd) {
4032 .do_if = { .n_clauses = 0 }
4035 size_t allocated_clauses = 0;
4038 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4041 while (lex_match_phrase (s->lexer, "ELSE IF"));
4043 if (lex_match_id (s->lexer, "ELSE")
4044 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4047 if (!lex_match_phrase (s->lexer, "END IF"))
4052 matrix_cmd_destroy (cmd);
4057 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4059 for (size_t i = 0; i < cmd->n_clauses; i++)
4061 struct do_if_clause *c = &cmd->clauses[i];
4065 if (!matrix_expr_evaluate_scalar (c->condition,
4066 i ? "ELSE IF" : "DO IF",
4071 for (size_t j = 0; j < c->commands.n; j++)
4072 if (!matrix_cmd_execute (c->commands.commands[j]))
4079 static struct matrix_cmd *
4080 matrix_parse_loop (struct matrix_state *s)
4082 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4083 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4085 struct loop_command *loop = &cmd->loop;
4086 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4088 struct substring name = lex_tokss (s->lexer);
4089 loop->var = matrix_var_lookup (s, name);
4091 loop->var = matrix_var_create (s, name);
4096 loop->start = matrix_parse_expr (s);
4097 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4100 loop->end = matrix_parse_expr (s);
4104 if (lex_match (s->lexer, T_BY))
4106 loop->increment = matrix_parse_expr (s);
4107 if (!loop->increment)
4112 if (lex_match_id (s->lexer, "IF"))
4114 loop->top_condition = matrix_parse_expr (s);
4115 if (!loop->top_condition)
4119 bool was_in_loop = s->in_loop;
4121 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4123 s->in_loop = was_in_loop;
4127 if (!lex_match_phrase (s->lexer, "END LOOP"))
4130 if (lex_match_id (s->lexer, "IF"))
4132 loop->bottom_condition = matrix_parse_expr (s);
4133 if (!loop->bottom_condition)
4140 matrix_cmd_destroy (cmd);
4145 matrix_cmd_execute_loop (struct loop_command *cmd)
4149 long int increment = 1;
4152 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4153 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4155 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4159 if (increment > 0 ? value > end
4160 : increment < 0 ? value < end
4165 int mxloops = settings_get_mxloops ();
4166 for (int i = 0; i < mxloops; i++)
4171 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4173 gsl_matrix_free (cmd->var->value);
4174 cmd->var->value = NULL;
4176 if (!cmd->var->value)
4177 cmd->var->value = gsl_matrix_alloc (1, 1);
4179 gsl_matrix_set (cmd->var->value, 0, 0, value);
4182 if (cmd->top_condition)
4185 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4190 for (size_t j = 0; j < cmd->commands.n; j++)
4191 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4194 if (cmd->bottom_condition)
4197 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4198 "END LOOP IF", &d) || d > 0)
4205 if (increment > 0 ? value > end : value < end)
4211 static struct matrix_cmd *
4212 matrix_parse_break (struct matrix_state *s)
4216 msg (SE, _("BREAK not inside LOOP."));
4220 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4221 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4225 static struct matrix_cmd *
4226 matrix_parse_display (struct matrix_state *s)
4228 bool dictionary = false;
4229 bool status = false;
4230 while (lex_token (s->lexer) == T_ID)
4232 if (lex_match_id (s->lexer, "DICTIONARY"))
4234 else if (lex_match_id (s->lexer, "STATUS"))
4238 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4242 if (!dictionary && !status)
4243 dictionary = status = true;
4245 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4246 *cmd = (struct matrix_cmd) {
4247 .type = MCMD_DISPLAY,
4250 .dictionary = dictionary,
4258 compare_matrix_var_pointers (const void *a_, const void *b_)
4260 const struct matrix_var *const *ap = a_;
4261 const struct matrix_var *const *bp = b_;
4262 const struct matrix_var *a = *ap;
4263 const struct matrix_var *b = *bp;
4264 return strcmp (a->name, b->name);
4268 matrix_cmd_execute_display (struct display_command *cmd)
4270 const struct matrix_state *s = cmd->state;
4272 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4273 pivot_dimension_create (
4274 table, PIVOT_AXIS_COLUMN, N_("Property"),
4275 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4277 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4280 const struct matrix_var *var;
4281 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4283 vars[n_vars++] = var;
4284 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4286 struct pivot_dimension *rows = pivot_dimension_create (
4287 table, PIVOT_AXIS_ROW, N_("Variable"));
4288 for (size_t i = 0; i < n_vars; i++)
4290 const struct matrix_var *var = vars[i];
4291 pivot_category_create_leaf (
4292 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4294 size_t r = var->value->size1;
4295 size_t c = var->value->size2;
4296 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4297 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4298 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4300 pivot_table_submit (table);
4303 static struct matrix_cmd *
4304 matrix_parse_release (struct matrix_state *s)
4306 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4307 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4309 size_t allocated_vars = 0;
4310 while (lex_token (s->lexer) == T_ID)
4312 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4315 if (cmd->release.n_vars >= allocated_vars)
4316 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4317 sizeof *cmd->release.vars);
4318 cmd->release.vars[cmd->release.n_vars++] = var;
4321 lex_error (s->lexer, _("Variable name expected."));
4324 if (!lex_match (s->lexer, T_COMMA))
4332 matrix_cmd_execute_release (struct release_command *cmd)
4334 for (size_t i = 0; i < cmd->n_vars; i++)
4336 struct matrix_var *var = cmd->vars[i];
4337 gsl_matrix_free (var->value);
4342 static struct matrix_cmd *
4343 matrix_parse_save (struct matrix_state *s)
4345 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4346 *cmd = (struct matrix_cmd) {
4349 .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4353 struct save_command *save = &cmd->save;
4354 save->expression = matrix_parse_exp (s);
4355 if (!save->expression)
4358 while (lex_match (s->lexer, T_SLASH))
4360 if (lex_match_id (s->lexer, "OUTFILE"))
4362 lex_match (s->lexer, T_EQUALS);
4364 fh_unref (save->outfile);
4365 save->outfile = (lex_match (s->lexer, T_ASTERISK)
4367 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4371 else if (lex_match_id (s->lexer, "VARIABLES"))
4373 lex_match (s->lexer, T_EQUALS);
4377 struct dictionary *d = dict_create (get_default_encoding ());
4378 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4379 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4384 string_array_destroy (save->variables);
4385 if (!save->variables)
4386 save->variables = xmalloc (sizeof *save->variables);
4387 *save->variables = (struct string_array) {
4393 else if (lex_match_id (s->lexer, "NAMES"))
4395 lex_match (s->lexer, T_EQUALS);
4396 matrix_expr_destroy (save->names);
4397 save->names = matrix_parse_exp (s);
4401 else if (lex_match_id (s->lexer, "STRINGS"))
4403 lex_match (s->lexer, T_EQUALS);
4404 while (lex_token (s->lexer) == T_ID)
4406 stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4408 if (!lex_match (s->lexer, T_COMMA))
4414 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4422 if (s->prev_save_outfile)
4423 save->outfile = fh_ref (s->prev_save_outfile);
4426 lex_sbc_missing ("OUTFILE");
4430 fh_unref (s->prev_save_outfile);
4431 s->prev_save_outfile = fh_ref (save->outfile);
4433 if (save->variables && save->names)
4435 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4436 matrix_expr_destroy (save->names);
4443 matrix_cmd_destroy (cmd);
4448 matrix_cmd_execute_save (const struct save_command *save)
4450 assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4452 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4457 struct dictionary *dict = dict_create (get_default_encoding ());
4458 struct string_array names = { .n = 0 };
4459 if (save->variables)
4460 string_array_clone (&names, save->variables);
4461 else if (save->names)
4463 gsl_matrix *nm = matrix_expr_evaluate (save->names);
4464 if (nm && is_vector (nm))
4466 gsl_vector nv = to_vector (nm);
4467 for (size_t i = 0; i < nv.size; i++)
4469 char *name = trimmed_string (gsl_vector_get (&nv, i));
4470 if (dict_id_is_valid (dict, name, true))
4471 string_array_append_nocopy (&names, name);
4478 struct stringi_set strings;
4479 stringi_set_clone (&strings, &save->strings);
4481 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4486 name = names.strings[i];
4489 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4493 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4494 struct variable *var = dict_create_var (dict, name, width);
4497 msg (SE, _("Duplicate variable name %s in SAVE statement."),
4503 if (!stringi_set_is_empty (&strings))
4505 const char *example = stringi_set_node_get_string (
4506 stringi_set_first (&strings));
4507 msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4508 "with that name was found on SAVE.",
4509 "STRINGS specified %2$zu variables, including %1$s, "
4510 "whose names were not found on SAVE.",
4511 stringi_set_count (&strings)),
4512 example, stringi_set_count (&strings));
4515 stringi_set_destroy (&strings);
4516 string_array_destroy (&names);
4520 gsl_matrix_free (m);
4525 struct casewriter *writer = any_writer_open (save->outfile, dict);
4528 gsl_matrix_free (m);
4533 const struct caseproto *proto = dict_get_proto (dict);
4534 for (size_t y = 0; y < m->size1; y++)
4536 struct ccase *c = case_create (proto);
4537 for (size_t x = 0; x < m->size2; x++)
4539 double d = gsl_matrix_get (m, y, x);
4540 int width = caseproto_get_width (proto, x);
4541 union value *value = case_data_rw_idx (c, x);
4545 memcpy (value->s, &d, width);
4547 casewriter_write (writer, c);
4549 casewriter_destroy (writer);
4551 gsl_matrix_free (m);
4555 static struct matrix_cmd *
4556 matrix_parse_read (struct matrix_state *s)
4558 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4559 *cmd = (struct matrix_cmd) {
4561 .read = { .format = FMT_F },
4564 struct file_handle *fh = NULL;
4565 char *encoding = NULL;
4566 struct read_command *read = &cmd->read;
4567 read->dst = matrix_lvalue_parse (s);
4572 int repetitions = 0;
4573 int record_width = 0;
4574 bool seen_format = false;
4575 while (lex_match (s->lexer, T_SLASH))
4577 if (lex_match_id (s->lexer, "FILE"))
4579 lex_match (s->lexer, T_EQUALS);
4582 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4586 else if (lex_match_id (s->lexer, "ENCODING"))
4588 lex_match (s->lexer, T_EQUALS);
4589 if (!lex_force_string (s->lexer))
4593 encoding = ss_xstrdup (lex_tokss (s->lexer));
4597 else if (lex_match_id (s->lexer, "FIELD"))
4599 lex_match (s->lexer, T_EQUALS);
4601 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4603 read->c1 = lex_integer (s->lexer);
4605 if (!lex_force_match (s->lexer, T_TO)
4606 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4608 read->c2 = lex_integer (s->lexer) + 1;
4611 record_width = read->c2 - read->c1;
4612 if (lex_match (s->lexer, T_BY))
4614 if (!lex_force_int_range (s->lexer, "BY", 1,
4615 read->c2 - read->c1))
4617 by = lex_integer (s->lexer);
4620 if (record_width % by)
4622 msg (SE, _("BY %d does not evenly divide record width %d."),
4630 else if (lex_match_id (s->lexer, "SIZE"))
4632 lex_match (s->lexer, T_EQUALS);
4633 matrix_expr_destroy (read->size);
4634 read->size = matrix_parse_exp (s);
4638 else if (lex_match_id (s->lexer, "MODE"))
4640 lex_match (s->lexer, T_EQUALS);
4641 if (lex_match_id (s->lexer, "RECTANGULAR"))
4642 read->symmetric = false;
4643 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4644 read->symmetric = true;
4647 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4651 else if (lex_match_id (s->lexer, "REREAD"))
4652 read->reread = true;
4653 else if (lex_match_id (s->lexer, "FORMAT"))
4657 lex_sbc_only_once ("FORMAT");
4662 lex_match (s->lexer, T_EQUALS);
4664 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4667 const char *p = lex_tokcstr (s->lexer);
4668 if (c_isdigit (p[0]))
4670 repetitions = atoi (p);
4671 p += strspn (p, "0123456789");
4672 if (!fmt_from_name (p, &read->format))
4674 lex_error (s->lexer, _("Unknown format %s."), p);
4679 else if (!fmt_from_name (p, &read->format))
4681 struct fmt_spec format;
4682 if (!parse_format_specifier (s->lexer, &format))
4684 read->format = format.type;
4691 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4692 "REREAD", "FORMAT");
4699 lex_sbc_missing ("FIELD");
4703 if (!read->dst->n_indexes && !read->size)
4705 msg (SE, _("SIZE is required for reading data into a full matrix "
4706 "(as opposed to a submatrix)."));
4712 if (s->prev_read_file)
4713 fh = fh_ref (s->prev_read_file);
4716 lex_sbc_missing ("FILE");
4720 fh_unref (s->prev_read_file);
4721 s->prev_read_file = fh_ref (fh);
4723 read->rf = read_file_create (s, fh);
4727 free (read->rf->encoding);
4728 read->rf->encoding = encoding;
4732 /* Field width may be specified in multiple ways:
4735 2. The format on FORMAT.
4736 3. The repetition factor on FORMAT.
4738 (2) and (3) are mutually exclusive.
4740 If more than one of these is present, they must agree. If none of them is
4741 present, then free-field format is used.
4743 if (repetitions > record_width)
4745 msg (SE, _("%d repetitions cannot fit in record width %d."),
4746 repetitions, record_width);
4749 int w = (repetitions ? record_width / repetitions
4755 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4756 "which implies field width %d, "
4757 "but BY specifies field width %d."),
4758 repetitions, record_width, w, by);
4760 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4769 matrix_cmd_destroy (cmd);
4775 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4776 gsl_matrix *m, struct substring p, size_t y, size_t x)
4778 const char *input_encoding = dfm_reader_get_encoding (reader);
4780 char *error = data_in (p, input_encoding, read->format,
4781 settings_get_fmt_settings (), &v, 0, NULL);
4782 /* XXX report error if value is missing */
4784 msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4787 gsl_matrix_set (m, y, x, v.f);
4788 if (read->symmetric && x != y)
4789 gsl_matrix_set (m, x, y, v.f);
4794 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4795 struct substring *line)
4797 if (dfm_eof (reader))
4799 msg (SE, _("Unexpected end of file reading matrix data."));
4802 dfm_expand_tabs (reader);
4803 *line = ss_substr (dfm_get_record (reader),
4804 read->c1 - 1, read->c2 - read->c1);
4809 matrix_read (struct read_command *read, struct dfm_reader *reader,
4812 for (size_t y = 0; y < m->size1; y++)
4814 size_t nx = read->symmetric ? y + 1 : m->size2;
4816 struct substring line = ss_empty ();
4817 for (size_t x = 0; x < nx; x++)
4824 ss_ltrim (&line, ss_cstr (" ,"));
4825 if (!ss_is_empty (line))
4827 if (!matrix_read_line (read, reader, &line))
4829 dfm_forward_record (reader);
4832 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4836 if (!matrix_read_line (read, reader, &line))
4838 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4839 int f = x % fields_per_line;
4840 if (f == fields_per_line - 1)
4841 dfm_forward_record (reader);
4843 p = ss_substr (line, read->w * f, read->w);
4846 matrix_read_set_field (read, reader, m, p, y, x);
4850 dfm_forward_record (reader);
4853 ss_ltrim (&line, ss_cstr (" ,"));
4854 if (!ss_is_empty (line))
4855 msg (SW, _("Trailing garbage on line \"%.*s\""),
4856 (int) line.length, line.string);
4862 matrix_cmd_execute_read (struct read_command *read)
4864 struct index_vector iv0, iv1;
4865 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4868 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4871 gsl_matrix *m = matrix_expr_evaluate (read->size);
4877 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4878 "not a %zu×%zu matrix."), m->size1, m->size2);
4879 gsl_matrix_free (m);
4885 gsl_vector v = to_vector (m);
4889 d[0] = gsl_vector_get (&v, 0);
4892 else if (v.size == 2)
4894 d[0] = gsl_vector_get (&v, 0);
4895 d[1] = gsl_vector_get (&v, 1);
4899 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4900 "not a %zu×%zu matrix."),
4901 m->size1, m->size2),
4902 gsl_matrix_free (m);
4907 gsl_matrix_free (m);
4909 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4911 msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4921 if (read->dst->n_indexes)
4923 size_t submatrix_size[2];
4924 if (read->dst->n_indexes == 2)
4926 submatrix_size[0] = iv0.n;
4927 submatrix_size[1] = iv1.n;
4929 else if (read->dst->var->value->size1 == 1)
4931 submatrix_size[0] = 1;
4932 submatrix_size[1] = iv0.n;
4936 submatrix_size[0] = iv0.n;
4937 submatrix_size[1] = 1;
4942 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4944 msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4947 submatrix_size[0], submatrix_size[1]);
4955 size[0] = submatrix_size[0];
4956 size[1] = submatrix_size[1];
4960 struct dfm_reader *reader = read_file_open (read->rf);
4962 dfm_reread_record (reader, 1);
4964 if (read->symmetric && size[0] != size[1])
4966 msg (SE, _("Cannot read non-square %zu×%zu matrix "
4967 "using READ with MODE=SYMMETRIC."),
4973 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4974 matrix_read (read, reader, tmp);
4975 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4978 static struct matrix_cmd *
4979 matrix_parse_write (struct matrix_state *s)
4981 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4982 *cmd = (struct matrix_cmd) {
4984 .write = { .format = FMT_F },
4987 struct file_handle *fh = NULL;
4988 char *encoding = NULL;
4989 struct write_command *write = &cmd->write;
4990 write->expression = matrix_parse_exp (s);
4991 if (!write->expression)
4995 int repetitions = 0;
4996 int record_width = 0;
4997 bool seen_format = false;
4998 while (lex_match (s->lexer, T_SLASH))
5000 if (lex_match_id (s->lexer, "OUTFILE"))
5002 lex_match (s->lexer, T_EQUALS);
5005 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5009 else if (lex_match_id (s->lexer, "ENCODING"))
5011 lex_match (s->lexer, T_EQUALS);
5012 if (!lex_force_string (s->lexer))
5016 encoding = ss_xstrdup (lex_tokss (s->lexer));
5020 else if (lex_match_id (s->lexer, "FIELD"))
5022 lex_match (s->lexer, T_EQUALS);
5024 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5026 write->c1 = lex_integer (s->lexer);
5028 if (!lex_force_match (s->lexer, T_TO)
5029 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5031 write->c2 = lex_integer (s->lexer) + 1;
5034 record_width = write->c2 - write->c1;
5035 if (lex_match (s->lexer, T_BY))
5037 if (!lex_force_int_range (s->lexer, "BY", 1,
5038 write->c2 - write->c1))
5040 by = lex_integer (s->lexer);
5043 if (record_width % by)
5045 msg (SE, _("BY %d does not evenly divide record width %d."),
5053 else if (lex_match_id (s->lexer, "MODE"))
5055 lex_match (s->lexer, T_EQUALS);
5056 if (lex_match_id (s->lexer, "RECTANGULAR"))
5057 write->triangular = false;
5058 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5059 write->triangular = true;
5062 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5066 else if (lex_match_id (s->lexer, "HOLD"))
5068 else if (lex_match_id (s->lexer, "FORMAT"))
5072 lex_sbc_only_once ("FORMAT");
5077 lex_match (s->lexer, T_EQUALS);
5079 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5082 const char *p = lex_tokcstr (s->lexer);
5083 if (c_isdigit (p[0]))
5085 repetitions = atoi (p);
5086 p += strspn (p, "0123456789");
5087 if (!fmt_from_name (p, &write->format))
5089 lex_error (s->lexer, _("Unknown format %s."), p);
5094 else if (!fmt_from_name (p, &write->format))
5096 struct fmt_spec format;
5097 if (!parse_format_specifier (s->lexer, &format))
5099 write->format = format.type;
5100 write->w = format.w;
5101 write->d = format.d;
5106 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5114 lex_sbc_missing ("FIELD");
5120 if (s->prev_write_file)
5121 fh = fh_ref (s->prev_write_file);
5124 lex_sbc_missing ("OUTFILE");
5128 fh_unref (s->prev_write_file);
5129 s->prev_write_file = fh_ref (fh);
5131 write->wf = write_file_create (s, fh);
5135 free (write->wf->encoding);
5136 write->wf->encoding = encoding;
5140 /* Field width may be specified in multiple ways:
5143 2. The format on FORMAT.
5144 3. The repetition factor on FORMAT.
5146 (2) and (3) are mutually exclusive.
5148 If more than one of these is present, they must agree. If none of them is
5149 present, then free-field format is used.
5151 if (repetitions > record_width)
5153 msg (SE, _("%d repetitions cannot fit in record width %d."),
5154 repetitions, record_width);
5157 int w = (repetitions ? record_width / repetitions
5158 : write->w ? write->w
5163 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5164 "which implies field width %d, "
5165 "but BY specifies field width %d."),
5166 repetitions, record_width, w, by);
5168 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5177 matrix_cmd_destroy (cmd);
5182 matrix_cmd_execute_write (struct write_command *write)
5184 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5188 if (write->triangular && m->size1 != m->size2)
5190 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5191 "the matrix to be written has dimensions %zu×%zu."),
5192 m->size1, m->size2);
5193 gsl_matrix_free (m);
5197 struct dfm_writer *writer = write_file_open (write->wf);
5198 if (!writer || !m->size1)
5200 gsl_matrix_free (m);
5204 const struct fmt_settings *settings = settings_get_fmt_settings ();
5205 struct fmt_spec format = {
5206 .type = write->format,
5207 .w = write->w ? write->w : 40,
5210 struct u8_line *line = write->wf->held;
5211 for (size_t y = 0; y < m->size1; y++)
5215 line = xmalloc (sizeof *line);
5216 u8_line_init (line);
5218 size_t nx = write->triangular ? y + 1 : m->size2;
5220 for (size_t x = 0; x < nx; x++)
5222 /* XXX string values */
5223 union value v = { .f = gsl_matrix_get (m, y, x) };
5225 ? data_out (&v, NULL, &format, settings)
5226 : data_out_stretchy (&v, NULL, &format, settings, NULL));
5227 size_t len = strlen (s);
5228 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5229 if (width + x0 > write->c2)
5231 dfm_put_record_utf8 (writer, line->s.ss.string,
5233 u8_line_clear (line);
5236 u8_line_put (line, x0, x0 + width, s, len);
5239 x0 += write->w ? write->w : width + 1;
5242 if (y + 1 >= m->size1 && write->hold)
5244 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5245 u8_line_clear (line);
5249 u8_line_destroy (line);
5252 write->wf->held = line;
5253 dfm_close_writer (writer);
5255 gsl_matrix_free (m);
5258 static struct matrix_cmd *
5259 matrix_parse_get (struct matrix_state *s)
5261 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5262 *cmd = (struct matrix_cmd) {
5265 .user = { .treatment = MGET_ERROR },
5266 .system = { .treatment = MGET_ERROR },
5270 struct get_command *get = &cmd->get;
5271 get->dst = matrix_lvalue_parse (s);
5275 while (lex_match (s->lexer, T_SLASH))
5277 if (lex_match_id (s->lexer, "FILE"))
5279 if (get->variables.n)
5281 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5284 lex_match (s->lexer, T_EQUALS);
5286 fh_unref (get->file);
5287 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5291 else if (lex_match_id (s->lexer, "ENCODING"))
5293 if (get->variables.n)
5295 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5298 lex_match (s->lexer, T_EQUALS);
5299 if (!lex_force_string (s->lexer))
5302 free (get->encoding);
5303 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5307 else if (lex_match_id (s->lexer, "VARIABLES"))
5309 lex_match (s->lexer, T_EQUALS);
5311 struct dictionary *dict = NULL;
5314 dict = dataset_dict (s->dataset);
5315 if (dict_get_var_cnt (dict) == 0)
5317 lex_error (s->lexer, _("GET cannot read empty active file."));
5323 struct casereader *reader = any_reader_open_and_decode (
5324 get->file, get->encoding, &dict, NULL);
5327 casereader_destroy (reader);
5330 struct variable **vars;
5332 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5333 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5340 string_array_clear (&get->variables);
5341 for (size_t i = 0; i < n_vars; i++)
5342 string_array_append (&get->variables, var_get_name (vars[i]));
5346 else if (lex_match_id (s->lexer, "NAMES"))
5348 lex_match (s->lexer, T_EQUALS);
5349 if (!lex_force_id (s->lexer))
5352 struct substring name = lex_tokss (s->lexer);
5353 get->names = matrix_var_lookup (s, name);
5355 get->names = matrix_var_create (s, name);
5358 else if (lex_match_id (s->lexer, "MISSING"))
5360 lex_match (s->lexer, T_EQUALS);
5361 if (lex_match_id (s->lexer, "ACCEPT"))
5362 get->user.treatment = MGET_ACCEPT;
5363 else if (lex_match_id (s->lexer, "OMIT"))
5364 get->user.treatment = MGET_OMIT;
5365 else if (lex_is_number (s->lexer))
5367 get->user.treatment = MGET_RECODE;
5368 get->user.substitute = lex_number (s->lexer);
5373 lex_error (s->lexer, NULL);
5377 else if (lex_match_id (s->lexer, "SYSMIS"))
5379 lex_match (s->lexer, T_EQUALS);
5380 if (lex_match_id (s->lexer, "OMIT"))
5381 get->user.treatment = MGET_OMIT;
5382 else if (lex_is_number (s->lexer))
5384 get->user.treatment = MGET_RECODE;
5385 get->user.substitute = lex_number (s->lexer);
5390 lex_error (s->lexer, NULL);
5396 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5397 "MISSING", "SYSMIS");
5404 matrix_cmd_destroy (cmd);
5409 matrix_cmd_execute_get (struct get_command *get)
5411 assert (get->file); /* XXX */
5413 struct dictionary *dict;
5414 struct casereader *reader = any_reader_open_and_decode (
5415 get->file, get->encoding, &dict, NULL);
5419 const struct variable **vars = xnmalloc (
5420 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5424 if (get->variables.n)
5426 for (size_t i = 0; i < get->variables.n; i++)
5428 const char *name = get->variables.strings[i];
5429 const struct variable *var = dict_lookup_var (dict, name);
5432 msg (SE, _("GET: Data file does not contain variable %s."),
5438 if (!var_is_numeric (var))
5440 msg (SE, _("GET: Variable %s is not numeric."), name);
5445 vars[n_vars++] = var;
5450 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5452 const struct variable *var = dict_get_var (dict, i);
5453 if (!var_is_numeric (var))
5455 msg (SE, _("GET: Variable %s is not numeric."),
5456 var_get_name (var));
5461 vars[n_vars++] = var;
5466 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5467 long long int casenum = 1;
5469 for (struct ccase *c = casereader_read (reader); c;
5470 c = casereader_read (reader), casenum++)
5472 if (n_rows >= m->size1)
5474 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5475 for (size_t y = 0; y < n_rows; y++)
5476 for (size_t x = 0; x < n_vars; x++)
5477 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5478 gsl_matrix_free (m);
5483 for (size_t x = 0; x < n_vars; x++)
5485 const struct variable *var = vars[x];
5486 double d = case_num (c, var);
5489 if (get->system.treatment == MGET_RECODE)
5490 d = get->system.substitute;
5491 else if (get->system.treatment == MGET_OMIT)
5495 msg (SE, _("GET: Variable %s in case %lld "
5496 "is system-missing."),
5497 var_get_name (var), casenum);
5501 else if (var_is_num_missing (var, d, MV_USER))
5503 if (get->user.treatment == MGET_RECODE)
5504 d = get->user.substitute;
5505 else if (get->user.treatment == MGET_OMIT)
5507 else if (get->user.treatment != MGET_ACCEPT)
5509 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5511 var_get_name (var), casenum, d);
5515 gsl_matrix_set (m, n_rows, x, d);
5523 casereader_destroy (reader);
5527 matrix_lvalue_evaluate_and_assign (get->dst, m);
5530 gsl_matrix_free (m);
5536 match_rowtype (struct lexer *lexer)
5538 static const char *rowtypes[] = {
5539 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5541 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5543 for (size_t i = 0; i < n_rowtypes; i++)
5544 if (lex_match_id (lexer, rowtypes[i]))
5547 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5552 parse_var_names (struct lexer *lexer, struct string_array *sa)
5554 lex_match (lexer, T_EQUALS);
5556 string_array_clear (sa);
5558 struct dictionary *dict = dict_create (get_default_encoding ());
5559 dict_create_var_assert (dict, "ROWTYPE_", 8);
5560 dict_create_var_assert (dict, "VARNAME_", 8);
5563 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5564 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5569 string_array_clear (sa);
5570 sa->strings = names;
5571 sa->n = sa->allocated = n_names;
5577 msave_common_uninit (struct msave_common *common)
5581 fh_unref (common->outfile);
5582 string_array_destroy (&common->variables);
5583 string_array_destroy (&common->fnames);
5584 string_array_destroy (&common->snames);
5589 compare_variables (const char *keyword,
5590 const struct string_array *new,
5591 const struct string_array *old)
5597 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5598 "on the first MSAVE within MATRIX."), keyword);
5601 else if (!string_array_equal_case (old, new))
5603 msg (SE, _("%s must specify the same variables each time within "
5604 "a given MATRIX."), keyword);
5611 static struct matrix_cmd *
5612 matrix_parse_msave (struct matrix_state *s)
5614 struct msave_common common = { .outfile = NULL };
5615 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5616 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5618 struct msave_command *msave = &cmd->msave;
5619 if (lex_token (s->lexer) == T_ID)
5620 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5621 msave->expr = matrix_parse_exp (s);
5625 while (lex_match (s->lexer, T_SLASH))
5627 if (lex_match_id (s->lexer, "TYPE"))
5629 lex_match (s->lexer, T_EQUALS);
5631 msave->rowtype = match_rowtype (s->lexer);
5632 if (!msave->rowtype)
5635 else if (lex_match_id (s->lexer, "OUTFILE"))
5637 lex_match (s->lexer, T_EQUALS);
5639 fh_unref (common.outfile);
5640 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5641 if (!common.outfile)
5644 else if (lex_match_id (s->lexer, "VARIABLES"))
5646 if (!parse_var_names (s->lexer, &common.variables))
5649 else if (lex_match_id (s->lexer, "FNAMES"))
5651 if (!parse_var_names (s->lexer, &common.fnames))
5654 else if (lex_match_id (s->lexer, "SNAMES"))
5656 if (!parse_var_names (s->lexer, &common.snames))
5659 else if (lex_match_id (s->lexer, "SPLIT"))
5661 lex_match (s->lexer, T_EQUALS);
5663 matrix_expr_destroy (msave->splits);
5664 msave->splits = matrix_parse_exp (s);
5668 else if (lex_match_id (s->lexer, "FACTOR"))
5670 lex_match (s->lexer, T_EQUALS);
5672 matrix_expr_destroy (msave->factors);
5673 msave->factors = matrix_parse_exp (s);
5674 if (!msave->factors)
5679 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5680 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5684 if (!msave->rowtype)
5686 lex_sbc_missing ("TYPE");
5689 common.has_splits = msave->splits || common.snames.n;
5690 common.has_factors = msave->factors || common.fnames.n;
5692 struct msave_common *c = s->common ? s->common : &common;
5693 if (c->fnames.n && !msave->factors)
5695 msg (SE, _("FNAMES requires FACTOR."));
5698 if (c->snames.n && !msave->splits)
5700 msg (SE, _("SNAMES requires SPLIT."));
5703 if (c->has_factors && !common.has_factors)
5705 msg (SE, _("%s is required because it was present on the first "
5706 "MSAVE in this MATRIX command."), "FACTOR");
5709 if (c->has_splits && !common.has_splits)
5711 msg (SE, _("%s is required because it was present on the first "
5712 "MSAVE in this MATRIX command."), "SPLIT");
5718 if (!common.outfile)
5720 lex_sbc_missing ("OUTFILE");
5723 s->common = xmemdup (&common, sizeof common);
5729 bool same = common.outfile == s->common->outfile;
5730 fh_unref (common.outfile);
5733 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5734 "within a single MATRIX command."));
5738 if (!compare_variables ("VARIABLES",
5739 &common.variables, &s->common->variables)
5740 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5741 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5743 msave_common_uninit (&common);
5745 msave->common = s->common;
5746 if (!msave->varname_)
5747 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5751 msave_common_uninit (&common);
5752 matrix_cmd_destroy (cmd);
5757 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5759 gsl_matrix *m = matrix_expr_evaluate (e);
5765 msg (SE, _("%s expression must evaluate to vector, "
5766 "not a %zu×%zu matrix."),
5767 name, m->size1, m->size2);
5768 gsl_matrix_free (m);
5772 return matrix_to_vector (m);
5776 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5778 for (size_t i = 0; i < vars->n; i++)
5779 if (!dict_create_var (d, vars->strings[i], 0))
5780 return vars->strings[i];
5784 static struct dictionary *
5785 msave_create_dict (const struct msave_common *common)
5787 struct dictionary *dict = dict_create (get_default_encoding ());
5789 const char *dup_split = msave_add_vars (dict, &common->snames);
5792 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5796 dict_create_var_assert (dict, "ROWTYPE_", 8);
5798 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5801 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5805 dict_create_var_assert (dict, "VARNAME_", 8);
5807 const char *dup_var = msave_add_vars (dict, &common->variables);
5810 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5822 matrix_cmd_execute_msave (struct msave_command *msave)
5824 struct msave_common *common = msave->common;
5825 gsl_matrix *m = NULL;
5826 gsl_vector *factors = NULL;
5827 gsl_vector *splits = NULL;
5829 m = matrix_expr_evaluate (msave->expr);
5833 if (!common->variables.n)
5834 for (size_t i = 0; i < m->size2; i++)
5835 string_array_append_nocopy (&common->variables,
5836 xasprintf ("COL%zu", i + 1));
5838 if (m->size2 != common->variables.n)
5840 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5841 m->size2, common->variables.n);
5847 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5851 if (!common->fnames.n)
5852 for (size_t i = 0; i < factors->size; i++)
5853 string_array_append_nocopy (&common->fnames,
5854 xasprintf ("FAC%zu", i + 1));
5856 if (factors->size != common->fnames.n)
5858 msg (SE, _("There are %zu factor variables, "
5859 "but %zu split values were supplied."),
5860 common->fnames.n, factors->size);
5867 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5871 if (!common->fnames.n)
5872 for (size_t i = 0; i < splits->size; i++)
5873 string_array_append_nocopy (&common->fnames,
5874 xasprintf ("SPL%zu", i + 1));
5876 if (splits->size != common->fnames.n)
5878 msg (SE, _("There are %zu split variables, "
5879 "but %zu split values were supplied."),
5880 common->fnames.n, splits->size);
5885 if (!common->writer)
5887 struct dictionary *dict = msave_create_dict (common);
5891 common->writer = any_writer_open (common->outfile, dict);
5892 if (!common->writer)
5898 common->dict = dict;
5901 for (size_t y = 0; y < m->size1; y++)
5903 struct ccase *c = case_create (dict_get_proto (common->dict));
5906 /* Split variables */
5908 for (size_t i = 0; i < splits->size; i++)
5909 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5912 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5913 msave->rowtype, ' ');
5917 for (size_t i = 0; i < factors->size; i++)
5918 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5921 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5922 msave->varname_, ' ');
5924 /* Continuous variables. */
5925 for (size_t x = 0; x < m->size2; x++)
5926 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5927 casewriter_write (common->writer, c);
5931 gsl_matrix_free (m);
5932 gsl_vector_free (factors);
5933 gsl_vector_free (splits);
5936 static struct matrix_cmd *
5937 matrix_parse_mget (struct matrix_state *s)
5939 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5940 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5942 struct mget_command *mget = &cmd->mget;
5946 if (lex_match_id (s->lexer, "FILE"))
5948 lex_match (s->lexer, T_EQUALS);
5950 fh_unref (mget->file);
5951 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5955 else if (lex_match_id (s->lexer, "TYPE"))
5957 lex_match (s->lexer, T_EQUALS);
5958 while (lex_token (s->lexer) != T_SLASH
5959 && lex_token (s->lexer) != T_ENDCMD)
5961 const char *rowtype = match_rowtype (s->lexer);
5965 stringi_set_insert (&mget->rowtypes, rowtype);
5970 lex_error_expecting (s->lexer, "FILE", "TYPE");
5973 if (lex_token (s->lexer) == T_ENDCMD)
5976 if (!lex_force_match (s->lexer, T_SLASH))
5982 matrix_cmd_destroy (cmd);
5986 static const struct variable *
5987 get_a8_var (const struct dictionary *d, const char *name)
5989 const struct variable *v = dict_lookup_var (d, name);
5992 msg (SE, _("Matrix data file lacks %s variable."), name);
5995 if (var_get_width (v) != 8)
5997 msg (SE, _("%s variable in matrix data file must be "
5998 "8-byte string, but it has width %d."),
5999 name, var_get_width (v));
6006 is_rowtype (const union value *v, const char *rowtype)
6008 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6009 ss_rtrim (&vs, ss_cstr (" "));
6010 return ss_equals_case (vs, ss_cstr (rowtype));
6014 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6015 const struct dictionary *d,
6016 const struct variable *rowtype_var,
6017 struct matrix_state *s, size_t si, size_t fi,
6018 size_t cs, size_t cn)
6023 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6024 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6025 : is_rowtype (rowtype_, "CORR") ? "CR"
6026 : is_rowtype (rowtype_, "MEAN") ? "MN"
6027 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6028 : is_rowtype (rowtype_, "N") ? "NC"
6031 struct string name = DS_EMPTY_INITIALIZER;
6032 ds_put_cstr (&name, name_prefix);
6034 ds_put_format (&name, "F%zu", fi);
6036 ds_put_format (&name, "S%zu", si);
6038 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6040 mv = matrix_var_create (s, ds_ss (&name));
6043 msg (SW, _("Matrix data file contains variable with existing name %s."),
6048 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6049 size_t n_missing = 0;
6050 for (size_t y = 0; y < n_rows; y++)
6052 for (size_t x = 0; x < cn; x++)
6054 struct variable *var = dict_get_var (d, cs + x);
6055 double value = case_num (rows[y], var);
6056 if (var_is_num_missing (var, value, MV_ANY))
6061 gsl_matrix_set (m, y, x, value);
6066 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6067 "value, which was treated as zero.",
6068 "Matrix data file variable %s contains %zu missing "
6069 "values, which were treated as zero.", n_missing),
6070 ds_cstr (&name), n_missing);
6075 for (size_t y = 0; y < n_rows; y++)
6076 case_unref (rows[y]);
6080 var_changed (const struct ccase *ca, const struct ccase *cb,
6081 const struct variable *var)
6083 return !value_equal (case_data (ca, var), case_data (cb, var),
6084 var_get_width (var));
6088 vars_changed (const struct ccase *ca, const struct ccase *cb,
6089 const struct dictionary *d,
6090 size_t first_var, size_t n_vars)
6092 for (size_t i = 0; i < n_vars; i++)
6094 const struct variable *v = dict_get_var (d, first_var + i);
6095 if (var_changed (ca, cb, v))
6102 matrix_cmd_execute_mget (struct mget_command *mget)
6104 struct dictionary *d;
6105 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6110 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6111 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6112 if (!rowtype_ || !varname_)
6115 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6117 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6120 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6122 msg (SE, _("Matrix data file contains no continuous variables."));
6126 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6128 const struct variable *v = dict_get_var (d, i);
6129 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6132 _("Matrix data file contains unexpected string variable %s."),
6138 /* SPLIT variables. */
6140 size_t sn = var_get_dict_index (rowtype_);
6141 struct ccase *sc = NULL;
6144 /* FACTOR variables. */
6145 size_t fs = var_get_dict_index (rowtype_) + 1;
6146 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6147 struct ccase *fc = NULL;
6150 /* Continuous variables. */
6151 size_t cs = var_get_dict_index (varname_) + 1;
6152 size_t cn = dict_get_var_cnt (d) - cs;
6153 struct ccase *cc = NULL;
6156 struct ccase **rows = NULL;
6157 size_t allocated_rows = 0;
6161 while ((c = casereader_read (r)) != NULL)
6163 bool sd = vars_changed (sc, c, d, ss, sn);
6164 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6165 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6180 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6181 mget->state, si, fi, cs, cn);
6187 if (n_rows >= allocated_rows)
6188 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6191 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6192 mget->state, si, fi, cs, cn);
6196 casereader_destroy (r);
6200 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6202 if (!lex_force_id (s->lexer))
6205 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6207 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6212 static struct matrix_cmd *
6213 matrix_parse_eigen (struct matrix_state *s)
6215 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6216 *cmd = (struct matrix_cmd) {
6218 .eigen = { .expr = NULL }
6221 struct eigen_command *eigen = &cmd->eigen;
6222 if (!lex_force_match (s->lexer, T_LPAREN))
6224 eigen->expr = matrix_parse_expr (s);
6226 || !lex_force_match (s->lexer, T_COMMA)
6227 || !matrix_parse_dst_var (s, &eigen->evec)
6228 || !lex_force_match (s->lexer, T_COMMA)
6229 || !matrix_parse_dst_var (s, &eigen->eval)
6230 || !lex_force_match (s->lexer, T_RPAREN))
6236 matrix_cmd_destroy (cmd);
6241 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6243 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6246 if (!is_symmetric (A))
6248 msg (SE, _("Argument of EIGEN must be symmetric."));
6249 gsl_matrix_free (A);
6253 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6254 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6255 gsl_vector v_eval = to_vector (eval);
6256 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6257 gsl_eigen_symmv (A, &v_eval, evec, w);
6258 gsl_eigen_symmv_free (w);
6260 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6262 gsl_matrix_free (A);
6264 gsl_matrix_free (eigen->eval->value);
6265 eigen->eval->value = eval;
6267 gsl_matrix_free (eigen->evec->value);
6268 eigen->evec->value = evec;
6271 static struct matrix_cmd *
6272 matrix_parse_setdiag (struct matrix_state *s)
6274 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6275 *cmd = (struct matrix_cmd) {
6276 .type = MCMD_SETDIAG,
6277 .setdiag = { .dst = NULL }
6280 struct setdiag_command *setdiag = &cmd->setdiag;
6281 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6284 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6287 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6292 if (!lex_force_match (s->lexer, T_COMMA))
6295 setdiag->expr = matrix_parse_expr (s);
6299 if (!lex_force_match (s->lexer, T_RPAREN))
6305 matrix_cmd_destroy (cmd);
6310 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6312 gsl_matrix *dst = setdiag->dst->value;
6315 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6316 setdiag->dst->name);
6320 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6324 size_t n = MIN (dst->size1, dst->size2);
6325 if (is_scalar (src))
6327 double d = to_scalar (src);
6328 for (size_t i = 0; i < n; i++)
6329 gsl_matrix_set (dst, i, i, d);
6331 else if (is_vector (src))
6333 gsl_vector v = to_vector (src);
6334 for (size_t i = 0; i < n && i < v.size; i++)
6335 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6338 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6339 "not a %zu×%zu matrix."),
6340 src->size1, src->size2);
6341 gsl_matrix_free (src);
6344 static struct matrix_cmd *
6345 matrix_parse_svd (struct matrix_state *s)
6347 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6348 *cmd = (struct matrix_cmd) {
6350 .svd = { .expr = NULL }
6353 struct svd_command *svd = &cmd->svd;
6354 if (!lex_force_match (s->lexer, T_LPAREN))
6356 svd->expr = matrix_parse_expr (s);
6358 || !lex_force_match (s->lexer, T_COMMA)
6359 || !matrix_parse_dst_var (s, &svd->u)
6360 || !lex_force_match (s->lexer, T_COMMA)
6361 || !matrix_parse_dst_var (s, &svd->s)
6362 || !lex_force_match (s->lexer, T_COMMA)
6363 || !matrix_parse_dst_var (s, &svd->v)
6364 || !lex_force_match (s->lexer, T_RPAREN))
6370 matrix_cmd_destroy (cmd);
6375 matrix_cmd_execute_svd (struct svd_command *svd)
6377 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6381 if (m->size1 >= m->size2)
6384 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6385 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6386 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6387 gsl_vector *work = gsl_vector_alloc (A->size2);
6388 gsl_linalg_SV_decomp (A, V, &Sv, work);
6389 gsl_vector_free (work);
6391 matrix_var_set (svd->u, A);
6392 matrix_var_set (svd->s, S);
6393 matrix_var_set (svd->v, V);
6397 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6398 gsl_matrix_transpose_memcpy (At, m);
6399 gsl_matrix_free (m);
6401 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6402 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6403 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6404 gsl_vector *work = gsl_vector_alloc (At->size2);
6405 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6406 gsl_vector_free (work);
6408 matrix_var_set (svd->v, At);
6409 matrix_var_set (svd->s, St);
6410 matrix_var_set (svd->u, Vt);
6415 matrix_cmd_execute (struct matrix_cmd *cmd)
6420 matrix_cmd_execute_compute (&cmd->compute);
6424 matrix_cmd_execute_print (&cmd->print);
6428 return matrix_cmd_execute_do_if (&cmd->do_if);
6431 matrix_cmd_execute_loop (&cmd->loop);
6438 matrix_cmd_execute_display (&cmd->display);
6442 matrix_cmd_execute_release (&cmd->release);
6446 matrix_cmd_execute_save (&cmd->save);
6450 matrix_cmd_execute_read (&cmd->read);
6454 matrix_cmd_execute_write (&cmd->write);
6458 matrix_cmd_execute_get (&cmd->get);
6462 matrix_cmd_execute_msave (&cmd->msave);
6466 matrix_cmd_execute_mget (&cmd->mget);
6470 matrix_cmd_execute_eigen (&cmd->eigen);
6474 matrix_cmd_execute_setdiag (&cmd->setdiag);
6478 matrix_cmd_execute_svd (&cmd->svd);
6486 matrix_cmds_uninit (struct matrix_cmds *cmds)
6488 for (size_t i = 0; i < cmds->n; i++)
6489 matrix_cmd_destroy (cmds->commands[i]);
6490 free (cmds->commands);
6494 matrix_cmd_destroy (struct matrix_cmd *cmd)
6502 matrix_lvalue_destroy (cmd->compute.lvalue);
6503 matrix_expr_destroy (cmd->compute.rvalue);
6507 matrix_expr_destroy (cmd->print.expression);
6508 free (cmd->print.title);
6509 print_labels_destroy (cmd->print.rlabels);
6510 print_labels_destroy (cmd->print.clabels);
6514 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6516 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6517 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6519 free (cmd->do_if.clauses);
6523 matrix_expr_destroy (cmd->loop.start);
6524 matrix_expr_destroy (cmd->loop.end);
6525 matrix_expr_destroy (cmd->loop.increment);
6526 matrix_expr_destroy (cmd->loop.top_condition);
6527 matrix_expr_destroy (cmd->loop.bottom_condition);
6528 matrix_cmds_uninit (&cmd->loop.commands);
6538 free (cmd->release.vars);
6542 matrix_expr_destroy (cmd->save.expression);
6543 fh_unref (cmd->save.outfile);
6544 string_array_destroy (cmd->save.variables);
6545 matrix_expr_destroy (cmd->save.names);
6546 stringi_set_destroy (&cmd->save.strings);
6550 matrix_lvalue_destroy (cmd->read.dst);
6551 matrix_expr_destroy (cmd->read.size);
6555 matrix_expr_destroy (cmd->write.expression);
6559 matrix_lvalue_destroy (cmd->get.dst);
6560 fh_unref (cmd->get.file);
6561 free (cmd->get.encoding);
6562 string_array_destroy (&cmd->get.variables);
6566 free (cmd->msave.varname_);
6567 matrix_expr_destroy (cmd->msave.expr);
6568 matrix_expr_destroy (cmd->msave.factors);
6569 matrix_expr_destroy (cmd->msave.splits);
6573 fh_unref (cmd->mget.file);
6574 stringi_set_destroy (&cmd->mget.rowtypes);
6578 matrix_expr_destroy (cmd->eigen.expr);
6582 matrix_expr_destroy (cmd->setdiag.expr);
6586 matrix_expr_destroy (cmd->svd.expr);
6592 struct matrix_command_name
6595 struct matrix_cmd *(*parse) (struct matrix_state *);
6598 static const struct matrix_command_name *
6599 matrix_parse_command_name (struct lexer *lexer)
6601 static const struct matrix_command_name commands[] = {
6602 { "COMPUTE", matrix_parse_compute },
6603 { "CALL EIGEN", matrix_parse_eigen },
6604 { "CALL SETDIAG", matrix_parse_setdiag },
6605 { "CALL SVD", matrix_parse_svd },
6606 { "PRINT", matrix_parse_print },
6607 { "DO IF", matrix_parse_do_if },
6608 { "LOOP", matrix_parse_loop },
6609 { "BREAK", matrix_parse_break },
6610 { "READ", matrix_parse_read },
6611 { "WRITE", matrix_parse_write },
6612 { "GET", matrix_parse_get },
6613 { "SAVE", matrix_parse_save },
6614 { "MGET", matrix_parse_mget },
6615 { "MSAVE", matrix_parse_msave },
6616 { "DISPLAY", matrix_parse_display },
6617 { "RELEASE", matrix_parse_release },
6619 static size_t n = sizeof commands / sizeof *commands;
6621 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6622 if (lex_match_phrase (lexer, c->name))
6627 static struct matrix_cmd *
6628 matrix_parse_command (struct matrix_state *s)
6630 size_t nesting_level = SIZE_MAX;
6632 struct matrix_cmd *c = NULL;
6633 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6635 lex_error (s->lexer, _("Unknown matrix command."));
6636 else if (!cmd->parse)
6637 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6641 nesting_level = output_open_group (
6642 group_item_create_nocopy (utf8_to_title (cmd->name),
6643 utf8_to_title (cmd->name)));
6648 lex_end_of_command (s->lexer);
6649 lex_discard_rest_of_command (s->lexer);
6650 if (nesting_level != SIZE_MAX)
6651 output_close_groups (nesting_level);
6657 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6659 if (!lex_force_match (lexer, T_ENDCMD))
6662 struct matrix_state state = {
6663 .session = dataset_session (ds),
6665 .vars = HMAP_INITIALIZER (state.vars),
6670 while (lex_match (lexer, T_ENDCMD))
6672 if (lex_token (lexer) == T_STOP)
6674 msg (SE, _("Unexpected end of input expecting matrix command."));
6678 if (lex_match_phrase (lexer, "END MATRIX"))
6681 struct matrix_cmd *c = matrix_parse_command (&state);
6684 matrix_cmd_execute (c);
6685 matrix_cmd_destroy (c);
6689 struct matrix_var *var, *next;
6690 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6693 gsl_matrix_free (var->value);
6694 hmap_delete (&state.vars, &var->hmap_node);
6697 hmap_destroy (&state.vars);
6698 fh_unref (state.prev_save_outfile);
6701 dict_unref (state.common->dict);
6702 casewriter_destroy (state.common->writer);
6703 free (state.common);
6705 fh_unref (state.prev_read_file);
6706 for (size_t i = 0; i < state.n_read_files; i++)
6707 read_file_destroy (state.read_files[i]);
6708 free (state.read_files);
6709 fh_unref (state.prev_write_file);
6710 for (size_t i = 0; i < state.n_write_files; i++)
6711 write_file_destroy (state.write_files[i]);
6712 free (state.write_files);