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 dataset *dataset;
102 struct session *session;
106 struct file_handle *prev_save_outfile;
107 struct file_handle *prev_write_outfile;
108 struct msave_common *common;
110 struct file_handle *prev_read_file;
111 struct read_file **read_files;
115 static struct matrix_var *
116 matrix_var_lookup (struct matrix_state *s, struct substring name)
118 struct matrix_var *var;
120 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
121 utf8_hash_case_substring (name, 0), &s->vars)
122 if (!utf8_sscasecmp (ss_cstr (var->name), name))
128 static struct matrix_var *
129 matrix_var_create (struct matrix_state *s, struct substring name)
131 struct matrix_var *var = xmalloc (sizeof *var);
132 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
133 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
138 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
140 gsl_matrix_free (var->value);
144 #define MATRIX_FUNCTIONS \
200 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
201 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
202 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
203 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
204 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
205 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
206 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
207 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
208 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
209 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
210 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
212 typedef gsl_matrix *matrix_proto_m_d (double);
213 typedef gsl_matrix *matrix_proto_m_dd (double, double);
214 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
215 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
216 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
217 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
218 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
219 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
220 typedef double matrix_proto_d_m (gsl_matrix *);
221 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
222 typedef gsl_matrix *matrix_proto_IDENT (double, double);
224 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
228 /* Matrix expressions. */
235 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
239 /* Elementwise and scalar arithmetic. */
240 MOP_NEGATE, /* unary - */
241 MOP_ADD_ELEMS, /* + */
242 MOP_SUB_ELEMS, /* - */
243 MOP_MUL_ELEMS, /* &* */
244 MOP_DIV_ELEMS, /* / and &/ */
245 MOP_EXP_ELEMS, /* &** */
247 MOP_SEQ_BY, /* a:b:c */
249 /* Matrix arithmetic. */
251 MOP_EXP_MAT, /* ** */
268 MOP_PASTE_HORZ, /* a, b, c, ... */
269 MOP_PASTE_VERT, /* a; b; c; ... */
273 MOP_VEC_INDEX, /* x(y) */
274 MOP_VEC_ALL, /* x(:) */
275 MOP_MAT_INDEX, /* x(y,z) */
276 MOP_ROW_INDEX, /* x(y,:) */
277 MOP_COL_INDEX, /* x(:,z) */
284 MOP_EOF, /* EOF('file') */
292 struct matrix_expr **subs;
297 struct matrix_var *variable;
298 struct read_file *eof;
303 matrix_expr_destroy (struct matrix_expr *e)
310 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
341 for (size_t i = 0; i < e->n_subs; i++)
342 matrix_expr_destroy (e->subs[i]);
353 static struct matrix_expr *
354 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
357 struct matrix_expr *e = xmalloc (sizeof *e);
358 *e = (struct matrix_expr) {
360 .subs = xmemdup (subs, n_subs * sizeof *subs),
366 static struct matrix_expr *
367 matrix_expr_create_0 (enum matrix_op op)
369 struct matrix_expr *sub;
370 return matrix_expr_create_subs (op, &sub, 0);
373 static struct matrix_expr *
374 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
376 return matrix_expr_create_subs (op, &sub, 1);
379 static struct matrix_expr *
380 matrix_expr_create_2 (enum matrix_op op,
381 struct matrix_expr *sub0, struct matrix_expr *sub1)
383 struct matrix_expr *subs[] = { sub0, sub1 };
384 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
387 static struct matrix_expr *
388 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
389 struct matrix_expr *sub1, struct matrix_expr *sub2)
391 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
392 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
395 static struct matrix_expr *
396 matrix_expr_create_number (double number)
398 struct matrix_expr *e = xmalloc (sizeof *e);
399 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
403 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
405 static struct matrix_expr *
406 matrix_parse_curly_comma (struct matrix_state *s)
408 struct matrix_expr *lhs = matrix_parse_or_xor (s);
412 while (lex_match (s->lexer, T_COMMA))
414 struct matrix_expr *rhs = matrix_parse_or_xor (s);
417 matrix_expr_destroy (lhs);
420 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
425 static struct matrix_expr *
426 matrix_parse_curly_semi (struct matrix_state *s)
428 if (lex_token (s->lexer) == T_RCURLY)
429 return matrix_expr_create_0 (MOP_EMPTY);
431 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
435 while (lex_match (s->lexer, T_SEMICOLON))
437 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
440 matrix_expr_destroy (lhs);
443 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
448 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
449 for (size_t Y = 0; Y < (M)->size1; Y++) \
450 for (size_t X = 0; X < (M)->size2; X++) \
451 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
454 is_vector (gsl_matrix *m)
456 return m->size1 <= 1 || m->size2 <= 1;
460 to_vector (gsl_matrix *m)
462 return (m->size1 == 1
463 ? gsl_matrix_row (m, 0).vector
464 : gsl_matrix_column (m, 0).vector);
469 matrix_eval_ABS (gsl_matrix *m)
471 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
477 matrix_eval_ALL (gsl_matrix *m)
479 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
486 matrix_eval_ANY (gsl_matrix *m)
488 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
495 matrix_eval_ARSIN (gsl_matrix *m)
497 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
503 matrix_eval_ARTAN (gsl_matrix *m)
505 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
511 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
515 for (size_t i = 0; i < n; i++)
520 gsl_matrix *block = gsl_matrix_calloc (r, c);
522 for (size_t i = 0; i < n; i++)
524 for (size_t y = 0; y < m[i]->size1; y++)
525 for (size_t x = 0; x < m[i]->size2; x++)
526 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
534 matrix_eval_CHOL (gsl_matrix *m)
536 if (!gsl_linalg_cholesky_decomp1 (m))
538 for (size_t y = 0; y < m->size1; y++)
539 for (size_t x = y + 1; x < m->size2; x++)
540 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
542 for (size_t y = 0; y < m->size1; y++)
543 for (size_t x = 0; x < y; x++)
544 gsl_matrix_set (m, y, x, 0);
549 msg (SE, _("Input to CHOL function is not positive-definite."));
555 matrix_eval_col_extremum (gsl_matrix *m, bool min)
560 return gsl_matrix_alloc (1, 0);
562 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
563 for (size_t x = 0; x < m->size2; x++)
565 double ext = gsl_matrix_get (m, 0, x);
566 for (size_t y = 1; y < m->size1; y++)
568 double value = gsl_matrix_get (m, y, x);
569 if (min ? value < ext : value > ext)
572 gsl_matrix_set (cext, 0, x, ext);
578 matrix_eval_CMAX (gsl_matrix *m)
580 return matrix_eval_col_extremum (m, false);
584 matrix_eval_CMIN (gsl_matrix *m)
586 return matrix_eval_col_extremum (m, true);
590 matrix_eval_COS (gsl_matrix *m)
592 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
598 matrix_eval_col_sum (gsl_matrix *m, bool square)
603 return gsl_matrix_alloc (1, 0);
605 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
606 for (size_t x = 0; x < m->size2; x++)
609 for (size_t y = 0; y < m->size1; y++)
611 double d = gsl_matrix_get (m, y, x);
612 sum += square ? pow2 (d) : d;
614 gsl_matrix_set (result, 0, x, sum);
620 matrix_eval_CSSQ (gsl_matrix *m)
622 return matrix_eval_col_sum (m, true);
626 matrix_eval_CSUM (gsl_matrix *m)
628 return matrix_eval_col_sum (m, false);
632 compare_double_3way (const void *a_, const void *b_)
634 const double *a = a_;
635 const double *b = b_;
636 return *a < *b ? -1 : *a > *b;
640 matrix_eval_DESIGN (gsl_matrix *m)
642 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
643 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
644 gsl_matrix_transpose_memcpy (&m2, m);
646 for (size_t y = 0; y < m2.size1; y++)
647 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
649 size_t *n = xcalloc (m2.size1, sizeof *n);
651 for (size_t i = 0; i < m2.size1; i++)
653 double *row = tmp + m2.size2 * i;
654 for (size_t j = 0; j < m2.size2; )
657 for (k = j + 1; k < m2.size2; k++)
658 if (row[j] != row[k])
660 row[n[i]++] = row[j];
665 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
670 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
672 for (size_t i = 0; i < m->size2; i++)
677 const double *unique = tmp + m2.size2 * i;
678 for (size_t j = 0; j < n[i]; j++, x++)
680 double value = unique[j];
681 for (size_t y = 0; y < m->size1; y++)
682 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
693 matrix_eval_DET (gsl_matrix *m)
695 gsl_permutation *p = gsl_permutation_alloc (m->size1);
697 gsl_linalg_LU_decomp (m, p, &signum);
698 gsl_permutation_free (p);
699 return gsl_linalg_LU_det (m, signum);
703 matrix_eval_DIAG (gsl_matrix *m)
705 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
706 for (size_t i = 0; i < diag->size1; i++)
707 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
712 is_symmetric (const gsl_matrix *m)
714 if (m->size1 != m->size2)
717 for (size_t y = 0; y < m->size1; y++)
718 for (size_t x = 0; x < y; x++)
719 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
726 compare_double_desc (const void *a_, const void *b_)
728 const double *a = a_;
729 const double *b = b_;
730 return *a > *b ? -1 : *a < *b;
734 matrix_eval_EVAL (gsl_matrix *m)
736 if (!is_symmetric (m))
738 msg (SE, _("Argument of EVAL must be symmetric."));
742 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
743 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
744 gsl_vector v_eval = to_vector (eval);
745 gsl_eigen_symm (m, &v_eval, w);
746 gsl_eigen_symm_free (w);
748 assert (v_eval.stride == 1);
749 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
755 matrix_eval_EXP (gsl_matrix *m)
757 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
762 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
765 Charl Linssen <charl@itfromb.it>
769 matrix_eval_GINV (gsl_matrix *A)
774 gsl_matrix *tmp_mat = NULL;
777 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
778 tmp_mat = gsl_matrix_alloc (m, n);
779 gsl_matrix_transpose_memcpy (tmp_mat, A);
787 gsl_matrix *V = gsl_matrix_alloc (m, m);
788 gsl_vector *u = gsl_vector_alloc (m);
790 gsl_vector *tmp_vec = gsl_vector_alloc (m);
791 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
792 gsl_vector_free (tmp_vec);
795 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
796 gsl_matrix_set_zero (Sigma_pinv);
797 double cutoff = 1e-15 * gsl_vector_max (u);
799 for (size_t i = 0; i < m; ++i)
801 double x = gsl_vector_get (u, i);
802 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
805 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
806 gsl_matrix *U = gsl_matrix_calloc (n, n);
807 for (size_t i = 0; i < n; i++)
808 for (size_t j = 0; j < m; j++)
809 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
811 /* two dot products to obtain pseudoinverse */
812 tmp_mat = gsl_matrix_alloc (m, n);
813 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat);
818 A_pinv = gsl_matrix_alloc (n, m);
819 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat, 0., A_pinv);
823 A_pinv = gsl_matrix_alloc (m, n);
824 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat, U, 0., A_pinv);
827 gsl_matrix_free (tmp_mat);
829 gsl_matrix_free (Sigma_pinv);
843 grade_compare_3way (const void *a_, const void *b_)
845 const struct grade *a = a_;
846 const struct grade *b = b_;
848 return (a->value < b->value ? -1
849 : a->value > b->value ? 1
857 matrix_eval_GRADE (gsl_matrix *m)
859 size_t n = m->size1 * m->size2;
860 struct grade *grades = xmalloc (n * sizeof *grades);
863 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
864 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
865 qsort (grades, n, sizeof *grades, grade_compare_3way);
867 for (size_t i = 0; i < n; i++)
868 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
876 dot (gsl_vector *a, gsl_vector *b)
879 for (size_t i = 0; i < a->size; i++)
880 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
885 norm2 (gsl_vector *v)
888 for (size_t i = 0; i < v->size; i++)
889 result += pow2 (gsl_vector_get (v, i));
896 return sqrt (norm2 (v));
900 matrix_eval_GSCH (gsl_matrix *v)
902 if (v->size2 < v->size1)
904 msg (SE, _("GSCH requires its argument to have at least as many columns "
905 "as rows, but it has dimensions (%zu,%zu)."),
909 if (!v->size1 || !v->size2)
912 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
914 for (size_t vx = 0; vx < v->size2; vx++)
916 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
917 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
919 gsl_vector_memcpy (&u_i, &v_i);
920 for (size_t j = 0; j < ux; j++)
922 gsl_vector u_j = gsl_matrix_column (u, j).vector;
923 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
924 for (size_t k = 0; k < u_i.size; k++)
925 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
926 - scale * gsl_vector_get (&u_j, k)));
929 double len = norm (&u_i);
932 gsl_vector_scale (&u_i, 1.0 / len);
933 if (++ux >= v->size1)
940 msg (SE, _("Argument to GSCH with dimensions (%zu,%zu) contains only "
941 "%zu linearly independent columns."),
942 v->size1, v->size2, ux);
951 matrix_eval_IDENT (double s1, double s2)
953 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
955 msg (SE, _("Arguments to IDENT must be non-negative integers."));
959 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
960 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
966 invert_matrix (gsl_matrix *x)
968 gsl_permutation *p = gsl_permutation_alloc (x->size1);
970 gsl_linalg_LU_decomp (x, p, &signum);
971 gsl_linalg_LU_invx (x, p);
972 gsl_permutation_free (p);
976 matrix_eval_INV (gsl_matrix *m)
983 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
985 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
986 a->size2 * b->size2);
988 for (size_t ar = 0; ar < a->size1; ar++)
989 for (size_t br = 0; br < b->size1; br++, y++)
992 for (size_t ac = 0; ac < a->size2; ac++)
993 for (size_t bc = 0; bc < b->size2; bc++, x++)
995 double av = gsl_matrix_get (a, ar, ac);
996 double bv = gsl_matrix_get (b, br, bc);
997 gsl_matrix_set (k, y, x, av * bv);
1004 matrix_eval_LG10 (gsl_matrix *m)
1006 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1012 matrix_eval_LN (gsl_matrix *m)
1014 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1020 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1022 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1025 for (size_t i = 1; i <= n * n; i++)
1027 gsl_matrix_set (m, y, x, i);
1029 size_t y1 = !y ? n - 1 : y - 1;
1030 size_t x1 = x + 1 >= n ? 0 : x + 1;
1031 if (gsl_matrix_get (m, y1, x1) == 0)
1037 y = y + 1 >= n ? 0 : y + 1;
1042 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1044 double a = gsl_matrix_get (m, y1, x1);
1045 double b = gsl_matrix_get (m, y2, x2);
1046 gsl_matrix_set (m, y1, x1, b);
1047 gsl_matrix_set (m, y2, x2, a);
1051 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1055 /* A. Umar, "On the Construction of Even Order Magic Squares",
1056 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1058 for (size_t i = 1; i <= n * n / 2; i++)
1060 gsl_matrix_set (m, y, x, i);
1070 for (size_t i = n * n; i > n * n / 2; i--)
1072 gsl_matrix_set (m, y, x, i);
1080 for (size_t y = 0; y < n; y++)
1081 for (size_t x = 0; x < n / 2; x++)
1083 unsigned int d = gsl_matrix_get (m, y, x);
1084 if (d % 2 != (y < n / 2))
1085 magic_exchange (m, y, x, y, n - x - 1);
1090 size_t x1 = n / 2 - 1;
1092 magic_exchange (m, y1, x1, y2, x1);
1093 magic_exchange (m, y1, x2, y2, x2);
1097 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1099 /* A. Umar, "On the Construction of Even Order Magic Squares",
1100 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1104 for (size_t i = 1; ; i++)
1106 gsl_matrix_set (m, y, x, i);
1107 if (++y == n / 2 - 1)
1119 for (size_t i = n * n; ; i--)
1121 gsl_matrix_set (m, y, x, i);
1122 if (++y == n / 2 - 1)
1131 for (size_t y = 0; y < n; y++)
1132 if (y != n / 2 - 1 && y != n / 2)
1133 for (size_t x = 0; x < n / 2; x++)
1135 unsigned int d = gsl_matrix_get (m, y, x);
1136 if (d % 2 != (y < n / 2))
1137 magic_exchange (m, y, x, y, n - x - 1);
1140 size_t a0 = (n * n - 2 * n) / 2 + 1;
1141 for (size_t i = 0; i < n / 2; i++)
1144 gsl_matrix_set (m, n / 2, i, a);
1145 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1147 for (size_t i = 0; i < n / 2; i++)
1149 size_t a = a0 + i + n / 2;
1150 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1151 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1153 for (size_t x = 1; x < n / 2; x += 2)
1154 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1155 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1156 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1157 size_t x1 = n / 2 - 2;
1158 size_t x2 = n / 2 + 1;
1159 size_t y1 = n / 2 - 2;
1160 size_t y2 = n / 2 + 1;
1161 magic_exchange (m, y1, x1, y2, x1);
1162 magic_exchange (m, y1, x2, y2, x2);
1166 matrix_eval_MAGIC (double n_)
1168 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1170 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1175 gsl_matrix *m = gsl_matrix_calloc (n, n);
1177 matrix_eval_MAGIC_odd (m, n);
1179 matrix_eval_MAGIC_singly_even (m, n);
1181 matrix_eval_MAGIC_doubly_even (m, n);
1186 matrix_eval_MAKE (double r, double c, double value)
1188 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1190 msg (SE, _("Size arguments to MAKE must be integers."));
1194 gsl_matrix *m = gsl_matrix_alloc (r, c);
1195 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1201 matrix_eval_MDIAG (gsl_vector *v)
1203 gsl_matrix *m = gsl_matrix_alloc (v->size, v->size);
1204 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1205 gsl_vector_memcpy (&diagonal, v);
1210 matrix_eval_MMAX (gsl_matrix *m)
1212 return gsl_matrix_max (m);
1216 matrix_eval_MMIN (gsl_matrix *m)
1218 return gsl_matrix_min (m);
1222 matrix_eval_MOD (gsl_matrix *m, double divisor)
1226 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1230 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1231 *d = fmod (*d, divisor);
1236 matrix_eval_MSSQ (gsl_matrix *m)
1239 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1245 matrix_eval_MSUM (gsl_matrix *m)
1248 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1254 matrix_eval_NCOL (gsl_matrix *m)
1260 matrix_eval_NROW (gsl_matrix *m)
1266 matrix_eval_RANK (gsl_matrix *m)
1268 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1269 gsl_linalg_QR_decomp (m, tau);
1270 gsl_vector_free (tau);
1272 return gsl_linalg_QRPT_rank (m, -1);
1276 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1278 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1280 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1285 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1287 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1288 "product of matrix dimensions."));
1292 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1295 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1297 gsl_matrix_set (dst, y1, x1, *d);
1308 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1313 return gsl_matrix_alloc (0, 1);
1315 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1316 for (size_t y = 0; y < m->size1; y++)
1318 double ext = gsl_matrix_get (m, y, 0);
1319 for (size_t x = 1; x < m->size2; x++)
1321 double value = gsl_matrix_get (m, y, x);
1322 if (min ? value < ext : value > ext)
1325 gsl_matrix_set (rext, y, 0, ext);
1331 matrix_eval_RMAX (gsl_matrix *m)
1333 return matrix_eval_row_extremum (m, false);
1337 matrix_eval_RMIN (gsl_matrix *m)
1339 return matrix_eval_row_extremum (m, true);
1343 matrix_eval_RND (gsl_matrix *m)
1345 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1357 rank_compare_3way (const void *a_, const void *b_)
1359 const struct rank *a = a_;
1360 const struct rank *b = b_;
1362 return a->value < b->value ? -1 : a->value > b->value;
1366 matrix_eval_RNKORDER (gsl_matrix *m)
1368 size_t n = m->size1 * m->size2;
1369 struct rank *ranks = xmalloc (n * sizeof *ranks);
1371 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1372 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1373 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1375 for (size_t i = 0; i < n; )
1378 for (j = i + 1; j < n; j++)
1379 if (ranks[i].value != ranks[j].value)
1382 double rank = (i + j + 1.0) / 2.0;
1383 for (size_t k = i; k < j; k++)
1384 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1395 matrix_eval_row_sum (gsl_matrix *m, bool square)
1400 return gsl_matrix_alloc (0, 1);
1402 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1403 for (size_t y = 0; y < m->size1; y++)
1406 for (size_t x = 0; x < m->size2; x++)
1408 double d = gsl_matrix_get (m, y, x);
1409 sum += square ? pow2 (d) : d;
1411 gsl_matrix_set (result, y, 0, sum);
1417 matrix_eval_RSSQ (gsl_matrix *m)
1419 return matrix_eval_row_sum (m, true);
1423 matrix_eval_RSUM (gsl_matrix *m)
1425 return matrix_eval_row_sum (m, false);
1429 matrix_eval_SIN (gsl_matrix *m)
1431 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1437 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1439 if (m1->size1 != m2->size1)
1441 msg (SE, _("SOLVE requires its arguments to have the same number of "
1442 "rows, but the first argument has dimensions (%zu,%zu) and "
1443 "the second (%zu,%zu)."),
1444 m1->size1, m1->size2,
1445 m2->size1, m2->size2);
1449 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1450 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1452 gsl_linalg_LU_decomp (m1, p, &signum);
1453 for (size_t i = 0; i < m2->size2; i++)
1455 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1456 gsl_vector xi = gsl_matrix_column (x, i).vector;
1457 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1459 gsl_permutation_free (p);
1464 matrix_eval_SQRT (gsl_matrix *m)
1466 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1470 msg (SE, _("Argument to SQRT must be nonnegative."));
1479 matrix_eval_SSCP (gsl_matrix *m)
1481 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1482 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1487 matrix_eval_SVAL (gsl_matrix *m)
1489 gsl_matrix *tmp_mat = NULL;
1490 if (m->size2 > m->size1)
1492 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1493 gsl_matrix_transpose_memcpy (tmp_mat, m);
1498 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1499 gsl_vector *S = gsl_vector_alloc (m->size2);
1500 gsl_vector *work = gsl_vector_alloc (m->size2);
1501 gsl_linalg_SV_decomp (m, V, S, work);
1503 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1504 for (size_t i = 0; i < m->size2; i++)
1505 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1507 gsl_matrix_free (V);
1508 gsl_vector_free (S);
1509 gsl_vector_free (work);
1510 gsl_matrix_free (tmp_mat);
1516 matrix_eval_SWEEP (gsl_matrix *m, double d)
1518 if (d < 1 || d > SIZE_MAX)
1520 msg (SE, _("Scalar argument to SWEEP must be integer."));
1524 if (k >= MIN (m->size1, m->size2))
1526 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1527 "equal to the smaller of the matrix argument's rows and "
1532 double m_kk = gsl_matrix_get (m, k, k);
1533 if (fabs (m_kk) > 1e-19)
1535 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1536 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1538 double m_ij = gsl_matrix_get (m, i, j);
1539 double m_ik = gsl_matrix_get (m, i, k);
1540 double m_kj = gsl_matrix_get (m, k, j);
1541 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1550 for (size_t i = 0; i < m->size1; i++)
1552 gsl_matrix_set (m, i, k, 0);
1553 gsl_matrix_set (m, k, i, 0);
1560 matrix_eval_TRACE (gsl_matrix *m)
1563 size_t n = MIN (m->size1, m->size2);
1564 for (size_t i = 0; i < n; i++)
1565 sum += gsl_matrix_get (m, i, i);
1570 matrix_eval_T (gsl_matrix *m)
1572 return matrix_eval_TRANSPOS (m);
1576 matrix_eval_TRANSPOS (gsl_matrix *m)
1578 if (m->size1 == m->size2)
1580 gsl_matrix_transpose (m);
1585 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1586 gsl_matrix_transpose_memcpy (t, m);
1592 matrix_eval_TRUNC (gsl_matrix *m)
1594 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1600 matrix_eval_UNIFORM (double r_, double c_)
1602 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1604 msg (SE, _("Arguments to UNIFORM must be integers."));
1609 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1611 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1615 gsl_matrix *m = gsl_matrix_alloc (r, c);
1616 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1617 *d = gsl_ran_flat (get_rng (), 0, 1);
1621 struct matrix_function
1625 size_t min_args, max_args;
1628 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1630 static const struct matrix_function *
1631 matrix_parse_function_name (const char *token)
1633 static const struct matrix_function functions[] =
1635 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1639 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1641 for (size_t i = 0; i < N_FUNCTIONS; i++)
1643 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1644 return &functions[i];
1649 static struct read_file *
1650 read_file_create (struct matrix_state *s, struct file_handle *fh)
1652 for (size_t i = 0; i < s->n_read_files; i++)
1654 struct read_file *rf = s->read_files[i];
1662 struct read_file *rf = xmalloc (sizeof *rf);
1663 *rf = (struct read_file) { .file = fh };
1665 s->read_files = xrealloc (s->read_files,
1666 (s->n_read_files + 1) * sizeof *s->read_files);
1667 s->read_files[s->n_read_files++] = rf;
1671 static struct dfm_reader *
1672 read_file_open (struct read_file *rf)
1675 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1680 matrix_parse_function (struct matrix_state *s, const char *token,
1681 struct matrix_expr **exprp)
1684 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1687 if (lex_match_id (s->lexer, "EOF"))
1690 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1694 if (!lex_force_match (s->lexer, T_RPAREN))
1700 struct read_file *rf = read_file_create (s, fh);
1702 struct matrix_expr *e = xmalloc (sizeof *e);
1703 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1708 const struct matrix_function *f = matrix_parse_function_name (token);
1712 lex_get_n (s->lexer, 2);
1714 struct matrix_expr *e = xmalloc (sizeof *e);
1715 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1717 size_t allocated_subs = 0;
1720 struct matrix_expr *sub = matrix_parse_expr (s);
1724 if (e->n_subs >= allocated_subs)
1725 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1726 e->subs[e->n_subs++] = sub;
1728 while (lex_match (s->lexer, T_COMMA));
1729 if (!lex_force_match (s->lexer, T_RPAREN))
1736 matrix_expr_destroy (e);
1740 static struct matrix_expr *
1741 matrix_parse_primary (struct matrix_state *s)
1743 if (lex_is_number (s->lexer))
1745 double number = lex_number (s->lexer);
1747 return matrix_expr_create_number (number);
1749 else if (lex_is_string (s->lexer))
1751 char string[sizeof (double)];
1752 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1756 memcpy (&number, string, sizeof number);
1757 return matrix_expr_create_number (number);
1759 else if (lex_match (s->lexer, T_LPAREN))
1761 struct matrix_expr *e = matrix_parse_or_xor (s);
1762 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1764 matrix_expr_destroy (e);
1769 else if (lex_match (s->lexer, T_LCURLY))
1771 struct matrix_expr *e = matrix_parse_curly_semi (s);
1772 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1774 matrix_expr_destroy (e);
1779 else if (lex_token (s->lexer) == T_ID)
1781 struct matrix_expr *retval;
1782 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1785 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1788 lex_error (s->lexer, _("Unknown variable %s."),
1789 lex_tokcstr (s->lexer));
1794 struct matrix_expr *e = xmalloc (sizeof *e);
1795 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1798 else if (lex_token (s->lexer) == T_ALL)
1800 struct matrix_expr *retval;
1801 if (matrix_parse_function (s, "ALL", &retval))
1805 lex_error (s->lexer, NULL);
1809 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1812 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1814 if (lex_match (s->lexer, T_COLON))
1821 *indexp = matrix_parse_expr (s);
1822 return *indexp != NULL;
1826 static struct matrix_expr *
1827 matrix_parse_postfix (struct matrix_state *s)
1829 struct matrix_expr *lhs = matrix_parse_primary (s);
1830 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1833 struct matrix_expr *i0;
1834 if (!matrix_parse_index_expr (s, &i0))
1836 matrix_expr_destroy (lhs);
1839 if (lex_match (s->lexer, T_RPAREN))
1841 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1842 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1843 else if (lex_match (s->lexer, T_COMMA))
1845 struct matrix_expr *i1;
1846 if (!matrix_parse_index_expr (s, &i1)
1847 || !lex_force_match (s->lexer, T_RPAREN))
1849 matrix_expr_destroy (lhs);
1850 matrix_expr_destroy (i0);
1851 matrix_expr_destroy (i1);
1854 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1855 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1856 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1861 lex_error_expecting (s->lexer, "`)'", "`,'");
1866 static struct matrix_expr *
1867 matrix_parse_unary (struct matrix_state *s)
1869 if (lex_match (s->lexer, T_DASH))
1871 struct matrix_expr *lhs = matrix_parse_unary (s);
1872 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1874 else if (lex_match (s->lexer, T_PLUS))
1875 return matrix_parse_unary (s);
1877 return matrix_parse_postfix (s);
1880 static struct matrix_expr *
1881 matrix_parse_seq (struct matrix_state *s)
1883 struct matrix_expr *start = matrix_parse_unary (s);
1884 if (!start || !lex_match (s->lexer, T_COLON))
1887 struct matrix_expr *end = matrix_parse_unary (s);
1890 matrix_expr_destroy (start);
1894 if (lex_match (s->lexer, T_COLON))
1896 struct matrix_expr *increment = matrix_parse_unary (s);
1899 matrix_expr_destroy (start);
1900 matrix_expr_destroy (end);
1903 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1906 return matrix_expr_create_2 (MOP_SEQ, start, end);
1909 static struct matrix_expr *
1910 matrix_parse_exp (struct matrix_state *s)
1912 struct matrix_expr *lhs = matrix_parse_seq (s);
1919 if (lex_match (s->lexer, T_EXP))
1921 else if (lex_match_phrase (s->lexer, "&**"))
1926 struct matrix_expr *rhs = matrix_parse_seq (s);
1929 matrix_expr_destroy (lhs);
1932 lhs = matrix_expr_create_2 (op, lhs, rhs);
1936 static struct matrix_expr *
1937 matrix_parse_mul_div (struct matrix_state *s)
1939 struct matrix_expr *lhs = matrix_parse_exp (s);
1946 if (lex_match (s->lexer, T_ASTERISK))
1948 else if (lex_match (s->lexer, T_SLASH))
1950 else if (lex_match_phrase (s->lexer, "&*"))
1952 else if (lex_match_phrase (s->lexer, "&/"))
1957 struct matrix_expr *rhs = matrix_parse_exp (s);
1960 matrix_expr_destroy (lhs);
1963 lhs = matrix_expr_create_2 (op, lhs, rhs);
1967 static struct matrix_expr *
1968 matrix_parse_add_sub (struct matrix_state *s)
1970 struct matrix_expr *lhs = matrix_parse_mul_div (s);
1977 if (lex_match (s->lexer, T_PLUS))
1979 else if (lex_match (s->lexer, T_DASH))
1981 else if (lex_token (s->lexer) == T_NEG_NUM)
1986 struct matrix_expr *rhs = matrix_parse_mul_div (s);
1989 matrix_expr_destroy (lhs);
1992 lhs = matrix_expr_create_2 (op, lhs, rhs);
1996 static struct matrix_expr *
1997 matrix_parse_relational (struct matrix_state *s)
1999 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2006 if (lex_match (s->lexer, T_GT))
2008 else if (lex_match (s->lexer, T_GE))
2010 else if (lex_match (s->lexer, T_LT))
2012 else if (lex_match (s->lexer, T_LE))
2014 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2016 else if (lex_match (s->lexer, T_NE))
2021 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2024 matrix_expr_destroy (lhs);
2027 lhs = matrix_expr_create_2 (op, lhs, rhs);
2031 static struct matrix_expr *
2032 matrix_parse_not (struct matrix_state *s)
2034 if (lex_match (s->lexer, T_NOT))
2036 struct matrix_expr *lhs = matrix_parse_not (s);
2037 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2040 return matrix_parse_relational (s);
2043 static struct matrix_expr *
2044 matrix_parse_and (struct matrix_state *s)
2046 struct matrix_expr *lhs = matrix_parse_not (s);
2050 while (lex_match (s->lexer, T_AND))
2052 struct matrix_expr *rhs = matrix_parse_not (s);
2055 matrix_expr_destroy (lhs);
2058 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2063 static struct matrix_expr *
2064 matrix_parse_or_xor (struct matrix_state *s)
2066 struct matrix_expr *lhs = matrix_parse_and (s);
2073 if (lex_match (s->lexer, T_OR))
2075 else if (lex_match_id (s->lexer, "XOR"))
2080 struct matrix_expr *rhs = matrix_parse_and (s);
2083 matrix_expr_destroy (lhs);
2086 lhs = matrix_expr_create_2 (op, lhs, rhs);
2090 static struct matrix_expr *
2091 matrix_parse_expr (struct matrix_state *s)
2093 return matrix_parse_or_xor (s);
2096 /* Expression evaluation. */
2099 matrix_op_eval (enum matrix_op op, double a, double b)
2103 case MOP_ADD_ELEMS: return a + b;
2104 case MOP_SUB_ELEMS: return a - b;
2105 case MOP_MUL_ELEMS: return a * b;
2106 case MOP_DIV_ELEMS: return a / b;
2107 case MOP_EXP_ELEMS: return pow (a, b);
2108 case MOP_GT: return a > b;
2109 case MOP_GE: return a >= b;
2110 case MOP_LT: return a < b;
2111 case MOP_LE: return a <= b;
2112 case MOP_EQ: return a == b;
2113 case MOP_NE: return a != b;
2114 case MOP_AND: return (a > 0) && (b > 0);
2115 case MOP_OR: return (a > 0) || (b > 0);
2116 case MOP_XOR: return (a > 0) != (b > 0);
2118 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2127 case MOP_PASTE_HORZ:
2128 case MOP_PASTE_VERT:
2144 matrix_op_name (enum matrix_op op)
2148 case MOP_ADD_ELEMS: return "+";
2149 case MOP_SUB_ELEMS: return "-";
2150 case MOP_MUL_ELEMS: return "&*";
2151 case MOP_DIV_ELEMS: return "&/";
2152 case MOP_EXP_ELEMS: return "&**";
2153 case MOP_GT: return ">";
2154 case MOP_GE: return ">=";
2155 case MOP_LT: return "<";
2156 case MOP_LE: return "<=";
2157 case MOP_EQ: return "=";
2158 case MOP_NE: return "<>";
2159 case MOP_AND: return "AND";
2160 case MOP_OR: return "OR";
2161 case MOP_XOR: return "XOR";
2163 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2172 case MOP_PASTE_HORZ:
2173 case MOP_PASTE_VERT:
2189 is_scalar (const gsl_matrix *m)
2191 return m->size1 == 1 && m->size2 == 1;
2195 to_scalar (const gsl_matrix *m)
2197 assert (is_scalar (m));
2198 return gsl_matrix_get (m, 0, 0);
2202 matrix_expr_evaluate_elementwise (enum matrix_op op,
2203 gsl_matrix *a, gsl_matrix *b)
2207 double be = to_scalar (b);
2208 for (size_t r = 0; r < a->size1; r++)
2209 for (size_t c = 0; c < a->size2; c++)
2211 double *ae = gsl_matrix_ptr (a, r, c);
2212 *ae = matrix_op_eval (op, *ae, be);
2216 else if (is_scalar (a))
2218 double ae = to_scalar (a);
2219 for (size_t r = 0; r < b->size1; r++)
2220 for (size_t c = 0; c < b->size2; c++)
2222 double *be = gsl_matrix_ptr (b, r, c);
2223 *be = matrix_op_eval (op, ae, *be);
2227 else if (a->size1 == b->size1 && a->size2 == b->size2)
2229 for (size_t r = 0; r < a->size1; r++)
2230 for (size_t c = 0; c < a->size2; c++)
2232 double *ae = gsl_matrix_ptr (a, r, c);
2233 double be = gsl_matrix_get (b, r, c);
2234 *ae = matrix_op_eval (op, *ae, be);
2240 msg (SE, _("Operands to %s must have the same dimensions or one "
2241 "must be a scalar, not matrices with dimensions (%zu,%zu) "
2243 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2249 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2251 if (is_scalar (a) || is_scalar (b))
2252 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2254 if (a->size2 != b->size1)
2256 msg (SE, _("Matrices with dimensions (%zu,%zu) and (%zu,%zu) are "
2257 "not conformable for multiplication."),
2258 a->size1, a->size2, b->size1, b->size2);
2262 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2263 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2268 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2270 gsl_matrix *tmp = *a;
2276 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2279 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2280 swap_matrix (z, tmp);
2284 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2286 mul_matrix (x, *x, *x, tmp);
2290 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2293 if (x->size1 != x->size2)
2295 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2296 "the left-hand size, not one with dimensions (%zu,%zu)."),
2297 x->size1, x->size2);
2302 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2303 "right-hand side, not a matrix with dimensions (%zu,%zu)."),
2304 b->size1, b->size2);
2307 double bf = to_scalar (b);
2308 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2310 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2311 "or outside the valid range."), bf);
2316 gsl_matrix *tmp = gsl_matrix_alloc (x->size1, x->size2);
2317 gsl_matrix *y = gsl_matrix_alloc (x->size1, x->size2);
2318 gsl_matrix_set_identity (y);
2322 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2325 mul_matrix (&y, x, y, &tmp);
2326 square_matrix (&x, &tmp);
2329 square_matrix (&x, &tmp);
2331 mul_matrix (&y, x, y, &tmp);
2336 gsl_matrix_free (tmp);
2341 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2344 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2346 msg (SE, _("All operands of : operator must be scalars."));
2350 long int start = to_scalar (start_);
2351 long int end = to_scalar (end_);
2352 long int by = by_ ? to_scalar (by_) : 1;
2356 msg (SE, _("The increment operand to : must be nonzero."));
2360 long int n = (end > start && by > 0 ? (end - start + by) / by
2361 : end < start && by < 0 ? (start - end - by) / -by
2363 gsl_matrix *m = gsl_matrix_alloc (1, n);
2364 for (long int i = 0; i < n; i++)
2365 gsl_matrix_set (m, 0, i, start + i * by);
2370 matrix_expr_evaluate_not (gsl_matrix *a)
2372 for (size_t r = 0; r < a->size1; r++)
2373 for (size_t c = 0; c < a->size2; c++)
2375 double *ae = gsl_matrix_ptr (a, r, c);
2382 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2384 if (a->size1 != b->size1)
2386 if (!a->size1 || !a->size2)
2388 else if (!b->size1 || !b->size2)
2391 msg (SE, _("All columns in a matrix must have the same number of rows, "
2392 "but this tries to paste matrices with %zu and %zu rows."),
2393 a->size1, b->size1);
2397 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2398 for (size_t y = 0; y < a->size1; y++)
2400 for (size_t x = 0; x < a->size2; x++)
2401 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2402 for (size_t x = 0; x < b->size2; x++)
2403 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2409 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2411 if (a->size2 != b->size2)
2413 if (!a->size1 || !a->size2)
2415 else if (!b->size1 || !b->size2)
2418 msg (SE, _("All rows in a matrix must have the same number of columns, "
2419 "but this tries to stack matrices with %zu and %zu columns."),
2420 a->size2, b->size2);
2424 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2425 for (size_t x = 0; x < a->size2; x++)
2427 for (size_t y = 0; y < a->size1; y++)
2428 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2429 for (size_t y = 0; y < b->size1; y++)
2430 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2436 matrix_to_vector (gsl_matrix *m)
2439 gsl_vector v = to_vector (m);
2440 assert (v.block == m->block || !v.block);
2444 return xmemdup (&v, sizeof v);
2458 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2461 matrix_normalize_index_vector (gsl_matrix *m, size_t size,
2462 enum index_type index_type, size_t other_size,
2463 struct index_vector *iv)
2472 msg (SE, _("Vector index must be scalar or vector, not a "
2473 "matrix with dimensions (%zu,%zu)."),
2474 m->size1, m->size2);
2478 msg (SE, _("Matrix row index must be scalar or vector, not a "
2479 "matrix with dimensions (%zu,%zu)."),
2480 m->size1, m->size2);
2484 msg (SE, _("Matrix column index must be scalar or vector, not a "
2485 "matrix with dimensions (%zu,%zu)."),
2486 m->size1, m->size2);
2492 gsl_vector v = to_vector (m);
2493 *iv = (struct index_vector) {
2494 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2497 for (size_t i = 0; i < v.size; i++)
2499 double index = gsl_vector_get (&v, i);
2500 if (index < 1 || index >= size + 1)
2505 msg (SE, _("Index %g is out of range for vector "
2506 "with %zu elements."), index, size);
2510 msg (SE, _("%g is not a valid row index for a matrix "
2511 "with dimensions (%zu,%zu)."),
2512 index, size, other_size);
2516 msg (SE, _("%g is not a valid column index for a matrix "
2517 "with dimensions (%zu,%zu)."),
2518 index, other_size, size);
2525 iv->indexes[i] = index - 1;
2531 *iv = (struct index_vector) {
2532 .indexes = xnmalloc (size, sizeof *iv->indexes),
2535 for (size_t i = 0; i < size; i++)
2542 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2544 if (!is_vector (sm))
2546 msg (SE, _("Vector index operator must be applied to vector, "
2547 "not a matrix with dimensions (%zu,%zu)."),
2548 sm->size1, sm->size2);
2556 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2558 if (!matrix_expr_evaluate_vec_all (sm))
2561 gsl_vector sv = to_vector (sm);
2562 struct index_vector iv;
2563 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2566 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2567 sm->size1 == 1 ? iv.n : 1);
2568 gsl_vector dv = to_vector (dm);
2569 for (size_t dx = 0; dx < iv.n; dx++)
2571 size_t sx = iv.indexes[dx];
2572 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2580 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2583 struct index_vector iv0;
2584 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2587 struct index_vector iv1;
2588 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2595 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2596 for (size_t dy = 0; dy < iv0.n; dy++)
2598 size_t sy = iv0.indexes[dy];
2600 for (size_t dx = 0; dx < iv1.n; dx++)
2602 size_t sx = iv1.indexes[dx];
2603 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2611 #define F(NAME, PROTOTYPE) \
2612 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2613 const char *name, gsl_matrix *[], size_t, \
2614 matrix_proto_##PROTOTYPE *);
2619 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2621 if (!is_scalar (subs[index]))
2623 msg (SE, _("Function %s argument %zu must be a scalar, "
2624 "but it has dimensions (%zu,%zu)."),
2625 name, index + 1, subs[index]->size1, subs[index]->size2);
2632 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2634 if (!is_vector (subs[index]))
2636 msg (SE, _("Function %s argument %zu must be a vector, "
2637 "but it has dimensions (%zu,%zu)."),
2638 name, index + 1, subs[index]->size1, subs[index]->size2);
2645 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2647 for (size_t i = 0; i < n_subs; i++)
2649 if (!check_scalar_arg (name, subs, i))
2651 d[i] = to_scalar (subs[i]);
2657 matrix_expr_evaluate_m_d (const char *name,
2658 gsl_matrix *subs[], size_t n_subs,
2659 matrix_proto_m_d *f)
2661 assert (n_subs == 1);
2664 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2668 matrix_expr_evaluate_m_dd (const char *name,
2669 gsl_matrix *subs[], size_t n_subs,
2670 matrix_proto_m_dd *f)
2672 assert (n_subs == 2);
2675 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2679 matrix_expr_evaluate_m_ddd (const char *name,
2680 gsl_matrix *subs[], size_t n_subs,
2681 matrix_proto_m_ddd *f)
2683 assert (n_subs == 3);
2686 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2690 matrix_expr_evaluate_m_m (const char *name UNUSED,
2691 gsl_matrix *subs[], size_t n_subs,
2692 matrix_proto_m_m *f)
2694 assert (n_subs == 1);
2699 matrix_expr_evaluate_m_md (const char *name UNUSED,
2700 gsl_matrix *subs[], size_t n_subs,
2701 matrix_proto_m_md *f)
2703 assert (n_subs == 2);
2704 if (!check_scalar_arg (name, subs, 1))
2706 return f (subs[0], to_scalar (subs[1]));
2710 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2711 gsl_matrix *subs[], size_t n_subs,
2712 matrix_proto_m_mdd *f)
2714 assert (n_subs == 3);
2715 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2717 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2721 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2722 gsl_matrix *subs[], size_t n_subs,
2723 matrix_proto_m_mm *f)
2725 assert (n_subs == 2);
2726 return f (subs[0], subs[1]);
2730 matrix_expr_evaluate_m_v (const char *name,
2731 gsl_matrix *subs[], size_t n_subs,
2732 matrix_proto_m_v *f)
2734 assert (n_subs == 1);
2735 if (!check_vector_arg (name, subs, 0))
2737 gsl_vector v = to_vector (subs[0]);
2742 matrix_expr_evaluate_d_m (const char *name UNUSED,
2743 gsl_matrix *subs[], size_t n_subs,
2744 matrix_proto_d_m *f)
2746 assert (n_subs == 1);
2748 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2749 gsl_matrix_set (m, 0, 0, f (subs[0]));
2754 matrix_expr_evaluate_m_any (const char *name UNUSED,
2755 gsl_matrix *subs[], size_t n_subs,
2756 matrix_proto_m_any *f)
2758 return f (subs, n_subs);
2762 matrix_expr_evaluate_IDENT (const char *name,
2763 gsl_matrix *subs[], size_t n_subs,
2764 matrix_proto_IDENT *f)
2766 assert (n_subs <= 2);
2769 if (!to_scalar_args (name, subs, n_subs, d))
2772 return f (d[0], d[n_subs - 1]);
2776 matrix_expr_evaluate (const struct matrix_expr *e)
2778 if (e->op == MOP_NUMBER)
2780 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2781 gsl_matrix_set (m, 0, 0, e->number);
2784 else if (e->op == MOP_VARIABLE)
2786 const gsl_matrix *src = e->variable->value;
2789 msg (SE, _("Uninitialized variable %s used in expression."),
2794 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2795 gsl_matrix_memcpy (dst, src);
2798 else if (e->op == MOP_EOF)
2800 struct dfm_reader *reader = read_file_open (e->eof);
2801 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2802 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2806 enum { N_LOCAL = 3 };
2807 gsl_matrix *local_subs[N_LOCAL];
2808 gsl_matrix **subs = (e->n_subs < N_LOCAL
2810 : xmalloc (e->n_subs * sizeof *subs));
2812 for (size_t i = 0; i < e->n_subs; i++)
2814 subs[i] = matrix_expr_evaluate (e->subs[i]);
2817 for (size_t j = 0; j < i; j++)
2818 gsl_matrix_free (subs[j]);
2819 if (subs != local_subs)
2825 gsl_matrix *result = NULL;
2828 #define F(NAME, PROTOTYPE) \
2829 case MOP_F_##NAME: \
2830 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2832 matrix_eval_##NAME); \
2838 gsl_matrix_scale (subs[0], -1.0);
2856 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2860 result = matrix_expr_evaluate_not (subs[0]);
2864 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2868 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2872 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2876 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2879 case MOP_PASTE_HORZ:
2880 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2883 case MOP_PASTE_VERT:
2884 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2888 result = gsl_matrix_alloc (0, 0);
2892 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2896 result = matrix_expr_evaluate_vec_all (subs[0]);
2900 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2904 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2908 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
2917 for (size_t i = 0; i < e->n_subs; i++)
2918 if (subs[i] != result)
2919 gsl_matrix_free (subs[i]);
2920 if (subs != local_subs)
2926 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
2929 gsl_matrix *m = matrix_expr_evaluate (e);
2935 msg (SE, _("Expression for %s must evaluate to scalar, "
2936 "not a matrix with dimensions (%zu,%zu)."),
2937 context, m->size1, m->size2);
2938 gsl_matrix_free (m);
2943 gsl_matrix_free (m);
2948 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
2952 if (!matrix_expr_evaluate_scalar (e, context, &d))
2956 if (d < LONG_MIN || d > LONG_MAX)
2958 msg (SE, _("Expression for %s is outside the integer range."), context);
2965 struct matrix_lvalue
2967 struct matrix_var *var;
2968 struct matrix_expr *indexes[2];
2973 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
2976 for (size_t i = 0; i < lvalue->n_indexes; i++)
2977 matrix_expr_destroy (lvalue->indexes[i]);
2980 static struct matrix_lvalue *
2981 matrix_lvalue_parse (struct matrix_state *s)
2983 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
2984 if (!lex_force_id (s->lexer))
2986 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
2987 if (lex_next_token (s->lexer, 1) == T_LPAREN)
2991 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
2995 lex_get_n (s->lexer, 2);
2997 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
2999 lvalue->n_indexes++;
3001 if (lex_match (s->lexer, T_COMMA))
3003 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3005 lvalue->n_indexes++;
3007 if (!lex_force_match (s->lexer, T_RPAREN))
3013 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3019 matrix_lvalue_destroy (lvalue);
3024 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3025 enum index_type index_type, size_t other_size,
3026 struct index_vector *iv)
3031 m = matrix_expr_evaluate (e);
3038 return matrix_normalize_index_vector (m, size, index_type, other_size, iv);
3042 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3043 struct index_vector *iv, gsl_matrix *sm)
3045 gsl_vector dv = to_vector (lvalue->var->value);
3047 /* Convert source matrix 'sm' to source vector 'sv'. */
3048 if (!is_vector (sm))
3050 msg (SE, _("Can't assign matrix with dimensions (%zu,%zu) to subvector."),
3051 sm->size1, sm->size2);
3054 gsl_vector sv = to_vector (sm);
3055 if (iv->n != sv.size)
3057 msg (SE, _("Can't assign vector with %zu elements "
3058 "to subvector with %zu."), sv.size, iv->n);
3062 /* Assign elements. */
3063 for (size_t x = 0; x < iv->n; x++)
3064 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3069 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3070 struct index_vector *iv0,
3071 struct index_vector *iv1, gsl_matrix *sm)
3073 gsl_matrix *dm = lvalue->var->value;
3075 /* Convert source matrix 'sm' to source vector 'sv'. */
3076 if (iv0->n != sm->size1)
3078 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3079 "but source matrix has %zu rows."),
3080 lvalue->var->name, iv0->n, sm->size1);
3083 if (iv1->n != sm->size2)
3085 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3086 "but source matrix has %zu columns."),
3087 lvalue->var->name, iv1->n, sm->size2);
3091 /* Assign elements. */
3092 for (size_t y = 0; y < iv0->n; y++)
3094 size_t f0 = iv0->indexes[y];
3095 for (size_t x = 0; x < iv1->n; x++)
3097 size_t f1 = iv1->indexes[x];
3098 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3104 /* Takes ownership of M. */
3106 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3107 struct index_vector *iv0, struct index_vector *iv1,
3110 if (!lvalue->n_indexes)
3112 gsl_matrix_free (lvalue->var->value);
3113 lvalue->var->value = sm;
3118 bool ok = (lvalue->n_indexes == 1
3119 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3120 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3121 free (iv0->indexes);
3122 free (iv1->indexes);
3123 gsl_matrix_free (sm);
3129 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3130 struct index_vector *iv0,
3131 struct index_vector *iv1)
3133 *iv0 = INDEX_VECTOR_INIT;
3134 *iv1 = INDEX_VECTOR_INIT;
3135 if (!lvalue->n_indexes)
3138 /* Validate destination matrix exists and has the right shape. */
3139 gsl_matrix *dm = lvalue->var->value;
3142 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3145 else if (dm->size1 == 0 || dm->size2 == 0)
3147 msg (SE, _("Cannot index matrix with dimensions (%zu,%zu)."),
3148 dm->size1, dm->size2);
3151 else if (lvalue->n_indexes == 1)
3153 if (!is_vector (dm))
3155 msg (SE, _("Can't use vector indexing on matrix %s with "
3156 "dimensions (%zu,%zu)."),
3157 lvalue->var->name, dm->size1, dm->size2);
3160 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3161 MAX (dm->size1, dm->size2),
3166 assert (lvalue->n_indexes == 2);
3167 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3168 IV_ROW, dm->size2, iv0))
3171 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3172 IV_COLUMN, dm->size1, iv1))
3174 free (iv0->indexes);
3181 /* Takes ownership of M. */
3183 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3185 struct index_vector iv0, iv1;
3186 return (matrix_lvalue_evaluate (lvalue, &iv0, &iv1)
3187 && matrix_lvalue_assign (lvalue, &iv0, &iv1, sm));
3191 /* Matrix command. */
3195 struct matrix_cmd **commands;
3201 enum matrix_cmd_type
3224 struct compute_command
3226 struct matrix_lvalue *lvalue;
3227 struct matrix_expr *rvalue;
3231 struct print_command
3233 struct matrix_expr *expression;
3234 bool use_default_format;
3235 struct fmt_spec format;
3237 int space; /* -1 means NEWPAGE. */
3241 struct string_array literals; /* CLABELS/RLABELS. */
3242 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3248 struct do_if_command
3252 struct matrix_expr *condition;
3253 struct matrix_cmds commands;
3263 struct matrix_var *var;
3264 struct matrix_expr *start;
3265 struct matrix_expr *end;
3266 struct matrix_expr *increment;
3268 /* Loop conditions. */
3269 struct matrix_expr *top_condition;
3270 struct matrix_expr *bottom_condition;
3273 struct matrix_cmds commands;
3279 struct matrix_expr *expression;
3280 struct file_handle *outfile;
3281 struct string_array *variables;
3282 struct matrix_expr *names;
3283 struct stringi_set strings;
3289 struct read_file *rf;
3290 struct matrix_lvalue *dst;
3291 struct matrix_expr *size;
3293 enum fmt_type format;
3301 struct write_command
3303 struct matrix_expr *expression;
3304 struct file_handle *outfile;
3307 enum fmt_type format;
3315 struct display_command
3317 struct matrix_state *state;
3323 struct release_command
3325 struct matrix_var **vars;
3332 struct matrix_lvalue *dst;
3333 struct file_handle *file;
3335 struct string_array variables;
3336 struct matrix_var *names;
3338 /* Treatment of missing values. */
3343 MGET_ERROR, /* Flag error on command. */
3344 MGET_ACCEPT, /* Accept without change, user-missing only. */
3345 MGET_OMIT, /* Drop this case. */
3346 MGET_RECODE /* Recode to 'substitute'. */
3349 double substitute; /* MGET_RECODE only. */
3355 struct msave_command
3357 struct msave_common *common;
3359 struct matrix_expr *expr;
3360 const char *rowtype;
3361 struct matrix_expr *factors;
3362 struct matrix_expr *splits;
3368 struct matrix_state *state;
3369 struct file_handle *file;
3370 struct stringi_set rowtypes;
3374 struct eigen_command
3376 struct matrix_expr *expr;
3377 struct matrix_var *evec;
3378 struct matrix_var *eval;
3382 struct setdiag_command
3384 struct matrix_var *dst;
3385 struct matrix_expr *expr;
3391 struct matrix_expr *expr;
3392 struct matrix_var *u;
3393 struct matrix_var *s;
3394 struct matrix_var *v;
3400 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3401 static bool matrix_cmd_execute (struct matrix_cmd *);
3404 matrix_cmd_destroy (struct matrix_cmd *cmd)
3410 static struct matrix_cmd *
3411 matrix_parse_compute (struct matrix_state *s)
3413 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3414 *cmd = (struct matrix_cmd) {
3415 .type = MCMD_COMPUTE,
3416 .compute = { .lvalue = NULL },
3419 struct compute_command *compute = &cmd->compute;
3420 compute->lvalue = matrix_lvalue_parse (s);
3421 if (!compute->lvalue)
3424 if (!lex_force_match (s->lexer, T_EQUALS))
3427 compute->rvalue = matrix_parse_expr (s);
3428 if (!compute->rvalue)
3434 matrix_cmd_destroy (cmd);
3439 matrix_cmd_execute_compute (struct compute_command *compute)
3441 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3445 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3449 print_labels_destroy (struct print_labels *labels)
3453 string_array_destroy (&labels->literals);
3454 matrix_expr_destroy (labels->expr);
3460 parse_literal_print_labels (struct matrix_state *s,
3461 struct print_labels **labelsp)
3463 lex_match (s->lexer, T_EQUALS);
3464 print_labels_destroy (*labelsp);
3465 *labelsp = xzalloc (sizeof **labelsp);
3466 while (lex_token (s->lexer) != T_SLASH
3467 && lex_token (s->lexer) != T_ENDCMD
3468 && lex_token (s->lexer) != T_STOP)
3470 struct string label = DS_EMPTY_INITIALIZER;
3471 while (lex_token (s->lexer) != T_COMMA
3472 && lex_token (s->lexer) != T_SLASH
3473 && lex_token (s->lexer) != T_ENDCMD
3474 && lex_token (s->lexer) != T_STOP)
3476 if (!ds_is_empty (&label))
3477 ds_put_byte (&label, ' ');
3479 if (lex_token (s->lexer) == T_STRING)
3480 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3483 char *rep = lex_next_representation (s->lexer, 0, 0);
3484 ds_put_cstr (&label, rep);
3489 string_array_append_nocopy (&(*labelsp)->literals,
3490 ds_steal_cstr (&label));
3492 if (!lex_match (s->lexer, T_COMMA))
3498 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3500 lex_match (s->lexer, T_EQUALS);
3501 struct matrix_expr *e = matrix_parse_exp (s);
3505 print_labels_destroy (*labelsp);
3506 *labelsp = xzalloc (sizeof **labelsp);
3507 (*labelsp)->expr = e;
3511 static struct matrix_cmd *
3512 matrix_parse_print (struct matrix_state *s)
3514 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3515 *cmd = (struct matrix_cmd) {
3518 .use_default_format = true,
3522 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3525 for (size_t i = 0; ; i++)
3527 enum token_type t = lex_next_token (s->lexer, i);
3528 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3530 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3532 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3535 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3540 cmd->print.expression = matrix_parse_exp (s);
3541 if (!cmd->print.expression)
3545 while (lex_match (s->lexer, T_SLASH))
3547 if (lex_match_id (s->lexer, "FORMAT"))
3549 lex_match (s->lexer, T_EQUALS);
3550 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3552 cmd->print.use_default_format = false;
3554 else if (lex_match_id (s->lexer, "TITLE"))
3556 lex_match (s->lexer, T_EQUALS);
3557 if (!lex_force_string (s->lexer))
3559 free (cmd->print.title);
3560 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3563 else if (lex_match_id (s->lexer, "SPACE"))
3565 lex_match (s->lexer, T_EQUALS);
3566 if (lex_match_id (s->lexer, "NEWPAGE"))
3567 cmd->print.space = -1;
3568 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3570 cmd->print.space = lex_integer (s->lexer);
3576 else if (lex_match_id (s->lexer, "RLABELS"))
3577 parse_literal_print_labels (s, &cmd->print.rlabels);
3578 else if (lex_match_id (s->lexer, "CLABELS"))
3579 parse_literal_print_labels (s, &cmd->print.clabels);
3580 else if (lex_match_id (s->lexer, "RNAMES"))
3582 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3585 else if (lex_match_id (s->lexer, "CNAMES"))
3587 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3592 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3593 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3601 matrix_cmd_destroy (cmd);
3606 matrix_is_integer (const gsl_matrix *m)
3608 for (size_t y = 0; y < m->size1; y++)
3609 for (size_t x = 0; x < m->size2; x++)
3611 double d = gsl_matrix_get (m, y, x);
3619 matrix_max_magnitude (const gsl_matrix *m)
3622 for (size_t y = 0; y < m->size1; y++)
3623 for (size_t x = 0; x < m->size2; x++)
3625 double d = fabs (gsl_matrix_get (m, y, x));
3633 format_fits (struct fmt_spec format, double x)
3635 char *s = data_out (&(union value) { .f = x }, NULL,
3636 &format, settings_get_fmt_settings ());
3637 bool fits = *s != '*' && !strchr (s, 'E');
3642 static struct fmt_spec
3643 default_format (const gsl_matrix *m, int *log_scale)
3647 double max = matrix_max_magnitude (m);
3649 if (matrix_is_integer (m))
3650 for (int w = 1; w <= 12; w++)
3652 struct fmt_spec format = { .type = FMT_F, .w = w };
3653 if (format_fits (format, -max))
3657 if (max >= 1e9 || max <= 1e-4)
3660 snprintf (s, sizeof s, "%.1e", max);
3662 const char *e = strchr (s, 'e');
3664 *log_scale = atoi (e + 1);
3667 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3671 trimmed_string (double d)
3673 struct substring s = ss_buffer ((char *) &d, sizeof d);
3674 ss_rtrim (&s, ss_cstr (" "));
3675 return ss_xstrdup (s);
3678 static struct string_array *
3679 print_labels_get (const struct print_labels *labels, size_t n,
3680 const char *prefix, bool truncate)
3685 struct string_array *out = xzalloc (sizeof *out);
3686 if (labels->literals.n)
3687 string_array_clone (out, &labels->literals);
3688 else if (labels->expr)
3690 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3691 if (m && is_vector (m))
3693 gsl_vector v = to_vector (m);
3694 for (size_t i = 0; i < v.size; i++)
3695 string_array_append_nocopy (out, trimmed_string (
3696 gsl_vector_get (&v, i)));
3698 gsl_matrix_free (m);
3704 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3707 string_array_append (out, "");
3710 string_array_delete (out, out->n - 1);
3713 for (size_t i = 0; i < out->n; i++)
3715 char *s = out->strings[i];
3716 s[strnlen (s, 8)] = '\0';
3723 matrix_cmd_print_space (int space)
3726 output_item_submit (page_break_item_create ());
3727 for (int i = 0; i < space; i++)
3728 output_log ("%s", "");
3732 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3733 struct fmt_spec format, int log_scale)
3735 matrix_cmd_print_space (print->space);
3737 output_log ("%s", print->title);
3739 output_log (" 10 ** %d X", log_scale);
3741 struct string_array *clabels = print_labels_get (print->clabels,
3742 m->size2, "col", true);
3743 if (clabels && format.w < 8)
3745 struct string_array *rlabels = print_labels_get (print->rlabels,
3746 m->size1, "row", true);
3750 struct string line = DS_EMPTY_INITIALIZER;
3752 ds_put_byte_multiple (&line, ' ', 8);
3753 for (size_t x = 0; x < m->size2; x++)
3754 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3755 output_log_nocopy (ds_steal_cstr (&line));
3758 double scale = pow (10.0, log_scale);
3759 bool numeric = fmt_is_numeric (format.type);
3760 for (size_t y = 0; y < m->size1; y++)
3762 struct string line = DS_EMPTY_INITIALIZER;
3764 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3766 for (size_t x = 0; x < m->size2; x++)
3768 double f = gsl_matrix_get (m, y, x);
3770 ? data_out (&(union value) { .f = f / scale}, NULL,
3771 &format, settings_get_fmt_settings ())
3772 : trimmed_string (f));
3773 ds_put_format (&line, " %s", s);
3776 output_log_nocopy (ds_steal_cstr (&line));
3779 string_array_destroy (rlabels);
3780 string_array_destroy (clabels);
3784 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3785 const struct print_labels *print_labels, size_t n,
3786 const char *name, const char *prefix)
3788 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3790 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3791 for (size_t i = 0; i < n; i++)
3792 pivot_category_create_leaf (
3794 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3795 : pivot_value_new_integer (i + 1)));
3797 d->hide_all_labels = true;
3798 string_array_destroy (labels);
3802 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3803 struct fmt_spec format, int log_scale)
3805 struct pivot_table *table = pivot_table_create__ (
3806 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3809 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3811 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3812 N_("Columns"), "col");
3814 struct pivot_footnote *footnote = NULL;
3817 char *s = xasprintf ("× 10**%d", log_scale);
3818 footnote = pivot_table_create_footnote (
3819 table, pivot_value_new_user_text_nocopy (s));
3822 double scale = pow (10.0, log_scale);
3823 bool numeric = fmt_is_numeric (format.type);
3824 for (size_t y = 0; y < m->size1; y++)
3825 for (size_t x = 0; x < m->size2; x++)
3827 double f = gsl_matrix_get (m, y, x);
3828 struct pivot_value *v;
3831 v = pivot_value_new_number (f / scale);
3832 v->numeric.format = format;
3835 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3837 pivot_value_add_footnote (v, footnote);
3838 pivot_table_put2 (table, y, x, v);
3841 pivot_table_submit (table);
3845 matrix_cmd_execute_print (const struct print_command *print)
3847 if (print->expression)
3849 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3854 struct fmt_spec format = (print->use_default_format
3855 ? default_format (m, &log_scale)
3858 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3859 matrix_cmd_print_text (print, m, format, log_scale);
3861 matrix_cmd_print_tables (print, m, format, log_scale);
3863 gsl_matrix_free (m);
3867 matrix_cmd_print_space (print->space);
3869 output_log ("%s", print->title);
3876 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3877 const char *command_name,
3878 const char *stop1, const char *stop2)
3880 lex_end_of_command (s->lexer);
3881 lex_discard_rest_of_command (s->lexer);
3883 size_t allocated = 0;
3886 while (lex_token (s->lexer) == T_ENDCMD)
3889 if (lex_at_phrase (s->lexer, stop1)
3890 || (stop2 && lex_at_phrase (s->lexer, stop2)))
3893 if (lex_at_phrase (s->lexer, "END MATRIX"))
3895 msg (SE, _("Premature END MATRIX within %s."), command_name);
3899 struct matrix_cmd *cmd = matrix_parse_command (s);
3903 if (c->n >= allocated)
3904 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
3905 c->commands[c->n++] = cmd;
3910 matrix_parse_do_if_clause (struct matrix_state *s,
3911 struct do_if_command *ifc,
3912 bool parse_condition,
3913 size_t *allocated_clauses)
3915 if (ifc->n_clauses >= *allocated_clauses)
3916 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
3917 sizeof *ifc->clauses);
3918 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
3919 *c = (struct do_if_clause) { .condition = NULL };
3921 if (parse_condition)
3923 c->condition = matrix_parse_expr (s);
3928 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
3931 static struct matrix_cmd *
3932 matrix_parse_do_if (struct matrix_state *s)
3934 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3935 *cmd = (struct matrix_cmd) {
3937 .do_if = { .n_clauses = 0 }
3940 size_t allocated_clauses = 0;
3943 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
3946 while (lex_match_phrase (s->lexer, "ELSE IF"));
3948 if (lex_match_id (s->lexer, "ELSE")
3949 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
3952 if (!lex_match_phrase (s->lexer, "END IF"))
3957 matrix_cmd_destroy (cmd);
3962 matrix_cmd_execute_do_if (struct do_if_command *cmd)
3964 for (size_t i = 0; i < cmd->n_clauses; i++)
3966 struct do_if_clause *c = &cmd->clauses[i];
3970 if (!matrix_expr_evaluate_scalar (c->condition,
3971 i ? "ELSE IF" : "DO IF",
3976 for (size_t j = 0; j < c->commands.n; j++)
3977 if (!matrix_cmd_execute (c->commands.commands[j]))
3984 static struct matrix_cmd *
3985 matrix_parse_loop (struct matrix_state *s)
3987 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3988 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
3990 struct loop_command *loop = &cmd->loop;
3991 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
3993 struct substring name = lex_tokss (s->lexer);
3994 loop->var = matrix_var_lookup (s, name);
3996 loop->var = matrix_var_create (s, name);
4001 loop->start = matrix_parse_expr (s);
4002 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4005 loop->end = matrix_parse_expr (s);
4009 if (lex_match (s->lexer, T_BY))
4011 loop->increment = matrix_parse_expr (s);
4012 if (!loop->increment)
4017 if (lex_match_id (s->lexer, "IF"))
4019 loop->top_condition = matrix_parse_expr (s);
4020 if (!loop->top_condition)
4024 bool was_in_loop = s->in_loop;
4026 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4028 s->in_loop = was_in_loop;
4032 if (!lex_match_phrase (s->lexer, "END LOOP"))
4035 if (lex_match_id (s->lexer, "IF"))
4037 loop->bottom_condition = matrix_parse_expr (s);
4038 if (!loop->bottom_condition)
4045 matrix_cmd_destroy (cmd);
4050 matrix_cmd_execute_loop (struct loop_command *cmd)
4054 long int increment = 1;
4057 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4058 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4060 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4064 if (increment > 0 ? value > end
4065 : increment < 0 ? value < end
4070 int mxloops = settings_get_mxloops ();
4071 for (int i = 0; i < mxloops; i++)
4076 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4078 gsl_matrix_free (cmd->var->value);
4079 cmd->var->value = NULL;
4081 if (!cmd->var->value)
4082 cmd->var->value = gsl_matrix_alloc (1, 1);
4084 gsl_matrix_set (cmd->var->value, 0, 0, value);
4087 if (cmd->top_condition)
4090 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4095 for (size_t j = 0; j < cmd->commands.n; j++)
4096 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4099 if (cmd->bottom_condition)
4102 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4103 "END LOOP IF", &d) || d > 0)
4110 if (increment > 0 ? value > end : value < end)
4116 static struct matrix_cmd *
4117 matrix_parse_break (struct matrix_state *s)
4121 msg (SE, _("BREAK not inside LOOP."));
4125 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4126 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4130 static struct matrix_cmd *
4131 matrix_parse_display (struct matrix_state *s)
4133 bool dictionary = false;
4134 bool status = false;
4135 while (lex_token (s->lexer) == T_ID)
4137 if (lex_match_id (s->lexer, "DICTIONARY"))
4139 else if (lex_match_id (s->lexer, "STATUS"))
4143 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4147 if (!dictionary && !status)
4148 dictionary = status = true;
4150 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4151 *cmd = (struct matrix_cmd) {
4152 .type = MCMD_DISPLAY,
4155 .dictionary = dictionary,
4163 compare_matrix_var_pointers (const void *a_, const void *b_)
4165 const struct matrix_var *const *ap = a_;
4166 const struct matrix_var *const *bp = b_;
4167 const struct matrix_var *a = *ap;
4168 const struct matrix_var *b = *bp;
4169 return strcmp (a->name, b->name);
4173 matrix_cmd_execute_display (struct display_command *cmd)
4175 const struct matrix_state *s = cmd->state;
4177 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4178 pivot_dimension_create (
4179 table, PIVOT_AXIS_COLUMN, N_("Property"),
4180 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4182 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4185 const struct matrix_var *var;
4186 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4188 vars[n_vars++] = var;
4189 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4191 struct pivot_dimension *rows = pivot_dimension_create (
4192 table, PIVOT_AXIS_ROW, N_("Variable"));
4193 for (size_t i = 0; i < n_vars; i++)
4195 const struct matrix_var *var = vars[i];
4196 pivot_category_create_leaf (
4197 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4199 size_t r = var->value->size1;
4200 size_t c = var->value->size2;
4201 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4202 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4203 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4205 pivot_table_submit (table);
4208 static struct matrix_cmd *
4209 matrix_parse_release (struct matrix_state *s)
4211 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4212 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4214 size_t allocated_vars = 0;
4215 while (lex_token (s->lexer) == T_ID)
4217 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4220 if (cmd->release.n_vars >= allocated_vars)
4221 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4222 sizeof *cmd->release.vars);
4223 cmd->release.vars[cmd->release.n_vars++] = var;
4226 lex_error (s->lexer, _("Variable name expected."));
4229 if (!lex_match (s->lexer, T_COMMA))
4237 matrix_cmd_execute_release (struct release_command *cmd)
4239 for (size_t i = 0; i < cmd->n_vars; i++)
4241 struct matrix_var *var = cmd->vars[i];
4242 gsl_matrix_free (var->value);
4247 static struct matrix_cmd *
4248 matrix_parse_save (struct matrix_state *s)
4250 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4251 *cmd = (struct matrix_cmd) {
4254 .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4258 struct save_command *save = &cmd->save;
4259 save->expression = matrix_parse_exp (s);
4260 if (!save->expression)
4263 while (lex_match (s->lexer, T_SLASH))
4265 if (lex_match_id (s->lexer, "OUTFILE"))
4267 lex_match (s->lexer, T_EQUALS);
4269 fh_unref (save->outfile);
4270 save->outfile = (lex_match (s->lexer, T_ASTERISK)
4272 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4276 else if (lex_match_id (s->lexer, "VARIABLES"))
4278 lex_match (s->lexer, T_EQUALS);
4282 struct dictionary *d = dict_create (get_default_encoding ());
4283 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4284 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4289 string_array_destroy (save->variables);
4290 if (!save->variables)
4291 save->variables = xmalloc (sizeof *save->variables);
4292 *save->variables = (struct string_array) {
4298 else if (lex_match_id (s->lexer, "NAMES"))
4300 lex_match (s->lexer, T_EQUALS);
4301 matrix_expr_destroy (save->names);
4302 save->names = matrix_parse_exp (s);
4306 else if (lex_match_id (s->lexer, "STRINGS"))
4308 lex_match (s->lexer, T_EQUALS);
4309 while (lex_token (s->lexer) == T_ID)
4311 stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4313 if (!lex_match (s->lexer, T_COMMA))
4319 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4327 if (s->prev_save_outfile)
4328 save->outfile = fh_ref (s->prev_save_outfile);
4331 lex_sbc_missing ("OUTFILE");
4335 fh_unref (s->prev_save_outfile);
4336 s->prev_save_outfile = fh_ref (save->outfile);
4338 if (save->variables && save->names)
4340 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4341 matrix_expr_destroy (save->names);
4348 matrix_cmd_destroy (cmd);
4353 matrix_cmd_execute_save (const struct save_command *save)
4355 assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4357 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4362 struct dictionary *dict = dict_create (get_default_encoding ());
4363 struct string_array names = { .n = 0 };
4364 if (save->variables)
4365 string_array_clone (&names, save->variables);
4366 else if (save->names)
4368 gsl_matrix *nm = matrix_expr_evaluate (save->names);
4369 if (nm && is_vector (nm))
4371 gsl_vector nv = to_vector (nm);
4372 for (size_t i = 0; i < nv.size; i++)
4374 char *name = trimmed_string (gsl_vector_get (&nv, i));
4375 if (dict_id_is_valid (dict, name, true))
4376 string_array_append_nocopy (&names, name);
4383 struct stringi_set strings;
4384 stringi_set_clone (&strings, &save->strings);
4386 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4391 name = names.strings[i];
4394 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4398 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4399 struct variable *var = dict_create_var (dict, name, width);
4402 msg (SE, _("Duplicate variable name %s in SAVE statement."),
4408 if (!stringi_set_is_empty (&strings))
4410 const char *example = stringi_set_node_get_string (
4411 stringi_set_first (&strings));
4412 msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4413 "with that name was found on SAVE.",
4414 "STRINGS specified %2$zu variables, including %1$s, "
4415 "whose names were not found on SAVE.",
4416 stringi_set_count (&strings)),
4417 example, stringi_set_count (&strings));
4420 stringi_set_destroy (&strings);
4421 string_array_destroy (&names);
4425 gsl_matrix_free (m);
4430 struct casewriter *writer = any_writer_open (save->outfile, dict);
4433 gsl_matrix_free (m);
4438 const struct caseproto *proto = dict_get_proto (dict);
4439 for (size_t y = 0; y < m->size1; y++)
4441 struct ccase *c = case_create (proto);
4442 for (size_t x = 0; x < m->size2; x++)
4444 double d = gsl_matrix_get (m, y, x);
4445 int width = caseproto_get_width (proto, x);
4446 union value *value = case_data_rw_idx (c, x);
4450 memcpy (value->s, &d, width);
4452 casewriter_write (writer, c);
4454 casewriter_destroy (writer);
4456 gsl_matrix_free (m);
4460 static struct matrix_cmd *
4461 matrix_parse_read (struct matrix_state *s)
4463 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4464 *cmd = (struct matrix_cmd) {
4466 .read = { .format = FMT_F },
4469 struct file_handle *fh = NULL;
4470 char *encoding = NULL;
4471 struct read_command *read = &cmd->read;
4472 read->dst = matrix_lvalue_parse (s);
4477 int repetitions = 0;
4478 int record_width = 0;
4479 bool seen_format = false;
4480 while (lex_match (s->lexer, T_SLASH))
4482 if (lex_match_id (s->lexer, "FILE"))
4484 lex_match (s->lexer, T_EQUALS);
4487 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4491 else if (lex_match_id (s->lexer, "ENCODING"))
4493 lex_match (s->lexer, T_EQUALS);
4494 if (!lex_force_string (s->lexer))
4498 encoding = ss_xstrdup (lex_tokss (s->lexer));
4502 else if (lex_match_id (s->lexer, "FIELD"))
4504 lex_match (s->lexer, T_EQUALS);
4506 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4508 read->c1 = lex_integer (s->lexer);
4510 if (!lex_force_match (s->lexer, T_TO)
4511 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4513 read->c2 = lex_integer (s->lexer) + 1;
4516 record_width = read->c2 - read->c1;
4517 if (lex_match (s->lexer, T_BY))
4519 if (!lex_force_int_range (s->lexer, "BY", 1,
4520 read->c2 - read->c1))
4522 by = lex_integer (s->lexer);
4525 if (record_width % by)
4527 msg (SE, _("BY %d does not evenly divide record width %d."),
4535 else if (lex_match_id (s->lexer, "SIZE"))
4537 lex_match (s->lexer, T_EQUALS);
4538 read->size = matrix_parse_exp (s);
4542 else if (lex_match_id (s->lexer, "MODE"))
4544 lex_match (s->lexer, T_EQUALS);
4545 if (lex_match_id (s->lexer, "RECTANGULAR"))
4546 read->symmetric = false;
4547 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4548 read->symmetric = true;
4551 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4555 else if (lex_match_id (s->lexer, "REREAD"))
4556 read->reread = true;
4557 else if (lex_match_id (s->lexer, "FORMAT"))
4561 lex_sbc_only_once ("FORMAT");
4566 lex_match (s->lexer, T_EQUALS);
4568 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4571 const char *p = lex_tokcstr (s->lexer);
4572 if (c_isdigit (p[0]))
4574 repetitions = atoi (p);
4575 p += strspn (p, "0123456789");
4576 if (!fmt_from_name (p, &read->format))
4578 lex_error (s->lexer, _("Unknown format %s."), p);
4583 else if (!fmt_from_name (p, &read->format))
4585 struct fmt_spec format;
4586 if (!parse_format_specifier (s->lexer, &format))
4588 read->format = format.type;
4595 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4596 "REREAD", "FORMAT");
4603 lex_sbc_missing ("FIELD");
4607 if (!read->dst->n_indexes && !read->size)
4609 msg (SE, _("SIZE is required for reading data into a full matrix "
4610 "(as opposed to a submatrix)."));
4616 if (s->prev_read_file)
4617 fh = fh_ref (s->prev_read_file);
4620 lex_sbc_missing ("FILE");
4624 fh_unref (s->prev_read_file);
4625 s->prev_read_file = fh_ref (fh);
4627 read->rf = read_file_create (s, fh);
4630 free (read->rf->encoding);
4631 read->rf->encoding = encoding;
4635 /* Field width may be specified in multiple ways:
4638 2. The format on FORMAT.
4639 3. The repetition factor on FORMAT.
4641 (2) and (3) are mutually exclusive.
4643 If more than one of these is present, they must agree. If none of them is
4644 present, then free-field format is used.
4646 if (repetitions > record_width)
4648 msg (SE, _("%d repetitions cannot fit in record width %d."),
4649 repetitions, record_width);
4652 int w = (repetitions ? record_width / repetitions
4658 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4659 "which implies field width %d, "
4660 "but BY specifies field width %d."),
4661 repetitions, record_width, w, by);
4663 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4671 matrix_cmd_destroy (cmd);
4677 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4678 gsl_matrix *m, struct substring p, size_t y, size_t x)
4680 const char *input_encoding = dfm_reader_get_encoding (reader);
4682 char *error = data_in (p, input_encoding, read->format,
4683 settings_get_fmt_settings (), &v, 0, NULL);
4684 /* XXX report error if value is missing */
4686 msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4689 gsl_matrix_set (m, y, x, v.f);
4690 if (read->symmetric && x != y)
4691 gsl_matrix_set (m, x, y, v.f);
4696 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4697 struct substring *line)
4699 if (dfm_eof (reader))
4701 msg (SE, _("Unexpected end of file reading matrix data."));
4704 dfm_expand_tabs (reader);
4705 *line = ss_substr (dfm_get_record (reader),
4706 read->c1 - 1, read->c2 - read->c1);
4711 matrix_read (struct read_command *read, struct dfm_reader *reader,
4714 for (size_t y = 0; y < m->size1; y++)
4716 size_t nx = read->symmetric ? y + 1 : m->size2;
4718 struct substring line = ss_empty ();
4719 for (size_t x = 0; x < nx; x++)
4726 ss_ltrim (&line, ss_cstr (" ,"));
4727 if (!ss_is_empty (line))
4729 if (!matrix_read_line (read, reader, &line))
4731 dfm_forward_record (reader);
4734 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4738 if (!matrix_read_line (read, reader, &line))
4740 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4741 int f = x % fields_per_line;
4742 if (f == fields_per_line - 1)
4743 dfm_forward_record (reader);
4745 p = ss_substr (line, read->w * f, read->w);
4748 matrix_read_set_field (read, reader, m, p, y, x);
4752 dfm_forward_record (reader);
4755 ss_ltrim (&line, ss_cstr (" ,"));
4756 if (!ss_is_empty (line))
4757 msg (SW, _("Trailing garbage on line \"%.*s\""),
4758 (int) line.length, line.string);
4764 matrix_cmd_execute_read (struct read_command *read)
4766 struct index_vector iv0, iv1;
4767 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4770 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4773 gsl_matrix *m = matrix_expr_evaluate (read->size);
4779 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4780 gsl_matrix_free (m);
4786 gsl_vector v = to_vector (m);
4790 d[0] = gsl_vector_get (&v, 0);
4793 else if (v.size == 2)
4795 d[0] = gsl_vector_get (&v, 0);
4796 d[1] = gsl_vector_get (&v, 1);
4800 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4801 gsl_matrix_free (m);
4806 gsl_matrix_free (m);
4808 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4810 msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4820 if (read->dst->n_indexes)
4822 size_t submatrix_size[2];
4823 if (read->dst->n_indexes == 2)
4825 submatrix_size[0] = iv0.n;
4826 submatrix_size[1] = iv1.n;
4828 else if (read->dst->var->value->size1 == 1)
4830 submatrix_size[0] = 1;
4831 submatrix_size[1] = iv0.n;
4835 submatrix_size[0] = iv0.n;
4836 submatrix_size[1] = 1;
4841 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4843 msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4846 submatrix_size[0], submatrix_size[1]);
4854 size[0] = submatrix_size[0];
4855 size[1] = submatrix_size[1];
4859 struct dfm_reader *reader = read_file_open (read->rf);
4861 dfm_reread_record (reader, 1);
4863 if (read->symmetric && size[0] != size[1])
4865 msg (SE, _("Cannot read matrix with non-square dimensions (%zu,%zu) "
4866 "using READ with MODE=SYMMETRIC."),
4872 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4873 matrix_read (read, reader, tmp);
4874 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4877 static struct matrix_cmd *
4878 matrix_parse_write (struct matrix_state *s)
4880 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4881 *cmd = (struct matrix_cmd) {
4883 .write = { .format = FMT_F },
4886 struct write_command *write = &cmd->write;
4887 write->expression = matrix_parse_exp (s);
4888 if (!write->expression)
4892 int repetitions = 0;
4893 int record_width = 0;
4894 bool seen_format = false;
4895 while (lex_match (s->lexer, T_SLASH))
4897 if (lex_match_id (s->lexer, "OUTFILE"))
4899 lex_match (s->lexer, T_EQUALS);
4901 fh_unref (write->outfile);
4902 write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
4903 if (!write->outfile)
4906 else if (lex_match_id (s->lexer, "ENCODING"))
4908 lex_match (s->lexer, T_EQUALS);
4909 if (!lex_force_string (s->lexer))
4912 free (write->encoding);
4913 write->encoding = ss_xstrdup (lex_tokss (s->lexer));
4917 else if (lex_match_id (s->lexer, "FIELD"))
4919 lex_match (s->lexer, T_EQUALS);
4921 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4923 write->c1 = lex_integer (s->lexer);
4925 if (!lex_force_match (s->lexer, T_TO)
4926 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
4928 write->c2 = lex_integer (s->lexer) + 1;
4931 record_width = write->c2 - write->c1;
4932 if (lex_match (s->lexer, T_BY))
4934 if (!lex_force_int_range (s->lexer, "BY", 1,
4935 write->c2 - write->c1))
4937 by = lex_integer (s->lexer);
4940 if (record_width % by)
4942 msg (SE, _("BY %d does not evenly divide record width %d."),
4950 else if (lex_match_id (s->lexer, "MODE"))
4952 lex_match (s->lexer, T_EQUALS);
4953 if (lex_match_id (s->lexer, "RECTANGULAR"))
4954 write->triangular = false;
4955 else if (lex_match_id (s->lexer, "TRIANGULAR"))
4956 write->triangular = true;
4959 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
4963 else if (lex_match_id (s->lexer, "HOLD"))
4965 else if (lex_match_id (s->lexer, "FORMAT"))
4969 lex_sbc_only_once ("FORMAT");
4974 lex_match (s->lexer, T_EQUALS);
4976 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4979 const char *p = lex_tokcstr (s->lexer);
4980 if (c_isdigit (p[0]))
4982 repetitions = atoi (p);
4983 p += strspn (p, "0123456789");
4984 if (!fmt_from_name (p, &write->format))
4986 lex_error (s->lexer, _("Unknown format %s."), p);
4991 else if (!fmt_from_name (p, &write->format))
4993 struct fmt_spec format;
4994 if (!parse_format_specifier (s->lexer, &format))
4996 write->format = format.type;
4997 write->w = format.w;
4998 write->d = format.d;
5003 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5011 lex_sbc_missing ("FIELD");
5015 if (!write->outfile)
5017 if (s->prev_write_outfile)
5018 write->outfile = fh_ref (s->prev_write_outfile);
5021 lex_sbc_missing ("OUTFILE");
5025 fh_unref (s->prev_write_outfile);
5026 s->prev_write_outfile = fh_ref (write->outfile);
5028 /* Field width may be specified in multiple ways:
5031 2. The format on FORMAT.
5032 3. The repetition factor on FORMAT.
5034 (2) and (3) are mutually exclusive.
5036 If more than one of these is present, they must agree. If none of them is
5037 present, then free-field format is used.
5039 if (repetitions > record_width)
5041 msg (SE, _("%d repetitions cannot fit in record width %d."),
5042 repetitions, record_width);
5045 int w = (repetitions ? record_width / repetitions
5046 : write->w ? write->w
5051 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5052 "which implies field width %d, "
5053 "but BY specifies field width %d."),
5054 repetitions, record_width, w, by);
5056 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5064 matrix_cmd_destroy (cmd);
5069 matrix_cmd_execute_write (struct write_command *write)
5071 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5075 if (write->triangular && m->size1 != m->size2)
5077 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5078 "the matrix to be written has dimensions (%zu,%zu)."),
5079 m->size1, m->size2);
5080 gsl_matrix_free (m);
5084 struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
5087 gsl_matrix_free (m);
5091 const struct fmt_settings *settings = settings_get_fmt_settings ();
5092 struct fmt_spec format = {
5093 .type = write->format,
5094 .w = write->w ? write->w : 40,
5097 struct u8_line line = U8_LINE_INITIALIZER;
5098 for (size_t y = 0; y < m->size1; y++)
5100 size_t nx = write->triangular ? y + 1 : m->size2;
5102 for (size_t x = 0; x < nx; x++)
5104 /* XXX string values */
5105 union value v = { .f = gsl_matrix_get (m, y, x) };
5107 ? data_out (&v, NULL, &format, settings)
5108 : data_out_stretchy (&v, NULL, &format, settings, NULL));
5109 size_t len = strlen (s);
5110 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5111 if (width + x0 > write->c2)
5113 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5114 u8_line_clear (&line);
5117 u8_line_put (&line, x0, x0 + width, s, len);
5120 x0 += write->w ? write->w : width + 1;
5123 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5124 u8_line_clear (&line);
5126 u8_line_destroy (&line);
5127 dfm_close_writer (writer);
5129 gsl_matrix_free (m);
5132 static struct matrix_cmd *
5133 matrix_parse_get (struct matrix_state *s)
5135 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5136 *cmd = (struct matrix_cmd) {
5139 .user = { .treatment = MGET_ERROR },
5140 .system = { .treatment = MGET_ERROR },
5144 struct get_command *get = &cmd->get;
5145 get->dst = matrix_lvalue_parse (s);
5149 while (lex_match (s->lexer, T_SLASH))
5151 if (lex_match_id (s->lexer, "FILE"))
5153 if (get->variables.n)
5155 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5158 lex_match (s->lexer, T_EQUALS);
5160 fh_unref (get->file);
5161 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5165 else if (lex_match_id (s->lexer, "ENCODING"))
5167 if (get->variables.n)
5169 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5172 lex_match (s->lexer, T_EQUALS);
5173 if (!lex_force_string (s->lexer))
5176 free (get->encoding);
5177 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5181 else if (lex_match_id (s->lexer, "VARIABLES"))
5183 lex_match (s->lexer, T_EQUALS);
5185 struct dictionary *dict = NULL;
5188 dict = dataset_dict (s->dataset);
5189 if (dict_get_var_cnt (dict) == 0)
5191 lex_error (s->lexer, _("GET cannot read empty active file."));
5197 struct casereader *reader = any_reader_open_and_decode (
5198 get->file, get->encoding, &dict, NULL);
5201 casereader_destroy (reader);
5204 struct variable **vars;
5206 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5207 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5214 string_array_clear (&get->variables);
5215 for (size_t i = 0; i < n_vars; i++)
5216 string_array_append (&get->variables, var_get_name (vars[i]));
5220 else if (lex_match_id (s->lexer, "NAMES"))
5222 lex_match (s->lexer, T_EQUALS);
5223 if (!lex_force_id (s->lexer))
5226 struct substring name = lex_tokss (s->lexer);
5227 get->names = matrix_var_lookup (s, name);
5229 get->names = matrix_var_create (s, name);
5232 else if (lex_match_id (s->lexer, "MISSING"))
5234 lex_match (s->lexer, T_EQUALS);
5235 if (lex_match_id (s->lexer, "ACCEPT"))
5236 get->user.treatment = MGET_ACCEPT;
5237 else if (lex_match_id (s->lexer, "OMIT"))
5238 get->user.treatment = MGET_OMIT;
5239 else if (lex_is_number (s->lexer))
5241 get->user.treatment = MGET_RECODE;
5242 get->user.substitute = lex_number (s->lexer);
5247 lex_error (s->lexer, NULL);
5251 else if (lex_match_id (s->lexer, "SYSMIS"))
5253 lex_match (s->lexer, T_EQUALS);
5254 if (lex_match_id (s->lexer, "OMIT"))
5255 get->user.treatment = MGET_OMIT;
5256 else if (lex_is_number (s->lexer))
5258 get->user.treatment = MGET_RECODE;
5259 get->user.substitute = lex_number (s->lexer);
5264 lex_error (s->lexer, NULL);
5270 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5271 "MISSING", "SYSMIS");
5278 matrix_cmd_destroy (cmd);
5283 matrix_cmd_execute_get (struct get_command *get)
5285 assert (get->file); /* XXX */
5287 struct dictionary *dict;
5288 struct casereader *reader = any_reader_open_and_decode (
5289 get->file, get->encoding, &dict, NULL);
5293 const struct variable **vars = xnmalloc (
5294 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5298 if (get->variables.n)
5300 for (size_t i = 0; i < get->variables.n; i++)
5302 const char *name = get->variables.strings[i];
5303 const struct variable *var = dict_lookup_var (dict, name);
5306 msg (SE, _("GET: Data file does not contain variable %s."),
5312 if (!var_is_numeric (var))
5314 msg (SE, _("GET: Variable %s is not numeric."), name);
5319 vars[n_vars++] = var;
5324 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5326 const struct variable *var = dict_get_var (dict, i);
5327 if (!var_is_numeric (var))
5329 msg (SE, _("GET: Variable %s is not numeric."),
5330 var_get_name (var));
5335 vars[n_vars++] = var;
5340 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5341 long long int casenum = 1;
5343 for (struct ccase *c = casereader_read (reader); c;
5344 c = casereader_read (reader), casenum++)
5346 if (n_rows >= m->size1)
5348 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5349 for (size_t y = 0; y < n_rows; y++)
5350 for (size_t x = 0; x < n_vars; x++)
5351 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5352 gsl_matrix_free (m);
5357 for (size_t x = 0; x < n_vars; x++)
5359 const struct variable *var = vars[x];
5360 double d = case_num (c, var);
5363 if (get->system.treatment == MGET_RECODE)
5364 d = get->system.substitute;
5365 else if (get->system.treatment == MGET_OMIT)
5369 msg (SE, _("GET: Variable %s in case %lld "
5370 "is system-missing."),
5371 var_get_name (var), casenum);
5375 else if (var_is_num_missing (var, d, MV_USER))
5377 if (get->user.treatment == MGET_RECODE)
5378 d = get->user.substitute;
5379 else if (get->user.treatment == MGET_OMIT)
5381 else if (get->user.treatment != MGET_ACCEPT)
5383 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5385 var_get_name (var), casenum, d);
5389 gsl_matrix_set (m, n_rows, x, d);
5397 casereader_destroy (reader);
5401 matrix_lvalue_evaluate_and_assign (get->dst, m);
5404 gsl_matrix_free (m);
5410 match_rowtype (struct lexer *lexer)
5412 static const char *rowtypes[] = {
5413 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5415 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5417 for (size_t i = 0; i < n_rowtypes; i++)
5418 if (lex_match_id (lexer, rowtypes[i]))
5421 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5426 parse_var_names (struct lexer *lexer, struct string_array *sa)
5428 lex_match (lexer, T_EQUALS);
5430 string_array_clear (sa);
5432 struct dictionary *dict = dict_create (get_default_encoding ());
5433 dict_create_var_assert (dict, "ROWTYPE_", 8);
5434 dict_create_var_assert (dict, "VARNAME_", 8);
5437 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5438 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5443 string_array_clear (sa);
5444 sa->strings = names;
5445 sa->n = sa->allocated = n_names;
5451 msave_common_uninit (struct msave_common *common)
5455 fh_unref (common->outfile);
5456 string_array_destroy (&common->variables);
5457 string_array_destroy (&common->fnames);
5458 string_array_destroy (&common->snames);
5463 compare_variables (const char *keyword,
5464 const struct string_array *new,
5465 const struct string_array *old)
5471 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5472 "on the first MSAVE within MATRIX."), keyword);
5475 else if (!string_array_equal_case (old, new))
5477 msg (SE, _("%s must specify the same variables each time within "
5478 "a given MATRIX."), keyword);
5485 static struct matrix_cmd *
5486 matrix_parse_msave (struct matrix_state *s)
5488 struct msave_common common = { .outfile = NULL };
5489 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5490 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5492 struct msave_command *msave = &cmd->msave;
5493 if (lex_token (s->lexer) == T_ID)
5494 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5495 msave->expr = matrix_parse_exp (s);
5499 while (lex_match (s->lexer, T_SLASH))
5501 if (lex_match_id (s->lexer, "TYPE"))
5503 lex_match (s->lexer, T_EQUALS);
5505 msave->rowtype = match_rowtype (s->lexer);
5506 if (!msave->rowtype)
5509 else if (lex_match_id (s->lexer, "OUTFILE"))
5511 lex_match (s->lexer, T_EQUALS);
5513 fh_unref (common.outfile);
5514 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5515 if (!common.outfile)
5518 else if (lex_match_id (s->lexer, "VARIABLES"))
5520 if (!parse_var_names (s->lexer, &common.variables))
5523 else if (lex_match_id (s->lexer, "FNAMES"))
5525 if (!parse_var_names (s->lexer, &common.fnames))
5528 else if (lex_match_id (s->lexer, "SNAMES"))
5530 if (!parse_var_names (s->lexer, &common.snames))
5533 else if (lex_match_id (s->lexer, "SPLIT"))
5535 lex_match (s->lexer, T_EQUALS);
5537 matrix_expr_destroy (msave->splits);
5538 msave->splits = matrix_parse_exp (s);
5542 else if (lex_match_id (s->lexer, "FACTOR"))
5544 lex_match (s->lexer, T_EQUALS);
5546 matrix_expr_destroy (msave->factors);
5547 msave->factors = matrix_parse_exp (s);
5548 if (!msave->factors)
5553 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5554 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5558 if (!msave->rowtype)
5560 lex_sbc_missing ("TYPE");
5563 common.has_splits = msave->splits || common.snames.n;
5564 common.has_factors = msave->factors || common.fnames.n;
5566 struct msave_common *c = s->common ? s->common : &common;
5567 if (c->fnames.n && !msave->factors)
5569 msg (SE, _("FNAMES requires FACTOR."));
5572 if (c->snames.n && !msave->splits)
5574 msg (SE, _("SNAMES requires SPLIT."));
5577 if (c->has_factors && !common.has_factors)
5579 msg (SE, _("%s is required because it was present on the first "
5580 "MSAVE in this MATRIX command."), "FACTOR");
5583 if (c->has_splits && !common.has_splits)
5585 msg (SE, _("%s is required because it was present on the first "
5586 "MSAVE in this MATRIX command."), "SPLIT");
5592 if (!common.outfile)
5594 lex_sbc_missing ("OUTFILE");
5597 s->common = xmemdup (&common, sizeof common);
5603 bool same = common.outfile == s->common->outfile;
5604 fh_unref (common.outfile);
5607 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5608 "within a single MATRIX command."));
5612 if (!compare_variables ("VARIABLES",
5613 &common.variables, &s->common->variables)
5614 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5615 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5617 msave_common_uninit (&common);
5619 msave->common = s->common;
5620 if (!msave->varname_)
5621 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5625 msave_common_uninit (&common);
5626 matrix_cmd_destroy (cmd);
5631 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5633 gsl_matrix *m = matrix_expr_evaluate (e);
5639 msg (SE, _("%s expression must evaluate to vector, not a matrix with "
5640 "dimensions (%zu,%zu)."),
5641 name, m->size1, m->size2);
5642 gsl_matrix_free (m);
5646 return matrix_to_vector (m);
5650 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5652 for (size_t i = 0; i < vars->n; i++)
5653 if (!dict_create_var (d, vars->strings[i], 0))
5654 return vars->strings[i];
5658 static struct dictionary *
5659 msave_create_dict (const struct msave_common *common)
5661 struct dictionary *dict = dict_create (get_default_encoding ());
5663 const char *dup_split = msave_add_vars (dict, &common->snames);
5666 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5670 dict_create_var_assert (dict, "ROWTYPE_", 8);
5672 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5675 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5679 dict_create_var_assert (dict, "VARNAME_", 8);
5681 const char *dup_var = msave_add_vars (dict, &common->variables);
5684 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5696 matrix_cmd_execute_msave (struct msave_command *msave)
5698 struct msave_common *common = msave->common;
5699 gsl_matrix *m = NULL;
5700 gsl_vector *factors = NULL;
5701 gsl_vector *splits = NULL;
5703 m = matrix_expr_evaluate (msave->expr);
5707 if (!common->variables.n)
5708 for (size_t i = 0; i < m->size2; i++)
5709 string_array_append_nocopy (&common->variables,
5710 xasprintf ("COL%zu", i + 1));
5712 if (m->size2 != common->variables.n)
5714 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5715 m->size2, common->variables.n);
5721 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5725 if (!common->fnames.n)
5726 for (size_t i = 0; i < factors->size; i++)
5727 string_array_append_nocopy (&common->fnames,
5728 xasprintf ("FAC%zu", i + 1));
5730 if (factors->size != common->fnames.n)
5732 msg (SE, _("There are %zu factor variables, "
5733 "but %zu split values were supplied."),
5734 common->fnames.n, factors->size);
5741 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5745 if (!common->fnames.n)
5746 for (size_t i = 0; i < splits->size; i++)
5747 string_array_append_nocopy (&common->fnames,
5748 xasprintf ("SPL%zu", i + 1));
5750 if (splits->size != common->fnames.n)
5752 msg (SE, _("There are %zu split variables, "
5753 "but %zu split values were supplied."),
5754 common->fnames.n, splits->size);
5759 if (!common->writer)
5761 struct dictionary *dict = msave_create_dict (common);
5765 common->writer = any_writer_open (common->outfile, dict);
5766 if (!common->writer)
5772 common->dict = dict;
5775 for (size_t y = 0; y < m->size1; y++)
5777 struct ccase *c = case_create (dict_get_proto (common->dict));
5780 /* Split variables */
5782 for (size_t i = 0; i < splits->size; i++)
5783 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5786 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5787 msave->rowtype, ' ');
5791 for (size_t i = 0; i < factors->size; i++)
5792 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5795 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5796 msave->varname_, ' ');
5798 /* Continuous variables. */
5799 for (size_t x = 0; x < m->size2; x++)
5800 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5801 casewriter_write (common->writer, c);
5805 gsl_matrix_free (m);
5806 gsl_vector_free (factors);
5807 gsl_vector_free (splits);
5810 static struct matrix_cmd *
5811 matrix_parse_mget (struct matrix_state *s)
5813 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5814 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5816 struct mget_command *mget = &cmd->mget;
5820 if (lex_match_id (s->lexer, "FILE"))
5822 lex_match (s->lexer, T_EQUALS);
5824 fh_unref (mget->file);
5825 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5829 else if (lex_match_id (s->lexer, "TYPE"))
5831 lex_match (s->lexer, T_EQUALS);
5832 while (lex_token (s->lexer) != T_SLASH
5833 && lex_token (s->lexer) != T_ENDCMD)
5835 const char *rowtype = match_rowtype (s->lexer);
5839 stringi_set_insert (&mget->rowtypes, rowtype);
5844 lex_error_expecting (s->lexer, "FILE", "TYPE");
5847 if (lex_token (s->lexer) == T_ENDCMD)
5850 if (!lex_force_match (s->lexer, T_SLASH))
5856 matrix_cmd_destroy (cmd);
5860 static const struct variable *
5861 get_a8_var (const struct dictionary *d, const char *name)
5863 const struct variable *v = dict_lookup_var (d, name);
5866 msg (SE, _("Matrix data file lacks %s variable."), name);
5869 if (var_get_width (v) != 8)
5871 msg (SE, _("%s variable in matrix data file must be "
5872 "8-byte string, but it has width %d."),
5873 name, var_get_width (v));
5880 is_rowtype (const union value *v, const char *rowtype)
5882 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
5883 ss_rtrim (&vs, ss_cstr (" "));
5884 return ss_equals_case (vs, ss_cstr (rowtype));
5888 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
5889 const struct dictionary *d,
5890 const struct variable *rowtype_var,
5891 struct matrix_state *s, size_t si, size_t fi,
5892 size_t cs, size_t cn)
5897 const union value *rowtype_ = case_data (rows[0], rowtype_var);
5898 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
5899 : is_rowtype (rowtype_, "CORR") ? "CR"
5900 : is_rowtype (rowtype_, "MEAN") ? "MN"
5901 : is_rowtype (rowtype_, "STDDEV") ? "SD"
5902 : is_rowtype (rowtype_, "N") ? "NC"
5905 struct string name = DS_EMPTY_INITIALIZER;
5906 ds_put_cstr (&name, name_prefix);
5908 ds_put_format (&name, "F%zu", fi);
5910 ds_put_format (&name, "S%zu", si);
5912 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
5914 mv = matrix_var_create (s, ds_ss (&name));
5917 msg (SW, _("Matrix data file contains variable with existing name %s."),
5922 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
5923 size_t n_missing = 0;
5924 for (size_t y = 0; y < n_rows; y++)
5926 for (size_t x = 0; x < cn; x++)
5928 struct variable *var = dict_get_var (d, cs + x);
5929 double value = case_num (rows[y], var);
5930 if (var_is_num_missing (var, value, MV_ANY))
5935 gsl_matrix_set (m, y, x, value);
5940 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
5941 "value, which was treated as zero.",
5942 "Matrix data file variable %s contains %zu missing "
5943 "values, which were treated as zero.", n_missing),
5944 ds_cstr (&name), n_missing);
5949 for (size_t y = 0; y < n_rows; y++)
5950 case_unref (rows[y]);
5954 var_changed (const struct ccase *ca, const struct ccase *cb,
5955 const struct variable *var)
5957 return !value_equal (case_data (ca, var), case_data (cb, var),
5958 var_get_width (var));
5962 vars_changed (const struct ccase *ca, const struct ccase *cb,
5963 const struct dictionary *d,
5964 size_t first_var, size_t n_vars)
5966 for (size_t i = 0; i < n_vars; i++)
5968 const struct variable *v = dict_get_var (d, first_var + i);
5969 if (var_changed (ca, cb, v))
5976 matrix_cmd_execute_mget (struct mget_command *mget)
5978 struct dictionary *d;
5979 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
5984 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
5985 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
5986 if (!rowtype_ || !varname_)
5989 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
5991 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
5994 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
5996 msg (SE, _("Matrix data file contains no continuous variables."));
6000 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6002 const struct variable *v = dict_get_var (d, i);
6003 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6006 _("Matrix data file contains unexpected string variable %s."),
6012 /* SPLIT variables. */
6014 size_t sn = var_get_dict_index (rowtype_);
6015 struct ccase *sc = NULL;
6018 /* FACTOR variables. */
6019 size_t fs = var_get_dict_index (rowtype_) + 1;
6020 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6021 struct ccase *fc = NULL;
6024 /* Continuous variables. */
6025 size_t cs = var_get_dict_index (varname_) + 1;
6026 size_t cn = dict_get_var_cnt (d) - cs;
6027 struct ccase *cc = NULL;
6030 struct ccase **rows = NULL;
6031 size_t allocated_rows = 0;
6035 while ((c = casereader_read (r)) != NULL)
6037 bool sd = vars_changed (sc, c, d, ss, sn);
6038 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6039 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6054 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6055 mget->state, si, fi, cs, cn);
6061 if (n_rows >= allocated_rows)
6062 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6065 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6066 mget->state, si, fi, cs, cn);
6070 casereader_destroy (r);
6074 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6076 if (!lex_force_id (s->lexer))
6079 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6081 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6086 static struct matrix_cmd *
6087 matrix_parse_eigen (struct matrix_state *s)
6089 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6090 *cmd = (struct matrix_cmd) {
6092 .eigen = { .expr = NULL }
6095 struct eigen_command *eigen = &cmd->eigen;
6096 if (!lex_force_match (s->lexer, T_LPAREN))
6098 eigen->expr = matrix_parse_expr (s);
6100 || !lex_force_match (s->lexer, T_COMMA)
6101 || !matrix_parse_dst_var (s, &eigen->evec)
6102 || !lex_force_match (s->lexer, T_COMMA)
6103 || !matrix_parse_dst_var (s, &eigen->eval)
6104 || !lex_force_match (s->lexer, T_RPAREN))
6110 matrix_cmd_destroy (cmd);
6115 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6117 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6120 if (!is_symmetric (A))
6122 msg (SE, _("Argument of EIGEN must be symmetric."));
6123 gsl_matrix_free (A);
6127 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6128 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6129 gsl_vector v_eval = to_vector (eval);
6130 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6131 gsl_eigen_symmv (A, &v_eval, evec, w);
6132 gsl_eigen_symmv_free (w);
6134 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6136 gsl_matrix_free (A);
6138 gsl_matrix_free (eigen->eval->value);
6139 eigen->eval->value = eval;
6141 gsl_matrix_free (eigen->evec->value);
6142 eigen->evec->value = evec;
6145 static struct matrix_cmd *
6146 matrix_parse_setdiag (struct matrix_state *s)
6148 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6149 *cmd = (struct matrix_cmd) {
6150 .type = MCMD_SETDIAG,
6151 .setdiag = { .dst = NULL }
6154 struct setdiag_command *setdiag = &cmd->setdiag;
6155 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6158 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6161 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6166 if (!lex_force_match (s->lexer, T_COMMA))
6169 setdiag->expr = matrix_parse_expr (s);
6173 if (!lex_force_match (s->lexer, T_RPAREN))
6179 matrix_cmd_destroy (cmd);
6184 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6186 gsl_matrix *dst = setdiag->dst->value;
6189 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6190 setdiag->dst->name);
6194 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6198 size_t n = MIN (dst->size1, dst->size2);
6199 if (is_scalar (src))
6201 double d = to_scalar (src);
6202 for (size_t i = 0; i < n; i++)
6203 gsl_matrix_set (dst, i, i, d);
6205 else if (is_vector (src))
6207 gsl_vector v = to_vector (src);
6208 for (size_t i = 0; i < n && i < v.size; i++)
6209 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6212 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector but it "
6213 "has dimensions (%zu,%zu)."),
6214 src->size1, src->size2);
6215 gsl_matrix_free (src);
6218 static struct matrix_cmd *
6219 matrix_parse_svd (struct matrix_state *s)
6221 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6222 *cmd = (struct matrix_cmd) {
6224 .svd = { .expr = NULL }
6227 struct svd_command *svd = &cmd->svd;
6228 if (!lex_force_match (s->lexer, T_LPAREN))
6230 svd->expr = matrix_parse_expr (s);
6232 || !lex_force_match (s->lexer, T_COMMA)
6233 || !matrix_parse_dst_var (s, &svd->u)
6234 || !lex_force_match (s->lexer, T_COMMA)
6235 || !matrix_parse_dst_var (s, &svd->s)
6236 || !lex_force_match (s->lexer, T_COMMA)
6237 || !matrix_parse_dst_var (s, &svd->v)
6238 || !lex_force_match (s->lexer, T_RPAREN))
6244 matrix_cmd_destroy (cmd);
6249 matrix_cmd_execute_svd (struct svd_command *svd)
6251 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6255 if (m->size1 >= m->size2)
6258 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6259 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6260 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6261 gsl_vector *work = gsl_vector_alloc (A->size2);
6262 gsl_linalg_SV_decomp (A, V, &Sv, work);
6263 gsl_vector_free (work);
6265 matrix_var_set (svd->u, A);
6266 matrix_var_set (svd->s, S);
6267 matrix_var_set (svd->v, V);
6271 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6272 gsl_matrix_transpose_memcpy (At, m);
6273 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6274 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6275 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6276 gsl_vector *work = gsl_vector_alloc (At->size2);
6277 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6278 gsl_vector_free (work);
6280 matrix_var_set (svd->v, At);
6281 matrix_var_set (svd->s, St);
6282 matrix_var_set (svd->u, Vt);
6287 matrix_cmd_execute (struct matrix_cmd *cmd)
6292 matrix_cmd_execute_compute (&cmd->compute);
6296 matrix_cmd_execute_print (&cmd->print);
6300 return matrix_cmd_execute_do_if (&cmd->do_if);
6303 matrix_cmd_execute_loop (&cmd->loop);
6310 matrix_cmd_execute_display (&cmd->display);
6314 matrix_cmd_execute_release (&cmd->release);
6318 matrix_cmd_execute_save (&cmd->save);
6322 matrix_cmd_execute_read (&cmd->read);
6326 matrix_cmd_execute_write (&cmd->write);
6330 matrix_cmd_execute_get (&cmd->get);
6334 matrix_cmd_execute_msave (&cmd->msave);
6338 matrix_cmd_execute_mget (&cmd->mget);
6342 matrix_cmd_execute_eigen (&cmd->eigen);
6346 matrix_cmd_execute_setdiag (&cmd->setdiag);
6350 matrix_cmd_execute_svd (&cmd->svd);
6357 struct matrix_command_name
6360 struct matrix_cmd *(*parse) (struct matrix_state *);
6363 static const struct matrix_command_name *
6364 matrix_parse_command_name (struct lexer *lexer)
6366 static const struct matrix_command_name commands[] = {
6367 { "COMPUTE", matrix_parse_compute },
6368 { "CALL EIGEN", matrix_parse_eigen },
6369 { "CALL SETDIAG", matrix_parse_setdiag },
6370 { "CALL SVD", matrix_parse_svd },
6371 { "PRINT", matrix_parse_print },
6372 { "DO IF", matrix_parse_do_if },
6373 { "LOOP", matrix_parse_loop },
6374 { "BREAK", matrix_parse_break },
6375 { "READ", matrix_parse_read },
6376 { "WRITE", matrix_parse_write },
6377 { "GET", matrix_parse_get },
6378 { "SAVE", matrix_parse_save },
6379 { "MGET", matrix_parse_mget },
6380 { "MSAVE", matrix_parse_msave },
6381 { "DISPLAY", matrix_parse_display },
6382 { "RELEASE", matrix_parse_release },
6384 static size_t n = sizeof commands / sizeof *commands;
6386 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6387 if (lex_match_phrase (lexer, c->name))
6392 static struct matrix_cmd *
6393 matrix_parse_command (struct matrix_state *s)
6395 size_t nesting_level = SIZE_MAX;
6397 struct matrix_cmd *c = NULL;
6398 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6400 lex_error (s->lexer, _("Unknown matrix command."));
6401 else if (!cmd->parse)
6402 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6406 nesting_level = output_open_group (
6407 group_item_create_nocopy (utf8_to_title (cmd->name),
6408 utf8_to_title (cmd->name)));
6413 lex_end_of_command (s->lexer);
6414 lex_discard_rest_of_command (s->lexer);
6415 if (nesting_level != SIZE_MAX)
6416 output_close_groups (nesting_level);
6422 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6424 if (!lex_force_match (lexer, T_ENDCMD))
6427 struct matrix_state state = {
6428 .session = dataset_session (ds),
6430 .vars = HMAP_INITIALIZER (state.vars),
6435 while (lex_match (lexer, T_ENDCMD))
6437 if (lex_token (lexer) == T_STOP)
6439 msg (SE, _("Unexpected end of input expecting matrix command."));
6443 if (lex_match_phrase (lexer, "END MATRIX"))
6446 struct matrix_cmd *c = matrix_parse_command (&state);
6449 matrix_cmd_execute (c);
6450 matrix_cmd_destroy (c);
6456 dict_unref (state.common->dict);
6457 casewriter_destroy (state.common->writer);
6458 free (state.common);