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 read->size = matrix_parse_exp (s);
4637 else if (lex_match_id (s->lexer, "MODE"))
4639 lex_match (s->lexer, T_EQUALS);
4640 if (lex_match_id (s->lexer, "RECTANGULAR"))
4641 read->symmetric = false;
4642 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4643 read->symmetric = true;
4646 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4650 else if (lex_match_id (s->lexer, "REREAD"))
4651 read->reread = true;
4652 else if (lex_match_id (s->lexer, "FORMAT"))
4656 lex_sbc_only_once ("FORMAT");
4661 lex_match (s->lexer, T_EQUALS);
4663 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4666 const char *p = lex_tokcstr (s->lexer);
4667 if (c_isdigit (p[0]))
4669 repetitions = atoi (p);
4670 p += strspn (p, "0123456789");
4671 if (!fmt_from_name (p, &read->format))
4673 lex_error (s->lexer, _("Unknown format %s."), p);
4678 else if (!fmt_from_name (p, &read->format))
4680 struct fmt_spec format;
4681 if (!parse_format_specifier (s->lexer, &format))
4683 read->format = format.type;
4690 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4691 "REREAD", "FORMAT");
4698 lex_sbc_missing ("FIELD");
4702 if (!read->dst->n_indexes && !read->size)
4704 msg (SE, _("SIZE is required for reading data into a full matrix "
4705 "(as opposed to a submatrix)."));
4711 if (s->prev_read_file)
4712 fh = fh_ref (s->prev_read_file);
4715 lex_sbc_missing ("FILE");
4719 fh_unref (s->prev_read_file);
4720 s->prev_read_file = fh_ref (fh);
4722 read->rf = read_file_create (s, fh);
4726 free (read->rf->encoding);
4727 read->rf->encoding = encoding;
4731 /* Field width may be specified in multiple ways:
4734 2. The format on FORMAT.
4735 3. The repetition factor on FORMAT.
4737 (2) and (3) are mutually exclusive.
4739 If more than one of these is present, they must agree. If none of them is
4740 present, then free-field format is used.
4742 if (repetitions > record_width)
4744 msg (SE, _("%d repetitions cannot fit in record width %d."),
4745 repetitions, record_width);
4748 int w = (repetitions ? record_width / repetitions
4754 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4755 "which implies field width %d, "
4756 "but BY specifies field width %d."),
4757 repetitions, record_width, w, by);
4759 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4768 matrix_cmd_destroy (cmd);
4774 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4775 gsl_matrix *m, struct substring p, size_t y, size_t x)
4777 const char *input_encoding = dfm_reader_get_encoding (reader);
4779 char *error = data_in (p, input_encoding, read->format,
4780 settings_get_fmt_settings (), &v, 0, NULL);
4781 /* XXX report error if value is missing */
4783 msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4786 gsl_matrix_set (m, y, x, v.f);
4787 if (read->symmetric && x != y)
4788 gsl_matrix_set (m, x, y, v.f);
4793 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4794 struct substring *line)
4796 if (dfm_eof (reader))
4798 msg (SE, _("Unexpected end of file reading matrix data."));
4801 dfm_expand_tabs (reader);
4802 *line = ss_substr (dfm_get_record (reader),
4803 read->c1 - 1, read->c2 - read->c1);
4808 matrix_read (struct read_command *read, struct dfm_reader *reader,
4811 for (size_t y = 0; y < m->size1; y++)
4813 size_t nx = read->symmetric ? y + 1 : m->size2;
4815 struct substring line = ss_empty ();
4816 for (size_t x = 0; x < nx; x++)
4823 ss_ltrim (&line, ss_cstr (" ,"));
4824 if (!ss_is_empty (line))
4826 if (!matrix_read_line (read, reader, &line))
4828 dfm_forward_record (reader);
4831 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4835 if (!matrix_read_line (read, reader, &line))
4837 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4838 int f = x % fields_per_line;
4839 if (f == fields_per_line - 1)
4840 dfm_forward_record (reader);
4842 p = ss_substr (line, read->w * f, read->w);
4845 matrix_read_set_field (read, reader, m, p, y, x);
4849 dfm_forward_record (reader);
4852 ss_ltrim (&line, ss_cstr (" ,"));
4853 if (!ss_is_empty (line))
4854 msg (SW, _("Trailing garbage on line \"%.*s\""),
4855 (int) line.length, line.string);
4861 matrix_cmd_execute_read (struct read_command *read)
4863 struct index_vector iv0, iv1;
4864 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4867 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4870 gsl_matrix *m = matrix_expr_evaluate (read->size);
4876 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4877 "not a %zu×%zu matrix."), m->size1, m->size2);
4878 gsl_matrix_free (m);
4884 gsl_vector v = to_vector (m);
4888 d[0] = gsl_vector_get (&v, 0);
4891 else if (v.size == 2)
4893 d[0] = gsl_vector_get (&v, 0);
4894 d[1] = gsl_vector_get (&v, 1);
4898 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4899 "not a %zu×%zu matrix."),
4900 m->size1, m->size2),
4901 gsl_matrix_free (m);
4906 gsl_matrix_free (m);
4908 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4910 msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4920 if (read->dst->n_indexes)
4922 size_t submatrix_size[2];
4923 if (read->dst->n_indexes == 2)
4925 submatrix_size[0] = iv0.n;
4926 submatrix_size[1] = iv1.n;
4928 else if (read->dst->var->value->size1 == 1)
4930 submatrix_size[0] = 1;
4931 submatrix_size[1] = iv0.n;
4935 submatrix_size[0] = iv0.n;
4936 submatrix_size[1] = 1;
4941 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4943 msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4946 submatrix_size[0], submatrix_size[1]);
4954 size[0] = submatrix_size[0];
4955 size[1] = submatrix_size[1];
4959 struct dfm_reader *reader = read_file_open (read->rf);
4961 dfm_reread_record (reader, 1);
4963 if (read->symmetric && size[0] != size[1])
4965 msg (SE, _("Cannot read non-square %zu×%zu matrix "
4966 "using READ with MODE=SYMMETRIC."),
4972 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4973 matrix_read (read, reader, tmp);
4974 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4977 static struct matrix_cmd *
4978 matrix_parse_write (struct matrix_state *s)
4980 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4981 *cmd = (struct matrix_cmd) {
4983 .write = { .format = FMT_F },
4986 struct file_handle *fh = NULL;
4987 char *encoding = NULL;
4988 struct write_command *write = &cmd->write;
4989 write->expression = matrix_parse_exp (s);
4990 if (!write->expression)
4994 int repetitions = 0;
4995 int record_width = 0;
4996 bool seen_format = false;
4997 while (lex_match (s->lexer, T_SLASH))
4999 if (lex_match_id (s->lexer, "OUTFILE"))
5001 lex_match (s->lexer, T_EQUALS);
5004 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5008 else if (lex_match_id (s->lexer, "ENCODING"))
5010 lex_match (s->lexer, T_EQUALS);
5011 if (!lex_force_string (s->lexer))
5015 encoding = ss_xstrdup (lex_tokss (s->lexer));
5019 else if (lex_match_id (s->lexer, "FIELD"))
5021 lex_match (s->lexer, T_EQUALS);
5023 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5025 write->c1 = lex_integer (s->lexer);
5027 if (!lex_force_match (s->lexer, T_TO)
5028 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5030 write->c2 = lex_integer (s->lexer) + 1;
5033 record_width = write->c2 - write->c1;
5034 if (lex_match (s->lexer, T_BY))
5036 if (!lex_force_int_range (s->lexer, "BY", 1,
5037 write->c2 - write->c1))
5039 by = lex_integer (s->lexer);
5042 if (record_width % by)
5044 msg (SE, _("BY %d does not evenly divide record width %d."),
5052 else if (lex_match_id (s->lexer, "MODE"))
5054 lex_match (s->lexer, T_EQUALS);
5055 if (lex_match_id (s->lexer, "RECTANGULAR"))
5056 write->triangular = false;
5057 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5058 write->triangular = true;
5061 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5065 else if (lex_match_id (s->lexer, "HOLD"))
5067 else if (lex_match_id (s->lexer, "FORMAT"))
5071 lex_sbc_only_once ("FORMAT");
5076 lex_match (s->lexer, T_EQUALS);
5078 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5081 const char *p = lex_tokcstr (s->lexer);
5082 if (c_isdigit (p[0]))
5084 repetitions = atoi (p);
5085 p += strspn (p, "0123456789");
5086 if (!fmt_from_name (p, &write->format))
5088 lex_error (s->lexer, _("Unknown format %s."), p);
5093 else if (!fmt_from_name (p, &write->format))
5095 struct fmt_spec format;
5096 if (!parse_format_specifier (s->lexer, &format))
5098 write->format = format.type;
5099 write->w = format.w;
5100 write->d = format.d;
5105 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5113 lex_sbc_missing ("FIELD");
5119 if (s->prev_write_file)
5120 fh = fh_ref (s->prev_write_file);
5123 lex_sbc_missing ("OUTFILE");
5127 fh_unref (s->prev_write_file);
5128 s->prev_write_file = fh_ref (fh);
5130 write->wf = write_file_create (s, fh);
5134 free (write->wf->encoding);
5135 write->wf->encoding = encoding;
5139 /* Field width may be specified in multiple ways:
5142 2. The format on FORMAT.
5143 3. The repetition factor on FORMAT.
5145 (2) and (3) are mutually exclusive.
5147 If more than one of these is present, they must agree. If none of them is
5148 present, then free-field format is used.
5150 if (repetitions > record_width)
5152 msg (SE, _("%d repetitions cannot fit in record width %d."),
5153 repetitions, record_width);
5156 int w = (repetitions ? record_width / repetitions
5157 : write->w ? write->w
5162 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5163 "which implies field width %d, "
5164 "but BY specifies field width %d."),
5165 repetitions, record_width, w, by);
5167 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5176 matrix_cmd_destroy (cmd);
5181 matrix_cmd_execute_write (struct write_command *write)
5183 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5187 if (write->triangular && m->size1 != m->size2)
5189 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5190 "the matrix to be written has dimensions %zu×%zu."),
5191 m->size1, m->size2);
5192 gsl_matrix_free (m);
5196 struct dfm_writer *writer = write_file_open (write->wf);
5197 if (!writer || !m->size1)
5199 gsl_matrix_free (m);
5203 const struct fmt_settings *settings = settings_get_fmt_settings ();
5204 struct fmt_spec format = {
5205 .type = write->format,
5206 .w = write->w ? write->w : 40,
5209 struct u8_line *line = write->wf->held;
5210 for (size_t y = 0; y < m->size1; y++)
5214 line = xmalloc (sizeof *line);
5215 u8_line_init (line);
5217 size_t nx = write->triangular ? y + 1 : m->size2;
5219 for (size_t x = 0; x < nx; x++)
5221 /* XXX string values */
5222 union value v = { .f = gsl_matrix_get (m, y, x) };
5224 ? data_out (&v, NULL, &format, settings)
5225 : data_out_stretchy (&v, NULL, &format, settings, NULL));
5226 size_t len = strlen (s);
5227 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5228 if (width + x0 > write->c2)
5230 dfm_put_record_utf8 (writer, line->s.ss.string,
5232 u8_line_clear (line);
5235 u8_line_put (line, x0, x0 + width, s, len);
5238 x0 += write->w ? write->w : width + 1;
5241 if (y + 1 >= m->size1 && write->hold)
5243 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5244 u8_line_clear (line);
5248 u8_line_destroy (line);
5251 write->wf->held = line;
5252 dfm_close_writer (writer);
5254 gsl_matrix_free (m);
5257 static struct matrix_cmd *
5258 matrix_parse_get (struct matrix_state *s)
5260 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5261 *cmd = (struct matrix_cmd) {
5264 .user = { .treatment = MGET_ERROR },
5265 .system = { .treatment = MGET_ERROR },
5269 struct get_command *get = &cmd->get;
5270 get->dst = matrix_lvalue_parse (s);
5274 while (lex_match (s->lexer, T_SLASH))
5276 if (lex_match_id (s->lexer, "FILE"))
5278 if (get->variables.n)
5280 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5283 lex_match (s->lexer, T_EQUALS);
5285 fh_unref (get->file);
5286 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5290 else if (lex_match_id (s->lexer, "ENCODING"))
5292 if (get->variables.n)
5294 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5297 lex_match (s->lexer, T_EQUALS);
5298 if (!lex_force_string (s->lexer))
5301 free (get->encoding);
5302 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5306 else if (lex_match_id (s->lexer, "VARIABLES"))
5308 lex_match (s->lexer, T_EQUALS);
5310 struct dictionary *dict = NULL;
5313 dict = dataset_dict (s->dataset);
5314 if (dict_get_var_cnt (dict) == 0)
5316 lex_error (s->lexer, _("GET cannot read empty active file."));
5322 struct casereader *reader = any_reader_open_and_decode (
5323 get->file, get->encoding, &dict, NULL);
5326 casereader_destroy (reader);
5329 struct variable **vars;
5331 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5332 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5339 string_array_clear (&get->variables);
5340 for (size_t i = 0; i < n_vars; i++)
5341 string_array_append (&get->variables, var_get_name (vars[i]));
5345 else if (lex_match_id (s->lexer, "NAMES"))
5347 lex_match (s->lexer, T_EQUALS);
5348 if (!lex_force_id (s->lexer))
5351 struct substring name = lex_tokss (s->lexer);
5352 get->names = matrix_var_lookup (s, name);
5354 get->names = matrix_var_create (s, name);
5357 else if (lex_match_id (s->lexer, "MISSING"))
5359 lex_match (s->lexer, T_EQUALS);
5360 if (lex_match_id (s->lexer, "ACCEPT"))
5361 get->user.treatment = MGET_ACCEPT;
5362 else if (lex_match_id (s->lexer, "OMIT"))
5363 get->user.treatment = MGET_OMIT;
5364 else if (lex_is_number (s->lexer))
5366 get->user.treatment = MGET_RECODE;
5367 get->user.substitute = lex_number (s->lexer);
5372 lex_error (s->lexer, NULL);
5376 else if (lex_match_id (s->lexer, "SYSMIS"))
5378 lex_match (s->lexer, T_EQUALS);
5379 if (lex_match_id (s->lexer, "OMIT"))
5380 get->user.treatment = MGET_OMIT;
5381 else if (lex_is_number (s->lexer))
5383 get->user.treatment = MGET_RECODE;
5384 get->user.substitute = lex_number (s->lexer);
5389 lex_error (s->lexer, NULL);
5395 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5396 "MISSING", "SYSMIS");
5403 matrix_cmd_destroy (cmd);
5408 matrix_cmd_execute_get (struct get_command *get)
5410 assert (get->file); /* XXX */
5412 struct dictionary *dict;
5413 struct casereader *reader = any_reader_open_and_decode (
5414 get->file, get->encoding, &dict, NULL);
5418 const struct variable **vars = xnmalloc (
5419 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5423 if (get->variables.n)
5425 for (size_t i = 0; i < get->variables.n; i++)
5427 const char *name = get->variables.strings[i];
5428 const struct variable *var = dict_lookup_var (dict, name);
5431 msg (SE, _("GET: Data file does not contain variable %s."),
5437 if (!var_is_numeric (var))
5439 msg (SE, _("GET: Variable %s is not numeric."), name);
5444 vars[n_vars++] = var;
5449 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5451 const struct variable *var = dict_get_var (dict, i);
5452 if (!var_is_numeric (var))
5454 msg (SE, _("GET: Variable %s is not numeric."),
5455 var_get_name (var));
5460 vars[n_vars++] = var;
5465 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5466 long long int casenum = 1;
5468 for (struct ccase *c = casereader_read (reader); c;
5469 c = casereader_read (reader), casenum++)
5471 if (n_rows >= m->size1)
5473 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5474 for (size_t y = 0; y < n_rows; y++)
5475 for (size_t x = 0; x < n_vars; x++)
5476 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5477 gsl_matrix_free (m);
5482 for (size_t x = 0; x < n_vars; x++)
5484 const struct variable *var = vars[x];
5485 double d = case_num (c, var);
5488 if (get->system.treatment == MGET_RECODE)
5489 d = get->system.substitute;
5490 else if (get->system.treatment == MGET_OMIT)
5494 msg (SE, _("GET: Variable %s in case %lld "
5495 "is system-missing."),
5496 var_get_name (var), casenum);
5500 else if (var_is_num_missing (var, d, MV_USER))
5502 if (get->user.treatment == MGET_RECODE)
5503 d = get->user.substitute;
5504 else if (get->user.treatment == MGET_OMIT)
5506 else if (get->user.treatment != MGET_ACCEPT)
5508 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5510 var_get_name (var), casenum, d);
5514 gsl_matrix_set (m, n_rows, x, d);
5522 casereader_destroy (reader);
5526 matrix_lvalue_evaluate_and_assign (get->dst, m);
5529 gsl_matrix_free (m);
5535 match_rowtype (struct lexer *lexer)
5537 static const char *rowtypes[] = {
5538 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5540 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5542 for (size_t i = 0; i < n_rowtypes; i++)
5543 if (lex_match_id (lexer, rowtypes[i]))
5546 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5551 parse_var_names (struct lexer *lexer, struct string_array *sa)
5553 lex_match (lexer, T_EQUALS);
5555 string_array_clear (sa);
5557 struct dictionary *dict = dict_create (get_default_encoding ());
5558 dict_create_var_assert (dict, "ROWTYPE_", 8);
5559 dict_create_var_assert (dict, "VARNAME_", 8);
5562 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5563 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5568 string_array_clear (sa);
5569 sa->strings = names;
5570 sa->n = sa->allocated = n_names;
5576 msave_common_uninit (struct msave_common *common)
5580 fh_unref (common->outfile);
5581 string_array_destroy (&common->variables);
5582 string_array_destroy (&common->fnames);
5583 string_array_destroy (&common->snames);
5588 compare_variables (const char *keyword,
5589 const struct string_array *new,
5590 const struct string_array *old)
5596 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5597 "on the first MSAVE within MATRIX."), keyword);
5600 else if (!string_array_equal_case (old, new))
5602 msg (SE, _("%s must specify the same variables each time within "
5603 "a given MATRIX."), keyword);
5610 static struct matrix_cmd *
5611 matrix_parse_msave (struct matrix_state *s)
5613 struct msave_common common = { .outfile = NULL };
5614 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5615 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5617 struct msave_command *msave = &cmd->msave;
5618 if (lex_token (s->lexer) == T_ID)
5619 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5620 msave->expr = matrix_parse_exp (s);
5624 while (lex_match (s->lexer, T_SLASH))
5626 if (lex_match_id (s->lexer, "TYPE"))
5628 lex_match (s->lexer, T_EQUALS);
5630 msave->rowtype = match_rowtype (s->lexer);
5631 if (!msave->rowtype)
5634 else if (lex_match_id (s->lexer, "OUTFILE"))
5636 lex_match (s->lexer, T_EQUALS);
5638 fh_unref (common.outfile);
5639 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5640 if (!common.outfile)
5643 else if (lex_match_id (s->lexer, "VARIABLES"))
5645 if (!parse_var_names (s->lexer, &common.variables))
5648 else if (lex_match_id (s->lexer, "FNAMES"))
5650 if (!parse_var_names (s->lexer, &common.fnames))
5653 else if (lex_match_id (s->lexer, "SNAMES"))
5655 if (!parse_var_names (s->lexer, &common.snames))
5658 else if (lex_match_id (s->lexer, "SPLIT"))
5660 lex_match (s->lexer, T_EQUALS);
5662 matrix_expr_destroy (msave->splits);
5663 msave->splits = matrix_parse_exp (s);
5667 else if (lex_match_id (s->lexer, "FACTOR"))
5669 lex_match (s->lexer, T_EQUALS);
5671 matrix_expr_destroy (msave->factors);
5672 msave->factors = matrix_parse_exp (s);
5673 if (!msave->factors)
5678 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5679 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5683 if (!msave->rowtype)
5685 lex_sbc_missing ("TYPE");
5688 common.has_splits = msave->splits || common.snames.n;
5689 common.has_factors = msave->factors || common.fnames.n;
5691 struct msave_common *c = s->common ? s->common : &common;
5692 if (c->fnames.n && !msave->factors)
5694 msg (SE, _("FNAMES requires FACTOR."));
5697 if (c->snames.n && !msave->splits)
5699 msg (SE, _("SNAMES requires SPLIT."));
5702 if (c->has_factors && !common.has_factors)
5704 msg (SE, _("%s is required because it was present on the first "
5705 "MSAVE in this MATRIX command."), "FACTOR");
5708 if (c->has_splits && !common.has_splits)
5710 msg (SE, _("%s is required because it was present on the first "
5711 "MSAVE in this MATRIX command."), "SPLIT");
5717 if (!common.outfile)
5719 lex_sbc_missing ("OUTFILE");
5722 s->common = xmemdup (&common, sizeof common);
5728 bool same = common.outfile == s->common->outfile;
5729 fh_unref (common.outfile);
5732 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5733 "within a single MATRIX command."));
5737 if (!compare_variables ("VARIABLES",
5738 &common.variables, &s->common->variables)
5739 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5740 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5742 msave_common_uninit (&common);
5744 msave->common = s->common;
5745 if (!msave->varname_)
5746 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5750 msave_common_uninit (&common);
5751 matrix_cmd_destroy (cmd);
5756 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5758 gsl_matrix *m = matrix_expr_evaluate (e);
5764 msg (SE, _("%s expression must evaluate to vector, "
5765 "not a %zu×%zu matrix."),
5766 name, m->size1, m->size2);
5767 gsl_matrix_free (m);
5771 return matrix_to_vector (m);
5775 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5777 for (size_t i = 0; i < vars->n; i++)
5778 if (!dict_create_var (d, vars->strings[i], 0))
5779 return vars->strings[i];
5783 static struct dictionary *
5784 msave_create_dict (const struct msave_common *common)
5786 struct dictionary *dict = dict_create (get_default_encoding ());
5788 const char *dup_split = msave_add_vars (dict, &common->snames);
5791 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5795 dict_create_var_assert (dict, "ROWTYPE_", 8);
5797 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5800 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5804 dict_create_var_assert (dict, "VARNAME_", 8);
5806 const char *dup_var = msave_add_vars (dict, &common->variables);
5809 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5821 matrix_cmd_execute_msave (struct msave_command *msave)
5823 struct msave_common *common = msave->common;
5824 gsl_matrix *m = NULL;
5825 gsl_vector *factors = NULL;
5826 gsl_vector *splits = NULL;
5828 m = matrix_expr_evaluate (msave->expr);
5832 if (!common->variables.n)
5833 for (size_t i = 0; i < m->size2; i++)
5834 string_array_append_nocopy (&common->variables,
5835 xasprintf ("COL%zu", i + 1));
5837 if (m->size2 != common->variables.n)
5839 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5840 m->size2, common->variables.n);
5846 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5850 if (!common->fnames.n)
5851 for (size_t i = 0; i < factors->size; i++)
5852 string_array_append_nocopy (&common->fnames,
5853 xasprintf ("FAC%zu", i + 1));
5855 if (factors->size != common->fnames.n)
5857 msg (SE, _("There are %zu factor variables, "
5858 "but %zu split values were supplied."),
5859 common->fnames.n, factors->size);
5866 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5870 if (!common->fnames.n)
5871 for (size_t i = 0; i < splits->size; i++)
5872 string_array_append_nocopy (&common->fnames,
5873 xasprintf ("SPL%zu", i + 1));
5875 if (splits->size != common->fnames.n)
5877 msg (SE, _("There are %zu split variables, "
5878 "but %zu split values were supplied."),
5879 common->fnames.n, splits->size);
5884 if (!common->writer)
5886 struct dictionary *dict = msave_create_dict (common);
5890 common->writer = any_writer_open (common->outfile, dict);
5891 if (!common->writer)
5897 common->dict = dict;
5900 for (size_t y = 0; y < m->size1; y++)
5902 struct ccase *c = case_create (dict_get_proto (common->dict));
5905 /* Split variables */
5907 for (size_t i = 0; i < splits->size; i++)
5908 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5911 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5912 msave->rowtype, ' ');
5916 for (size_t i = 0; i < factors->size; i++)
5917 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5920 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5921 msave->varname_, ' ');
5923 /* Continuous variables. */
5924 for (size_t x = 0; x < m->size2; x++)
5925 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5926 casewriter_write (common->writer, c);
5930 gsl_matrix_free (m);
5931 gsl_vector_free (factors);
5932 gsl_vector_free (splits);
5935 static struct matrix_cmd *
5936 matrix_parse_mget (struct matrix_state *s)
5938 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5939 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5941 struct mget_command *mget = &cmd->mget;
5945 if (lex_match_id (s->lexer, "FILE"))
5947 lex_match (s->lexer, T_EQUALS);
5949 fh_unref (mget->file);
5950 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5954 else if (lex_match_id (s->lexer, "TYPE"))
5956 lex_match (s->lexer, T_EQUALS);
5957 while (lex_token (s->lexer) != T_SLASH
5958 && lex_token (s->lexer) != T_ENDCMD)
5960 const char *rowtype = match_rowtype (s->lexer);
5964 stringi_set_insert (&mget->rowtypes, rowtype);
5969 lex_error_expecting (s->lexer, "FILE", "TYPE");
5972 if (lex_token (s->lexer) == T_ENDCMD)
5975 if (!lex_force_match (s->lexer, T_SLASH))
5981 matrix_cmd_destroy (cmd);
5985 static const struct variable *
5986 get_a8_var (const struct dictionary *d, const char *name)
5988 const struct variable *v = dict_lookup_var (d, name);
5991 msg (SE, _("Matrix data file lacks %s variable."), name);
5994 if (var_get_width (v) != 8)
5996 msg (SE, _("%s variable in matrix data file must be "
5997 "8-byte string, but it has width %d."),
5998 name, var_get_width (v));
6005 is_rowtype (const union value *v, const char *rowtype)
6007 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6008 ss_rtrim (&vs, ss_cstr (" "));
6009 return ss_equals_case (vs, ss_cstr (rowtype));
6013 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6014 const struct dictionary *d,
6015 const struct variable *rowtype_var,
6016 struct matrix_state *s, size_t si, size_t fi,
6017 size_t cs, size_t cn)
6022 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6023 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6024 : is_rowtype (rowtype_, "CORR") ? "CR"
6025 : is_rowtype (rowtype_, "MEAN") ? "MN"
6026 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6027 : is_rowtype (rowtype_, "N") ? "NC"
6030 struct string name = DS_EMPTY_INITIALIZER;
6031 ds_put_cstr (&name, name_prefix);
6033 ds_put_format (&name, "F%zu", fi);
6035 ds_put_format (&name, "S%zu", si);
6037 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6039 mv = matrix_var_create (s, ds_ss (&name));
6042 msg (SW, _("Matrix data file contains variable with existing name %s."),
6047 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6048 size_t n_missing = 0;
6049 for (size_t y = 0; y < n_rows; y++)
6051 for (size_t x = 0; x < cn; x++)
6053 struct variable *var = dict_get_var (d, cs + x);
6054 double value = case_num (rows[y], var);
6055 if (var_is_num_missing (var, value, MV_ANY))
6060 gsl_matrix_set (m, y, x, value);
6065 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6066 "value, which was treated as zero.",
6067 "Matrix data file variable %s contains %zu missing "
6068 "values, which were treated as zero.", n_missing),
6069 ds_cstr (&name), n_missing);
6074 for (size_t y = 0; y < n_rows; y++)
6075 case_unref (rows[y]);
6079 var_changed (const struct ccase *ca, const struct ccase *cb,
6080 const struct variable *var)
6082 return !value_equal (case_data (ca, var), case_data (cb, var),
6083 var_get_width (var));
6087 vars_changed (const struct ccase *ca, const struct ccase *cb,
6088 const struct dictionary *d,
6089 size_t first_var, size_t n_vars)
6091 for (size_t i = 0; i < n_vars; i++)
6093 const struct variable *v = dict_get_var (d, first_var + i);
6094 if (var_changed (ca, cb, v))
6101 matrix_cmd_execute_mget (struct mget_command *mget)
6103 struct dictionary *d;
6104 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6109 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6110 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6111 if (!rowtype_ || !varname_)
6114 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6116 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6119 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6121 msg (SE, _("Matrix data file contains no continuous variables."));
6125 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6127 const struct variable *v = dict_get_var (d, i);
6128 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6131 _("Matrix data file contains unexpected string variable %s."),
6137 /* SPLIT variables. */
6139 size_t sn = var_get_dict_index (rowtype_);
6140 struct ccase *sc = NULL;
6143 /* FACTOR variables. */
6144 size_t fs = var_get_dict_index (rowtype_) + 1;
6145 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6146 struct ccase *fc = NULL;
6149 /* Continuous variables. */
6150 size_t cs = var_get_dict_index (varname_) + 1;
6151 size_t cn = dict_get_var_cnt (d) - cs;
6152 struct ccase *cc = NULL;
6155 struct ccase **rows = NULL;
6156 size_t allocated_rows = 0;
6160 while ((c = casereader_read (r)) != NULL)
6162 bool sd = vars_changed (sc, c, d, ss, sn);
6163 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6164 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6179 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6180 mget->state, si, fi, cs, cn);
6186 if (n_rows >= allocated_rows)
6187 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6190 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6191 mget->state, si, fi, cs, cn);
6195 casereader_destroy (r);
6199 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6201 if (!lex_force_id (s->lexer))
6204 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6206 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6211 static struct matrix_cmd *
6212 matrix_parse_eigen (struct matrix_state *s)
6214 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6215 *cmd = (struct matrix_cmd) {
6217 .eigen = { .expr = NULL }
6220 struct eigen_command *eigen = &cmd->eigen;
6221 if (!lex_force_match (s->lexer, T_LPAREN))
6223 eigen->expr = matrix_parse_expr (s);
6225 || !lex_force_match (s->lexer, T_COMMA)
6226 || !matrix_parse_dst_var (s, &eigen->evec)
6227 || !lex_force_match (s->lexer, T_COMMA)
6228 || !matrix_parse_dst_var (s, &eigen->eval)
6229 || !lex_force_match (s->lexer, T_RPAREN))
6235 matrix_cmd_destroy (cmd);
6240 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6242 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6245 if (!is_symmetric (A))
6247 msg (SE, _("Argument of EIGEN must be symmetric."));
6248 gsl_matrix_free (A);
6252 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6253 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6254 gsl_vector v_eval = to_vector (eval);
6255 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6256 gsl_eigen_symmv (A, &v_eval, evec, w);
6257 gsl_eigen_symmv_free (w);
6259 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6261 gsl_matrix_free (A);
6263 gsl_matrix_free (eigen->eval->value);
6264 eigen->eval->value = eval;
6266 gsl_matrix_free (eigen->evec->value);
6267 eigen->evec->value = evec;
6270 static struct matrix_cmd *
6271 matrix_parse_setdiag (struct matrix_state *s)
6273 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6274 *cmd = (struct matrix_cmd) {
6275 .type = MCMD_SETDIAG,
6276 .setdiag = { .dst = NULL }
6279 struct setdiag_command *setdiag = &cmd->setdiag;
6280 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6283 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6286 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6291 if (!lex_force_match (s->lexer, T_COMMA))
6294 setdiag->expr = matrix_parse_expr (s);
6298 if (!lex_force_match (s->lexer, T_RPAREN))
6304 matrix_cmd_destroy (cmd);
6309 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6311 gsl_matrix *dst = setdiag->dst->value;
6314 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6315 setdiag->dst->name);
6319 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6323 size_t n = MIN (dst->size1, dst->size2);
6324 if (is_scalar (src))
6326 double d = to_scalar (src);
6327 for (size_t i = 0; i < n; i++)
6328 gsl_matrix_set (dst, i, i, d);
6330 else if (is_vector (src))
6332 gsl_vector v = to_vector (src);
6333 for (size_t i = 0; i < n && i < v.size; i++)
6334 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6337 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6338 "not a %zu×%zu matrix."),
6339 src->size1, src->size2);
6340 gsl_matrix_free (src);
6343 static struct matrix_cmd *
6344 matrix_parse_svd (struct matrix_state *s)
6346 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6347 *cmd = (struct matrix_cmd) {
6349 .svd = { .expr = NULL }
6352 struct svd_command *svd = &cmd->svd;
6353 if (!lex_force_match (s->lexer, T_LPAREN))
6355 svd->expr = matrix_parse_expr (s);
6357 || !lex_force_match (s->lexer, T_COMMA)
6358 || !matrix_parse_dst_var (s, &svd->u)
6359 || !lex_force_match (s->lexer, T_COMMA)
6360 || !matrix_parse_dst_var (s, &svd->s)
6361 || !lex_force_match (s->lexer, T_COMMA)
6362 || !matrix_parse_dst_var (s, &svd->v)
6363 || !lex_force_match (s->lexer, T_RPAREN))
6369 matrix_cmd_destroy (cmd);
6374 matrix_cmd_execute_svd (struct svd_command *svd)
6376 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6380 if (m->size1 >= m->size2)
6383 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6384 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6385 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6386 gsl_vector *work = gsl_vector_alloc (A->size2);
6387 gsl_linalg_SV_decomp (A, V, &Sv, work);
6388 gsl_vector_free (work);
6390 matrix_var_set (svd->u, A);
6391 matrix_var_set (svd->s, S);
6392 matrix_var_set (svd->v, V);
6396 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6397 gsl_matrix_transpose_memcpy (At, m);
6398 gsl_matrix_free (m);
6400 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6401 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6402 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6403 gsl_vector *work = gsl_vector_alloc (At->size2);
6404 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6405 gsl_vector_free (work);
6407 matrix_var_set (svd->v, At);
6408 matrix_var_set (svd->s, St);
6409 matrix_var_set (svd->u, Vt);
6414 matrix_cmd_execute (struct matrix_cmd *cmd)
6419 matrix_cmd_execute_compute (&cmd->compute);
6423 matrix_cmd_execute_print (&cmd->print);
6427 return matrix_cmd_execute_do_if (&cmd->do_if);
6430 matrix_cmd_execute_loop (&cmd->loop);
6437 matrix_cmd_execute_display (&cmd->display);
6441 matrix_cmd_execute_release (&cmd->release);
6445 matrix_cmd_execute_save (&cmd->save);
6449 matrix_cmd_execute_read (&cmd->read);
6453 matrix_cmd_execute_write (&cmd->write);
6457 matrix_cmd_execute_get (&cmd->get);
6461 matrix_cmd_execute_msave (&cmd->msave);
6465 matrix_cmd_execute_mget (&cmd->mget);
6469 matrix_cmd_execute_eigen (&cmd->eigen);
6473 matrix_cmd_execute_setdiag (&cmd->setdiag);
6477 matrix_cmd_execute_svd (&cmd->svd);
6485 matrix_cmds_uninit (struct matrix_cmds *cmds)
6487 for (size_t i = 0; i < cmds->n; i++)
6488 matrix_cmd_destroy (cmds->commands[i]);
6489 free (cmds->commands);
6493 matrix_cmd_destroy (struct matrix_cmd *cmd)
6501 matrix_lvalue_destroy (cmd->compute.lvalue);
6502 matrix_expr_destroy (cmd->compute.rvalue);
6506 matrix_expr_destroy (cmd->print.expression);
6507 free (cmd->print.title);
6508 print_labels_destroy (cmd->print.rlabels);
6509 print_labels_destroy (cmd->print.clabels);
6513 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6515 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6516 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6518 free (cmd->do_if.clauses);
6522 matrix_expr_destroy (cmd->loop.start);
6523 matrix_expr_destroy (cmd->loop.end);
6524 matrix_expr_destroy (cmd->loop.increment);
6525 matrix_expr_destroy (cmd->loop.top_condition);
6526 matrix_expr_destroy (cmd->loop.bottom_condition);
6527 matrix_cmds_uninit (&cmd->loop.commands);
6537 free (cmd->release.vars);
6541 matrix_expr_destroy (cmd->save.expression);
6542 fh_unref (cmd->save.outfile);
6543 string_array_destroy (cmd->save.variables);
6544 matrix_expr_destroy (cmd->save.names);
6545 stringi_set_destroy (&cmd->save.strings);
6549 matrix_lvalue_destroy (cmd->read.dst);
6550 matrix_expr_destroy (cmd->read.size);
6554 matrix_expr_destroy (cmd->write.expression);
6558 matrix_lvalue_destroy (cmd->get.dst);
6559 fh_unref (cmd->get.file);
6560 free (cmd->get.encoding);
6561 string_array_destroy (&cmd->get.variables);
6565 free (cmd->msave.varname_);
6566 matrix_expr_destroy (cmd->msave.expr);
6567 matrix_expr_destroy (cmd->msave.factors);
6568 matrix_expr_destroy (cmd->msave.splits);
6572 fh_unref (cmd->mget.file);
6573 stringi_set_destroy (&cmd->mget.rowtypes);
6577 matrix_expr_destroy (cmd->eigen.expr);
6581 matrix_expr_destroy (cmd->setdiag.expr);
6585 matrix_expr_destroy (cmd->svd.expr);
6591 struct matrix_command_name
6594 struct matrix_cmd *(*parse) (struct matrix_state *);
6597 static const struct matrix_command_name *
6598 matrix_parse_command_name (struct lexer *lexer)
6600 static const struct matrix_command_name commands[] = {
6601 { "COMPUTE", matrix_parse_compute },
6602 { "CALL EIGEN", matrix_parse_eigen },
6603 { "CALL SETDIAG", matrix_parse_setdiag },
6604 { "CALL SVD", matrix_parse_svd },
6605 { "PRINT", matrix_parse_print },
6606 { "DO IF", matrix_parse_do_if },
6607 { "LOOP", matrix_parse_loop },
6608 { "BREAK", matrix_parse_break },
6609 { "READ", matrix_parse_read },
6610 { "WRITE", matrix_parse_write },
6611 { "GET", matrix_parse_get },
6612 { "SAVE", matrix_parse_save },
6613 { "MGET", matrix_parse_mget },
6614 { "MSAVE", matrix_parse_msave },
6615 { "DISPLAY", matrix_parse_display },
6616 { "RELEASE", matrix_parse_release },
6618 static size_t n = sizeof commands / sizeof *commands;
6620 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6621 if (lex_match_phrase (lexer, c->name))
6626 static struct matrix_cmd *
6627 matrix_parse_command (struct matrix_state *s)
6629 size_t nesting_level = SIZE_MAX;
6631 struct matrix_cmd *c = NULL;
6632 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6634 lex_error (s->lexer, _("Unknown matrix command."));
6635 else if (!cmd->parse)
6636 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6640 nesting_level = output_open_group (
6641 group_item_create_nocopy (utf8_to_title (cmd->name),
6642 utf8_to_title (cmd->name)));
6647 lex_end_of_command (s->lexer);
6648 lex_discard_rest_of_command (s->lexer);
6649 if (nesting_level != SIZE_MAX)
6650 output_close_groups (nesting_level);
6656 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6658 if (!lex_force_match (lexer, T_ENDCMD))
6661 struct matrix_state state = {
6662 .session = dataset_session (ds),
6664 .vars = HMAP_INITIALIZER (state.vars),
6669 while (lex_match (lexer, T_ENDCMD))
6671 if (lex_token (lexer) == T_STOP)
6673 msg (SE, _("Unexpected end of input expecting matrix command."));
6677 if (lex_match_phrase (lexer, "END MATRIX"))
6680 struct matrix_cmd *c = matrix_parse_command (&state);
6683 matrix_cmd_execute (c);
6684 matrix_cmd_destroy (c);
6688 struct matrix_var *var, *next;
6689 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6692 gsl_matrix_free (var->value);
6693 hmap_delete (&state.vars, &var->hmap_node);
6696 hmap_destroy (&state.vars);
6697 fh_unref (state.prev_save_outfile);
6700 dict_unref (state.common->dict);
6701 casewriter_destroy (state.common->writer);
6702 free (state.common);
6704 fh_unref (state.prev_read_file);
6705 for (size_t i = 0; i < state.n_read_files; i++)
6706 read_file_destroy (state.read_files[i]);
6707 free (state.read_files);
6708 fh_unref (state.prev_write_file);
6709 for (size_t i = 0; i < state.n_write_files; i++)
6710 write_file_destroy (state.write_files[i]);
6711 free (state.write_files);