From: Ben Pfaff Date: Mon, 27 Sep 2021 05:35:44 +0000 (-0700) Subject: MATRIX X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=39add622928e4466ffbbfd6792a49e1a4013fed8;p=pspp MATRIX --- diff --git a/src/data/any-writer.c b/src/data/any-writer.c index b26440158c..ba4149d78b 100644 --- a/src/data/any-writer.c +++ b/src/data/any-writer.c @@ -37,7 +37,8 @@ #include "gettext.h" #define _(msgid) gettext (msgid) -/* Creates and returns a writer for HANDLE with the given DICT. */ +/* Creates and returns a writer for HANDLE with the given DICT. The caller + retains ownership of DICT, which might get modified to assign short names. */ struct casewriter * any_writer_open (struct file_handle *handle, struct dictionary *dict) { diff --git a/src/data/identifier.c b/src/data/identifier.c index 07c7675bd5..baa0abd1c5 100644 --- a/src/data/identifier.c +++ b/src/data/identifier.c @@ -99,9 +99,21 @@ token_type_to_string (enum token_type token) case T_RBRACK: return "]"; + case T_LCURLY: + return "{"; + + case T_RCURLY: + return "}"; + case T_COMMA: return ","; + case T_SEMICOLON: + return ";"; + + case T_COLON: + return ":"; + case T_AND: return "AND"; diff --git a/src/data/identifier.h b/src/data/identifier.h index d694a5201d..dcbce970cd 100644 --- a/src/data/identifier.h +++ b/src/data/identifier.h @@ -41,7 +41,11 @@ TOKEN_TYPE(RPAREN) /* ) */ \ TOKEN_TYPE(LBRACK) /* [ */ \ TOKEN_TYPE(RBRACK) /* ] */ \ + TOKEN_TYPE(LCURLY) /* { */ \ + TOKEN_TYPE(RCURLY) /* } */ \ TOKEN_TYPE(COMMA) /* , */ \ + TOKEN_TYPE(SEMICOLON) /* ; */ \ + TOKEN_TYPE(COLON) /* : */ \ \ TOKEN_TYPE(AND) /* AND */ \ TOKEN_TYPE(OR) /* OR */ \ diff --git a/src/data/settings.c b/src/data/settings.c index 313fed0baf..4bc36494fd 100644 --- a/src/data/settings.c +++ b/src/data/settings.c @@ -48,6 +48,9 @@ struct settings /* Format of reals in output (SET WRB). */ enum float_format output_float_format; + /* MATRIX...END MATRIX settings. */ + enum settings_mdisplay mdisplay; + int viewlength; int viewwidth; bool safer_mode; @@ -90,6 +93,7 @@ static struct settings the_settings = { .input_float_format = FLOAT_NATIVE_DOUBLE, .output_integer_format = INTEGER_NATIVE, .output_float_format = FLOAT_NATIVE_DOUBLE, + .mdisplay = SETTINGS_MDISPLAY_TEXT, .viewlength = 24, .viewwidth = 79, .safer_mode = false, @@ -701,3 +705,15 @@ settings_set_show_variables (enum settings_value_show s) { the_settings.show_variables = s; } + +enum settings_mdisplay +settings_get_mdisplay (void) +{ + return the_settings.mdisplay; +} + +void +settings_set_mdisplay (enum settings_mdisplay mdisplay) +{ + the_settings.mdisplay = mdisplay; +} diff --git a/src/data/settings.h b/src/data/settings.h index d7c616d9e0..9f6b94a888 100644 --- a/src/data/settings.h +++ b/src/data/settings.h @@ -187,4 +187,13 @@ void settings_set_output_routing (enum settings_output_type, enum settings_output_devices settings_get_output_routing ( enum settings_output_type); +enum settings_mdisplay + { + SETTINGS_MDISPLAY_TEXT, + SETTINGS_MDISPLAY_TABLES + }; + +enum settings_mdisplay settings_get_mdisplay (void); +void settings_set_mdisplay (enum settings_mdisplay); + #endif /* !settings_h */ diff --git a/src/language/command.def b/src/language/command.def index 9b3ff8a729..88cbaf5cbe 100644 --- a/src/language/command.def +++ b/src/language/command.def @@ -30,6 +30,7 @@ DEF_CMD (S_ANY, 0, "FINISH", cmd_finish) DEF_CMD (S_ANY, 0, "HOST", cmd_host) DEF_CMD (S_ANY, 0, "INCLUDE", cmd_include) DEF_CMD (S_ANY, 0, "INSERT", cmd_insert) +DEF_CMD (S_ANY, 0, "MATRIX", cmd_matrix) DEF_CMD (S_ANY, 0, "MCONVERT", cmd_mconvert) DEF_CMD (S_ANY, 0, "N OF CASES", cmd_n_of_cases) DEF_CMD (S_ANY, F_ABBREV, "N", cmd_n_of_cases) @@ -210,7 +211,6 @@ UNIMPL_CMD ("KM", "Kaplan-Meier") UNIMPL_CMD ("LOGLINEAR", "General model fitting") UNIMPL_CMD ("MANOVA", "Multivariate analysis of variance") UNIMPL_CMD ("MAPS", "Geographical display") -UNIMPL_CMD ("MATRIX", "Matrix processing") UNIMPL_CMD ("MIXED", "Mixed linear models") UNIMPL_CMD ("MODEL CLOSE", "Close server connection") UNIMPL_CMD ("MODEL HANDLE", "Define server connection") diff --git a/src/language/data-io/data-writer.c b/src/language/data-io/data-writer.c index bf9505e7dd..f38bc6896d 100644 --- a/src/language/data-io/data-writer.c +++ b/src/language/data-io/data-writer.c @@ -203,6 +203,25 @@ dfm_put_record (struct dfm_writer *w, const char *rec, size_t len) return !dfm_write_error (w); } +/* Writes record REC (which need not be null-terminated) having length LEN to + the file corresponding to HANDLE. REC is encoded in UTF-8, which this + function recodes to the correct encoding for W before writing. Adds any + needed formatting, such as a trailing new-line. Returns true on success, + false on failure. */ +bool +dfm_put_record_utf8 (struct dfm_writer *w, const char *rec, size_t len) +{ + if (is_encoding_utf8 (w->encoding)) + return dfm_put_record (w, rec, len); + else + { + char *recoded = recode_string (w->encoding, UTF8, rec, len); + bool ok = dfm_put_record (w, recoded, strlen (recoded)); + free (recoded); + return ok; + } +} + /* Closes data file writer W. */ bool dfm_close_writer (struct dfm_writer *w) diff --git a/src/language/data-io/data-writer.h b/src/language/data-io/data-writer.h index 10ad6cd656..573f4a0dd1 100644 --- a/src/language/data-io/data-writer.h +++ b/src/language/data-io/data-writer.h @@ -28,6 +28,7 @@ struct dfm_writer *dfm_open_writer (struct file_handle *, bool dfm_close_writer (struct dfm_writer *); bool dfm_write_error (const struct dfm_writer *); bool dfm_put_record (struct dfm_writer *, const char *rec, size_t len); +bool dfm_put_record_utf8 (struct dfm_writer *, const char *rec, size_t len); const char *dfm_writer_get_encoding (const struct dfm_writer *); #endif /* data-writer.h */ diff --git a/src/language/data-io/print.c b/src/language/data-io/print.c index 9e4d8a294d..c0fc3d84f1 100644 --- a/src/language/data-io/print.c +++ b/src/language/data-io/print.c @@ -575,14 +575,7 @@ print_text_flush_records (struct print_trns *trns, struct u8_line *line, len--; } - if (is_encoding_utf8 (trns->encoding)) - dfm_put_record (trns->writer, s, len); - else - { - char *recoded = recode_string (trns->encoding, UTF8, s, len); - dfm_put_record (trns->writer, recoded, strlen (recoded)); - free (recoded); - } + dfm_put_record_utf8 (trns->writer, s, len); } } } diff --git a/src/language/expressions/evaluate.c b/src/language/expressions/evaluate.c index 7543aed55e..d6d5738481 100644 --- a/src/language/expressions/evaluate.c +++ b/src/language/expressions/evaluate.c @@ -298,6 +298,5 @@ expr_debug_print_postfix (const struct expression *e) NOT_REACHED (); } } - output_log ("%s", ds_cstr (&s)); - ds_destroy (&s); + output_log_nocopy (ds_steal_cstr (&s)); } diff --git a/src/language/lexer/lexer.c b/src/language/lexer/lexer.c index 2e6232afdc..c9966e8c02 100644 --- a/src/language/lexer/lexer.c +++ b/src/language/lexer/lexer.c @@ -379,6 +379,14 @@ lex_get (struct lexer *lexer) return; } } + +/* Advances LEXER by N tokens. */ +void +lex_get_n (struct lexer *lexer, size_t n) +{ + while (n-- > 0) + lex_get (lexer); +} /* Issuing errors. */ @@ -564,7 +572,8 @@ lex_next_error_valist (struct lexer *lexer, int n0, int n1, ds_put_cstr (&s, ": "); ds_put_vformat (&s, format, args); } - ds_put_byte (&s, '.'); + if (ds_last (&s) != '.') + ds_put_byte (&s, '.'); msg (SE, "%s", ds_cstr (&s)); ds_destroy (&s); } @@ -1123,32 +1132,49 @@ lex_tokens_match (const struct token *actual, const struct token *expected) } } -/* If LEXER is positioned at the sequence of tokens that may be parsed from S, - skips it and returns true. Otherwise, returns false. - - S may consist of an arbitrary sequence of tokens, e.g. "KRUSKAL-WALLIS", - "2SLS", or "END INPUT PROGRAM". Identifiers may be abbreviated to their - first three letters. */ -bool -lex_match_phrase (struct lexer *lexer, const char *s) +static size_t +lex_at_phrase__ (struct lexer *lexer, const char *s) { struct string_lexer slex; struct token token; - int i; - i = 0; + size_t i = 0; string_lexer_init (&slex, s, strlen (s), SEG_MODE_INTERACTIVE, true); while (string_lexer_next (&slex, &token)) { bool match = lex_tokens_match (lex_next (lexer, i++), &token); token_uninit (&token); if (!match) - return false; + return 0; } + return i; +} - while (i-- > 0) - lex_get (lexer); - return true; +/* If LEXER is positioned at the sequence of tokens that may be parsed from S, + returns true. Otherwise, returns false. + + S may consist of an arbitrary sequence of tokens, e.g. "KRUSKAL-WALLIS", + "2SLS", or "END INPUT PROGRAM". Identifiers may be abbreviated to their + first three letters. */ +bool +lex_at_phrase (struct lexer *lexer, const char *s) +{ + return lex_at_phrase__ (lexer, s) > 0; +} + +/* If LEXER is positioned at the sequence of tokens that may be parsed from S, + skips it and returns true. Otherwise, returns false. + + S may consist of an arbitrary sequence of tokens, e.g. "KRUSKAL-WALLIS", + "2SLS", or "END INPUT PROGRAM". Identifiers may be abbreviated to their + first three letters. */ +bool +lex_match_phrase (struct lexer *lexer, const char *s) +{ + size_t n = lex_at_phrase__ (lexer, s); + if (n > 0) + lex_get_n (lexer, n); + return n > 0; } static int diff --git a/src/language/lexer/lexer.h b/src/language/lexer/lexer.h index 6aa900e8df..1282b6946b 100644 --- a/src/language/lexer/lexer.h +++ b/src/language/lexer/lexer.h @@ -100,6 +100,7 @@ void lex_append (struct lexer *, struct lex_reader *); /* Advancing. */ void lex_get (struct lexer *); +void lex_get_n (struct lexer *, size_t n); /* Token testing functions. */ bool lex_is_number (const struct lexer *); @@ -120,6 +121,7 @@ bool lex_match (struct lexer *, enum token_type); bool lex_match_id (struct lexer *, const char *); bool lex_match_id_n (struct lexer *, const char *, size_t n); bool lex_match_int (struct lexer *, int); +bool lex_at_phrase (struct lexer *, const char *s); bool lex_match_phrase (struct lexer *, const char *s); /* Forcible matching functions. */ diff --git a/src/language/lexer/macro.c b/src/language/lexer/macro.c index 814740bfac..4419306d1d 100644 --- a/src/language/lexer/macro.c +++ b/src/language/lexer/macro.c @@ -334,6 +334,8 @@ classify_token (enum token_type type) case T_RPAREN: case T_LBRACK: case T_RBRACK: + case T_LCURLY: + case T_RCURLY: return TC_PUNCT; case T_PLUS: @@ -341,6 +343,7 @@ classify_token (enum token_type type) case T_ASTERISK: case T_SLASH: case T_EQUALS: + case T_COLON: case T_AND: case T_OR: case T_NOT: @@ -359,6 +362,7 @@ classify_token (enum token_type type) return TC_BINOP; case T_COMMA: + case T_SEMICOLON: return TC_COMMA; } diff --git a/src/language/lexer/scan.c b/src/language/lexer/scan.c index b6baf5ce55..4e15cca6c6 100644 --- a/src/language/lexer/scan.c +++ b/src/language/lexer/scan.c @@ -188,6 +188,8 @@ scan_punct1__ (char c0) case '-': return T_DASH; case '[': return T_LBRACK; case ']': return T_RBRACK; + case '{': return T_LCURLY; + case '}': return T_RCURLY; case '&': return T_AND; case '|': return T_OR; case '+': return T_PLUS; @@ -196,6 +198,8 @@ scan_punct1__ (char c0) case '<': return T_LT; case '>': return T_GT; case '~': return T_NOT; + case ';': return T_SEMICOLON; + case ':': return T_COLON; default: return T_MACRO_PUNCT; } diff --git a/src/language/lexer/segment.c b/src/language/lexer/segment.c index 519f6ec9f2..39604d513e 100644 --- a/src/language/lexer/segment.c +++ b/src/language/lexer/segment.c @@ -956,8 +956,8 @@ segmenter_parse_mid_command__ (struct segmenter *s, return segmenter_parse_number__ (s, input, n, eof, type, ofs); } /* Fall through. */ - case '(': case ')': case ',': case '=': - case '[': case ']': case '&': case '|': case '+': + case '(': case ')': case '{': case ',': case '=': case ';': case ':': + case '[': case ']': case '}': case '&': case '|': case '+': *type = SEG_PUNCT; s->substate = 0; return 1; diff --git a/src/language/stats/automake.mk b/src/language/stats/automake.mk index 88c64abece..460e95e707 100644 --- a/src/language/stats/automake.mk +++ b/src/language/stats/automake.mk @@ -51,6 +51,7 @@ language_stats_sources = \ src/language/stats/jonckheere-terpstra.h \ src/language/stats/mann-whitney.c \ src/language/stats/mann-whitney.h \ + src/language/stats/matrix.c \ src/language/stats/means.c \ src/language/stats/means.h \ src/language/stats/means-calc.c \ diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c new file mode 100644 index 0000000000..c0c80eda77 --- /dev/null +++ b/src/language/stats/matrix.c @@ -0,0 +1,6462 @@ +/* PSPP - a program for statistical analysis. + Copyright (C) 2021 Free Software Foundation, Inc. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "data/any-reader.h" +#include "data/any-writer.h" +#include "data/casereader.h" +#include "data/casewriter.h" +#include "data/data-in.h" +#include "data/data-out.h" +#include "data/dataset.h" +#include "data/dictionary.h" +#include "data/file-handle-def.h" +#include "language/command.h" +#include "language/data-io/data-reader.h" +#include "language/data-io/data-writer.h" +#include "language/data-io/file-handle.h" +#include "language/lexer/format-parser.h" +#include "language/lexer/lexer.h" +#include "language/lexer/variable-parser.h" +#include "libpspp/array.h" +#include "libpspp/assertion.h" +#include "libpspp/compiler.h" +#include "libpspp/hmap.h" +#include "libpspp/i18n.h" +#include "libpspp/misc.h" +#include "libpspp/str.h" +#include "libpspp/string-array.h" +#include "libpspp/stringi-set.h" +#include "libpspp/u8-line.h" +#include "math/random.h" +#include "output/driver.h" +#include "output/output-item.h" +#include "output/pivot-table.h" + +#include "gl/c-ctype.h" +#include "gl/minmax.h" +#include "gl/xsize.h" + +#include "gettext.h" +#define _(msgid) gettext (msgid) +#define N_(msgid) (msgid) + +struct matrix_var + { + struct hmap_node hmap_node; + const char *name; + gsl_matrix *value; + }; + +struct msave_common + { + /* Configuration. */ + struct file_handle *outfile; + struct string_array variables; + struct string_array fnames; + struct string_array snames; + bool has_factors; + bool has_splits; + size_t n_varnames; + + /* Execution state. */ + struct dictionary *dict; + struct casewriter *writer; + }; + +struct read_file + { + struct file_handle *file; + struct dfm_reader *reader; + char *encoding; + }; + +struct matrix_state + { + struct dataset *dataset; + struct session *session; + struct lexer *lexer; + struct hmap vars; + bool in_loop; + struct file_handle *prev_save_outfile; + struct file_handle *prev_write_outfile; + struct msave_common *common; + + struct file_handle *prev_read_file; + struct read_file **read_files; + size_t n_read_files; + }; + +static struct matrix_var * +matrix_var_lookup (struct matrix_state *s, struct substring name) +{ + struct matrix_var *var; + + HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node, + utf8_hash_case_substring (name, 0), &s->vars) + if (!utf8_sscasecmp (ss_cstr (var->name), name)) + return var; + + return NULL; +} + +static struct matrix_var * +matrix_var_create (struct matrix_state *s, struct substring name) +{ + struct matrix_var *var = xmalloc (sizeof *var); + *var = (struct matrix_var) { .name = ss_xstrdup (name) }; + hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0)); + return var; +} + +static void +matrix_var_set (struct matrix_var *var, gsl_matrix *value) +{ + gsl_matrix_free (var->value); + var->value = value; +} + +#define MATRIX_FUNCTIONS \ + F(ABS, m_m) \ + F(ALL, d_m) \ + F(ANY, d_m) \ + F(ARSIN, m_m) \ + F(ARTAN, m_m) \ + F(BLOCK, m_any) \ + F(CHOL, m_m) \ + F(CMIN, m_m) \ + F(CMAX, m_m) \ + F(COS, m_m) \ + F(CSSQ, m_m) \ + F(CSUM, m_m) \ + F(DESIGN, m_m) \ + F(DET, d_m) \ + F(DIAG, m_m) \ + F(EVAL, m_m) \ + F(EXP, m_m) \ + F(GINV, m_m) \ + F(GRADE, m_m) \ + F(GSCH, m_m) \ + F(IDENT, IDENT) \ + F(INV, m_m) \ + F(KRONEKER, m_mm) \ + F(LG10, m_m) \ + F(LN, m_m) \ + F(MAGIC, m_d) \ + F(MAKE, m_ddd) \ + F(MDIAG, m_v) \ + F(MMAX, d_m) \ + F(MMIN, d_m) \ + F(MOD, m_md) \ + F(MSSQ, d_m) \ + F(MSUM, d_m) \ + F(NCOL, d_m) \ + F(NROW, d_m) \ + F(RANK, d_m) \ + F(RESHAPE, m_mdd) \ + F(RMAX, m_m) \ + F(RMIN, m_m) \ + F(RND, m_m) \ + F(RNKORDER, m_m) \ + F(RSSQ, m_m) \ + F(RSUM, m_m) \ + F(SIN, m_m) \ + F(SOLVE, m_mm) \ + F(SQRT, m_m) \ + F(SSCP, m_m) \ + F(SVAL, m_m) \ + F(SWEEP, m_md) \ + F(T, m_m) \ + F(TRACE, d_m) \ + F(TRANSPOS, m_m) \ + F(TRUNC, m_m) \ + F(UNIFORM, m_dd) + +enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 }; +enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 }; +enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 }; +enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 }; +enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 }; +enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 }; +enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 }; +enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 }; +enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 }; +enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX }; +enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 }; + +typedef gsl_matrix *matrix_proto_m_d (double); +typedef gsl_matrix *matrix_proto_m_dd (double, double); +typedef gsl_matrix *matrix_proto_m_ddd (double, double, double); +typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *); +typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double); +typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double); +typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *); +typedef gsl_matrix *matrix_proto_m_v (gsl_vector *); +typedef double matrix_proto_d_m (gsl_matrix *); +typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n); +typedef gsl_matrix *matrix_proto_IDENT (double, double); + +#define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME; +MATRIX_FUNCTIONS +#undef F + +/* Matrix expressions. */ + +struct matrix_expr + { + enum matrix_op + { + /* Functions. */ +#define F(NAME, PROTOTYPE) MOP_F_##NAME, + MATRIX_FUNCTIONS +#undef F + + /* Elementwise and scalar arithmetic. */ + MOP_NEGATE, /* unary - */ + MOP_ADD_ELEMS, /* + */ + MOP_SUB_ELEMS, /* - */ + MOP_MUL_ELEMS, /* &* */ + MOP_DIV_ELEMS, /* / and &/ */ + MOP_EXP_ELEMS, /* &** */ + MOP_SEQ, /* a:b */ + MOP_SEQ_BY, /* a:b:c */ + + /* Matrix arithmetic. */ + MOP_MUL_MAT, /* * */ + MOP_EXP_MAT, /* ** */ + + /* Relational. */ + MOP_GT, /* > */ + MOP_GE, /* >= */ + MOP_LT, /* < */ + MOP_LE, /* <= */ + MOP_EQ, /* = */ + MOP_NE, /* <> */ + + /* Logical. */ + MOP_NOT, /* NOT */ + MOP_AND, /* AND */ + MOP_OR, /* OR */ + MOP_XOR, /* XOR */ + + /* {}. */ + MOP_PASTE_HORZ, /* a, b, c, ... */ + MOP_PASTE_VERT, /* a; b; c; ... */ + MOP_EMPTY, /* {} */ + + /* Sub-matrices. */ + MOP_VEC_INDEX, /* x(y) */ + MOP_VEC_ALL, /* x(:) */ + MOP_MAT_INDEX, /* x(y,z) */ + MOP_ROW_INDEX, /* x(y,:) */ + MOP_COL_INDEX, /* x(:,z) */ + + /* Literals. */ + MOP_NUMBER, + MOP_VARIABLE, + + /* Oddball stuff. */ + MOP_EOF, /* EOF('file') */ + } + op; + + union + { + struct + { + struct matrix_expr **subs; + size_t n_subs; + }; + + double number; + struct matrix_var *variable; + struct read_file *eof; + }; + }; + +static void +matrix_expr_destroy (struct matrix_expr *e) +{ + if (!e) + return; + + switch (e->op) + { +#define F(NAME, PROTOTYPE) case MOP_F_##NAME: +MATRIX_FUNCTIONS +#undef F + case MOP_NEGATE: + case MOP_ADD_ELEMS: + case MOP_SUB_ELEMS: + case MOP_MUL_ELEMS: + case MOP_DIV_ELEMS: + case MOP_EXP_ELEMS: + case MOP_SEQ: + case MOP_SEQ_BY: + case MOP_MUL_MAT: + case MOP_EXP_MAT: + case MOP_GT: + case MOP_GE: + case MOP_LT: + case MOP_LE: + case MOP_EQ: + case MOP_NE: + case MOP_NOT: + case MOP_AND: + case MOP_OR: + case MOP_XOR: + case MOP_EMPTY: + case MOP_PASTE_HORZ: + case MOP_PASTE_VERT: + case MOP_VEC_INDEX: + case MOP_VEC_ALL: + case MOP_MAT_INDEX: + case MOP_ROW_INDEX: + case MOP_COL_INDEX: + for (size_t i = 0; i < e->n_subs; i++) + matrix_expr_destroy (e->subs[i]); + break; + + case MOP_NUMBER: + case MOP_VARIABLE: + case MOP_EOF: + break; + } + free (e); +} + +static struct matrix_expr * +matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs, + size_t n_subs) +{ + struct matrix_expr *e = xmalloc (sizeof *e); + *e = (struct matrix_expr) { + .op = op, + .subs = xmemdup (subs, n_subs * sizeof *subs), + .n_subs = n_subs + }; + return e; +} + +static struct matrix_expr * +matrix_expr_create_0 (enum matrix_op op) +{ + struct matrix_expr *sub; + return matrix_expr_create_subs (op, &sub, 0); +} + +static struct matrix_expr * +matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub) +{ + return matrix_expr_create_subs (op, &sub, 1); +} + +static struct matrix_expr * +matrix_expr_create_2 (enum matrix_op op, + struct matrix_expr *sub0, struct matrix_expr *sub1) +{ + struct matrix_expr *subs[] = { sub0, sub1 }; + return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs); +} + +static struct matrix_expr * +matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0, + struct matrix_expr *sub1, struct matrix_expr *sub2) +{ + struct matrix_expr *subs[] = { sub0, sub1, sub2 }; + return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs); +} + +static struct matrix_expr * +matrix_expr_create_number (double number) +{ + struct matrix_expr *e = xmalloc (sizeof *e); + *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number }; + return e; +} + +static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *); + +static struct matrix_expr * +matrix_parse_curly_comma (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_or_xor (s); + if (!lhs) + return NULL; + + while (lex_match (s->lexer, T_COMMA)) + { + struct matrix_expr *rhs = matrix_parse_or_xor (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs); + } + return lhs; +} + +static struct matrix_expr * +matrix_parse_curly_semi (struct matrix_state *s) +{ + if (lex_token (s->lexer) == T_RCURLY) + return matrix_expr_create_0 (MOP_EMPTY); + + struct matrix_expr *lhs = matrix_parse_curly_comma (s); + if (!lhs) + return NULL; + + while (lex_match (s->lexer, T_SEMICOLON)) + { + struct matrix_expr *rhs = matrix_parse_curly_comma (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs); + } + return lhs; +} + +#define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \ + for (size_t Y = 0; Y < (M)->size1; Y++) \ + for (size_t X = 0; X < (M)->size2; X++) \ + for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL) + +static bool +is_vector (gsl_matrix *m) +{ + return m->size1 <= 1 || m->size2 <= 1; +} + +static gsl_vector +to_vector (gsl_matrix *m) +{ + return (m->size1 == 1 + ? gsl_matrix_row (m, 0).vector + : gsl_matrix_column (m, 0).vector); +} + + +static gsl_matrix * +matrix_eval_ABS (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = fabs (*d); + return m; +} + +static double +matrix_eval_ALL (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + if (*d == 0.0) + return 0.0; + return 1.0; +} + +static double +matrix_eval_ANY (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + if (*d != 0.0) + return !.0; + return 0.0; +} + +static gsl_matrix * +matrix_eval_ARSIN (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = asin (*d); + return m; +} + +static gsl_matrix * +matrix_eval_ARTAN (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = atan (*d); + return m; +} + +static gsl_matrix * +matrix_eval_BLOCK (gsl_matrix *m[], size_t n) +{ + size_t r = 0; + size_t c = 0; + for (size_t i = 0; i < n; i++) + { + r += m[i]->size1; + c += m[i]->size2; + } + gsl_matrix *block = gsl_matrix_calloc (r, c); + r = c = 0; + for (size_t i = 0; i < n; i++) + { + for (size_t y = 0; y < m[i]->size1; y++) + for (size_t x = 0; x < m[i]->size2; x++) + gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x)); + r += m[i]->size1; + c += m[i]->size2; + } + return block; +} + +static gsl_matrix * +matrix_eval_CHOL (gsl_matrix *m) +{ + if (!gsl_linalg_cholesky_decomp1 (m)) + { + for (size_t y = 0; y < m->size1; y++) + for (size_t x = y + 1; x < m->size2; x++) + gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y)); + + for (size_t y = 0; y < m->size1; y++) + for (size_t x = 0; x < y; x++) + gsl_matrix_set (m, y, x, 0); + return m; + } + else + { + msg (SE, _("Input to CHOL function is not positive-definite.")); + return NULL; + } +} + +static gsl_matrix * +matrix_eval_col_extremum (gsl_matrix *m, bool min) +{ + if (m->size1 <= 1) + return m; + else if (!m->size2) + return gsl_matrix_alloc (1, 0); + + gsl_matrix *cext = gsl_matrix_alloc (1, m->size2); + for (size_t x = 0; x < m->size2; x++) + { + double ext = gsl_matrix_get (m, 0, x); + for (size_t y = 1; y < m->size1; y++) + { + double value = gsl_matrix_get (m, y, x); + if (min ? value < ext : value > ext) + ext = value; + } + gsl_matrix_set (cext, 0, x, ext); + } + return cext; +} + +static gsl_matrix * +matrix_eval_CMAX (gsl_matrix *m) +{ + return matrix_eval_col_extremum (m, false); +} + +static gsl_matrix * +matrix_eval_CMIN (gsl_matrix *m) +{ + return matrix_eval_col_extremum (m, true); +} + +static gsl_matrix * +matrix_eval_COS (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = cos (*d); + return m; +} + +static gsl_matrix * +matrix_eval_col_sum (gsl_matrix *m, bool square) +{ + if (m->size1 == 0) + return m; + else if (!m->size2) + return gsl_matrix_alloc (1, 0); + + gsl_matrix *result = gsl_matrix_alloc (1, m->size2); + for (size_t x = 0; x < m->size2; x++) + { + double sum = 0; + for (size_t y = 0; y < m->size1; y++) + { + double d = gsl_matrix_get (m, y, x); + sum += square ? pow2 (d) : d; + } + gsl_matrix_set (result, 0, x, sum); + } + return result; +} + +static gsl_matrix * +matrix_eval_CSSQ (gsl_matrix *m) +{ + return matrix_eval_col_sum (m, true); +} + +static gsl_matrix * +matrix_eval_CSUM (gsl_matrix *m) +{ + return matrix_eval_col_sum (m, false); +} + +static int +compare_double_3way (const void *a_, const void *b_) +{ + const double *a = a_; + const double *b = b_; + return *a < *b ? -1 : *a > *b; +} + +static gsl_matrix * +matrix_eval_DESIGN (gsl_matrix *m) +{ + double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp); + gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix; + gsl_matrix_transpose_memcpy (&m2, m); + + for (size_t y = 0; y < m2.size1; y++) + qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way); + + size_t *n = xcalloc (m2.size1, sizeof *n); + size_t n_total = 0; + for (size_t i = 0; i < m2.size1; i++) + { + double *row = tmp + m2.size2 * i; + for (size_t j = 0; j < m2.size2; ) + { + size_t k; + for (k = j + 1; k < m2.size2; k++) + if (row[j] != row[k]) + break; + row[n[i]++] = row[j]; + j = k; + } + + if (n[i] <= 1) + msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1); + else + n_total += n[i]; + } + + gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total); + size_t x = 0; + for (size_t i = 0; i < m->size2; i++) + { + if (n[i] <= 1) + continue; + + const double *unique = tmp + m2.size2 * i; + for (size_t j = 0; j < n[i]; j++, x++) + { + double value = unique[j]; + for (size_t y = 0; y < m->size1; y++) + gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value); + } + } + + free (n); + free (tmp); + + return result; +} + +static double +matrix_eval_DET (gsl_matrix *m) +{ + gsl_permutation *p = gsl_permutation_alloc (m->size1); + int signum; + gsl_linalg_LU_decomp (m, p, &signum); + gsl_permutation_free (p); + return gsl_linalg_LU_det (m, signum); +} + +static gsl_matrix * +matrix_eval_DIAG (gsl_matrix *m) +{ + gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1); + for (size_t i = 0; i < diag->size1; i++) + gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i)); + return diag; +} + +static bool +is_symmetric (const gsl_matrix *m) +{ + if (m->size1 != m->size2) + return false; + + for (size_t y = 0; y < m->size1; y++) + for (size_t x = 0; x < y; x++) + if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y)) + return false; + + return true; +} + +static int +compare_double_desc (const void *a_, const void *b_) +{ + const double *a = a_; + const double *b = b_; + return *a > *b ? -1 : *a < *b; +} + +static gsl_matrix * +matrix_eval_EVAL (gsl_matrix *m) +{ + if (!is_symmetric (m)) + { + msg (SE, _("Argument of EVAL must be symmetric.")); + return NULL; + } + + gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1); + gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1); + gsl_vector v_eval = to_vector (eval); + gsl_eigen_symm (m, &v_eval, w); + gsl_eigen_symm_free (w); + + assert (v_eval.stride == 1); + qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc); + + return eval; +} + +static gsl_matrix * +matrix_eval_EXP (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = exp (*d); + return m; +} + +/* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was + marked as: + + Charl Linssen + Feb 2016 + PUBLIC DOMAIN */ +static gsl_matrix * +matrix_eval_GINV (gsl_matrix *A) +{ + size_t n = A->size1; + size_t m = A->size2; + bool swap = m > n; + gsl_matrix *tmp_mat = NULL; + if (swap) + { + /* libgsl SVD can only handle the case m <= n, so transpose matrix. */ + tmp_mat = gsl_matrix_alloc (m, n); + gsl_matrix_transpose_memcpy (tmp_mat, A); + A = tmp_mat; + size_t i = m; + m = n; + n = i; + } + + /* Do SVD. */ + gsl_matrix *V = gsl_matrix_alloc (m, m); + gsl_vector *u = gsl_vector_alloc (m); + + gsl_vector *tmp_vec = gsl_vector_alloc (m); + gsl_linalg_SV_decomp (A, V, u, tmp_vec); + gsl_vector_free (tmp_vec); + + /* Compute Σ⁻¹. */ + gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n); + gsl_matrix_set_zero (Sigma_pinv); + double cutoff = 1e-15 * gsl_vector_max (u); + + for (size_t i = 0; i < m; ++i) + { + double x = gsl_vector_get (u, i); + gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0); + } + + /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */ + gsl_matrix *U = gsl_matrix_calloc (n, n); + for (size_t i = 0; i < n; i++) + for (size_t j = 0; j < m; j++) + gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j)); + + /* two dot products to obtain pseudoinverse */ + tmp_mat = gsl_matrix_alloc (m, n); + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat); + + gsl_matrix *A_pinv; + if (swap) + { + A_pinv = gsl_matrix_alloc (n, m); + gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat, 0., A_pinv); + } + else + { + A_pinv = gsl_matrix_alloc (m, n); + gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat, U, 0., A_pinv); + } + + gsl_matrix_free (tmp_mat); + gsl_matrix_free (U); + gsl_matrix_free (Sigma_pinv); + gsl_vector_free (u); + gsl_matrix_free (V); + + return A_pinv; +} + +struct grade + { + size_t y, x; + double value; + }; + +static int +grade_compare_3way (const void *a_, const void *b_) +{ + const struct grade *a = a_; + const struct grade *b = b_; + + return (a->value < b->value ? -1 + : a->value > b->value ? 1 + : a->y < b->y ? -1 + : a->y > b->y ? 1 + : a->x < b->x ? -1 + : a->x > b->x); +} + +static gsl_matrix * +matrix_eval_GRADE (gsl_matrix *m) +{ + size_t n = m->size1 * m->size2; + struct grade *grades = xmalloc (n * sizeof *grades); + + size_t i = 0; + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + grades[i++] = (struct grade) { .y = y, .x = x, .value = *d }; + qsort (grades, n, sizeof *grades, grade_compare_3way); + + for (size_t i = 0; i < n; i++) + gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1); + + free (grades); + + return m; +} + +static double +dot (gsl_vector *a, gsl_vector *b) +{ + double result = 0.0; + for (size_t i = 0; i < a->size; i++) + result += gsl_vector_get (a, i) * gsl_vector_get (b, i); + return result; +} + +static double +norm2 (gsl_vector *v) +{ + double result = 0.0; + for (size_t i = 0; i < v->size; i++) + result += pow2 (gsl_vector_get (v, i)); + return result; +} + +static double +norm (gsl_vector *v) +{ + return sqrt (norm2 (v)); +} + +static gsl_matrix * +matrix_eval_GSCH (gsl_matrix *v) +{ + if (v->size2 < v->size1) + { + msg (SE, _("GSCH requires its argument to have at least as many columns " + "as rows, but it has dimensions (%zu,%zu)."), + v->size1, v->size2); + return NULL; + } + if (!v->size1 || !v->size2) + return v; + + gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2); + size_t ux = 0; + for (size_t vx = 0; vx < v->size2; vx++) + { + gsl_vector u_i = gsl_matrix_column (u, ux).vector; + gsl_vector v_i = gsl_matrix_column (v, vx).vector; + + gsl_vector_memcpy (&u_i, &v_i); + for (size_t j = 0; j < ux; j++) + { + gsl_vector u_j = gsl_matrix_column (u, j).vector; + double scale = dot (&u_j, &u_i) / norm2 (&u_j); + for (size_t k = 0; k < u_i.size; k++) + gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k) + - scale * gsl_vector_get (&u_j, k))); + } + + double len = norm (&u_i); + if (len > 1e-15) + { + gsl_vector_scale (&u_i, 1.0 / len); + if (++ux >= v->size1) + break; + } + } + + if (ux < v->size1) + { + msg (SE, _("Argument to GSCH with dimensions (%zu,%zu) contains only " + "%zu linearly independent columns."), + v->size1, v->size2, ux); + return NULL; + } + + u->size2 = v->size1; + return u; +} + +static gsl_matrix * +matrix_eval_IDENT (double s1, double s2) +{ + if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX) + { + msg (SE, _("Arguments to IDENT must be non-negative integers.")); + return NULL; + } + + gsl_matrix *m = gsl_matrix_alloc (s1, s2); + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = x == y; + return m; +} + +static void +invert_matrix (gsl_matrix *x) +{ + gsl_permutation *p = gsl_permutation_alloc (x->size1); + int signum; + gsl_linalg_LU_decomp (x, p, &signum); + gsl_linalg_LU_invx (x, p); + gsl_permutation_free (p); +} + +static gsl_matrix * +matrix_eval_INV (gsl_matrix *m) +{ + invert_matrix (m); + return m; +} + +static gsl_matrix * +matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b) +{ + gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1, + a->size2 * b->size2); + size_t y = 0; + for (size_t ar = 0; ar < a->size1; ar++) + for (size_t br = 0; br < b->size1; br++, y++) + { + size_t x = 0; + for (size_t ac = 0; ac < a->size2; ac++) + for (size_t bc = 0; bc < b->size2; bc++, x++) + { + double av = gsl_matrix_get (a, ar, ac); + double bv = gsl_matrix_get (b, br, bc); + gsl_matrix_set (k, y, x, av * bv); + } + } + return k; +} + +static gsl_matrix * +matrix_eval_LG10 (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = log10 (*d); + return m; +} + +static gsl_matrix * +matrix_eval_LN (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = log (*d); + return m; +} + +static void +matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n) +{ + /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */ + size_t y = 0; + size_t x = n / 2; + for (size_t i = 1; i <= n * n; i++) + { + gsl_matrix_set (m, y, x, i); + + size_t y1 = !y ? n - 1 : y - 1; + size_t x1 = x + 1 >= n ? 0 : x + 1; + if (gsl_matrix_get (m, y1, x1) == 0) + { + y = y1; + x = x1; + } + else + y = y + 1 >= n ? 0 : y + 1; + } +} + +static void +magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2) +{ + double a = gsl_matrix_get (m, y1, x1); + double b = gsl_matrix_get (m, y2, x2); + gsl_matrix_set (m, y1, x1, b); + gsl_matrix_set (m, y2, x2, a); +} + +static void +matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n) +{ + size_t x, y; + + /* A. Umar, "On the Construction of Even Order Magic Squares", + https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */ + x = y = 0; + for (size_t i = 1; i <= n * n / 2; i++) + { + gsl_matrix_set (m, y, x, i); + if (++y >= n) + { + y = 0; + x++; + } + } + + x = n - 1; + y = 0; + for (size_t i = n * n; i > n * n / 2; i--) + { + gsl_matrix_set (m, y, x, i); + if (++y >= n) + { + y = 0; + x--; + } + } + + for (size_t y = 0; y < n; y++) + for (size_t x = 0; x < n / 2; x++) + { + unsigned int d = gsl_matrix_get (m, y, x); + if (d % 2 != (y < n / 2)) + magic_exchange (m, y, x, y, n - x - 1); + } + + size_t y1 = n / 2; + size_t y2 = n - 1; + size_t x1 = n / 2 - 1; + size_t x2 = n / 2; + magic_exchange (m, y1, x1, y2, x1); + magic_exchange (m, y1, x2, y2, x2); +} + +static void +matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n) +{ + /* A. Umar, "On the Construction of Even Order Magic Squares", + https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */ + size_t x, y; + + x = y = 0; + for (size_t i = 1; ; i++) + { + gsl_matrix_set (m, y, x, i); + if (++y == n / 2 - 1) + y += 2; + else if (y >= n) + { + y = 0; + if (++x >= n / 2) + break; + } + } + + x = n - 1; + y = 0; + for (size_t i = n * n; ; i--) + { + gsl_matrix_set (m, y, x, i); + if (++y == n / 2 - 1) + y += 2; + else if (y >= n) + { + y = 0; + if (--x < n / 2) + break; + } + } + for (size_t y = 0; y < n; y++) + if (y != n / 2 - 1 && y != n / 2) + for (size_t x = 0; x < n / 2; x++) + { + unsigned int d = gsl_matrix_get (m, y, x); + if (d % 2 != (y < n / 2)) + magic_exchange (m, y, x, y, n - x - 1); + } + + size_t a0 = (n * n - 2 * n) / 2 + 1; + for (size_t i = 0; i < n / 2; i++) + { + size_t a = a0 + i; + gsl_matrix_set (m, n / 2, i, a); + gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a); + } + for (size_t i = 0; i < n / 2; i++) + { + size_t a = a0 + i + n / 2; + gsl_matrix_set (m, n / 2 - 1, n - i - 1, a); + gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a); + } + for (size_t x = 1; x < n / 2; x += 2) + magic_exchange (m, n / 2, x, n / 2 - 1, x); + for (size_t x = n / 2 + 2; x <= n - 3; x += 2) + magic_exchange (m, n / 2, x, n / 2 - 1, x); + size_t x1 = n / 2 - 2; + size_t x2 = n / 2 + 1; + size_t y1 = n / 2 - 2; + size_t y2 = n / 2 + 1; + magic_exchange (m, y1, x1, y2, x1); + magic_exchange (m, y1, x2, y2, x2); +} + +static gsl_matrix * +matrix_eval_MAGIC (double n_) +{ + if (n_ < 3 || n_ >= sqrt (SIZE_MAX)) + { + msg (SE, _("MAGIC argument must be an integer 3 or greater.")); + return NULL; + } + size_t n = n_; + + gsl_matrix *m = gsl_matrix_calloc (n, n); + if (n % 2) + matrix_eval_MAGIC_odd (m, n); + else if (n % 4) + matrix_eval_MAGIC_singly_even (m, n); + else + matrix_eval_MAGIC_doubly_even (m, n); + return m; +} + +static gsl_matrix * +matrix_eval_MAKE (double r, double c, double value) +{ + if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX) + { + msg (SE, _("Size arguments to MAKE must be integers.")); + return NULL; + } + + gsl_matrix *m = gsl_matrix_alloc (r, c); + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = value; + return m; +} + +static gsl_matrix * +matrix_eval_MDIAG (gsl_vector *v) +{ + gsl_matrix *m = gsl_matrix_alloc (v->size, v->size); + gsl_vector diagonal = gsl_matrix_diagonal (m).vector; + gsl_vector_memcpy (&diagonal, v); + return m; +} + +static double +matrix_eval_MMAX (gsl_matrix *m) +{ + return gsl_matrix_max (m); +} + +static double +matrix_eval_MMIN (gsl_matrix *m) +{ + return gsl_matrix_min (m); +} + +static gsl_matrix * +matrix_eval_MOD (gsl_matrix *m, double divisor) +{ + if (divisor == 0.0) + { + msg (SE, _("Divisor argument to MOD function must be nonzero.")); + return NULL; + } + + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = fmod (*d, divisor); + return m; +} + +static double +matrix_eval_MSSQ (gsl_matrix *m) +{ + double mssq = 0.0; + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + mssq += *d * *d; + return mssq; +} + +static double +matrix_eval_MSUM (gsl_matrix *m) +{ + double msum = 0.0; + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + msum += *d; + return msum; +} + +static double +matrix_eval_NCOL (gsl_matrix *m) +{ + return m->size2; +} + +static double +matrix_eval_NROW (gsl_matrix *m) +{ + return m->size1; +} + +static double +matrix_eval_RANK (gsl_matrix *m) +{ + gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2)); + gsl_linalg_QR_decomp (m, tau); + gsl_vector_free (tau); + + return gsl_linalg_QRPT_rank (m, -1); +} + +static gsl_matrix * +matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_) +{ + if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX) + { + msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers.")); + return NULL; + } + size_t r = r_; + size_t c = c_; + if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2) + { + msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from " + "product of matrix dimensions.")); + return NULL; + } + + gsl_matrix *dst = gsl_matrix_alloc (r, c); + size_t y1 = 0; + size_t x1 = 0; + MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m) + { + gsl_matrix_set (dst, y1, x1, *d); + if (++x1 >= c) + { + x1 = 0; + y1++; + } + } + return dst; +} + +static gsl_matrix * +matrix_eval_row_extremum (gsl_matrix *m, bool min) +{ + if (m->size2 <= 1) + return m; + else if (!m->size1) + return gsl_matrix_alloc (0, 1); + + gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1); + for (size_t y = 0; y < m->size1; y++) + { + double ext = gsl_matrix_get (m, y, 0); + for (size_t x = 1; x < m->size2; x++) + { + double value = gsl_matrix_get (m, y, x); + if (min ? value < ext : value > ext) + ext = value; + } + gsl_matrix_set (rext, y, 0, ext); + } + return rext; +} + +static gsl_matrix * +matrix_eval_RMAX (gsl_matrix *m) +{ + return matrix_eval_row_extremum (m, false); +} + +static gsl_matrix * +matrix_eval_RMIN (gsl_matrix *m) +{ + return matrix_eval_row_extremum (m, true); +} + +static gsl_matrix * +matrix_eval_RND (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = rint (*d); + return m; +} + +struct rank + { + size_t y, x; + double value; + }; + +static int +rank_compare_3way (const void *a_, const void *b_) +{ + const struct rank *a = a_; + const struct rank *b = b_; + + return a->value < b->value ? -1 : a->value > b->value; +} + +static gsl_matrix * +matrix_eval_RNKORDER (gsl_matrix *m) +{ + size_t n = m->size1 * m->size2; + struct rank *ranks = xmalloc (n * sizeof *ranks); + size_t i = 0; + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d }; + qsort (ranks, n, sizeof *ranks, rank_compare_3way); + + for (size_t i = 0; i < n; ) + { + size_t j; + for (j = i + 1; j < n; j++) + if (ranks[i].value != ranks[j].value) + break; + + double rank = (i + j + 1.0) / 2.0; + for (size_t k = i; k < j; k++) + gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank); + + i = j; + } + + free (ranks); + + return m; +} + +static gsl_matrix * +matrix_eval_row_sum (gsl_matrix *m, bool square) +{ + if (m->size1 == 0) + return m; + else if (!m->size1) + return gsl_matrix_alloc (0, 1); + + gsl_matrix *result = gsl_matrix_alloc (m->size1, 1); + for (size_t y = 0; y < m->size1; y++) + { + double sum = 0; + for (size_t x = 0; x < m->size2; x++) + { + double d = gsl_matrix_get (m, y, x); + sum += square ? pow2 (d) : d; + } + gsl_matrix_set (result, y, 0, sum); + } + return result; +} + +static gsl_matrix * +matrix_eval_RSSQ (gsl_matrix *m) +{ + return matrix_eval_row_sum (m, true); +} + +static gsl_matrix * +matrix_eval_RSUM (gsl_matrix *m) +{ + return matrix_eval_row_sum (m, false); +} + +static gsl_matrix * +matrix_eval_SIN (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = sin (*d); + return m; +} + +static gsl_matrix * +matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2) +{ + if (m1->size1 != m2->size1) + { + msg (SE, _("SOLVE requires its arguments to have the same number of " + "rows, but the first argument has dimensions (%zu,%zu) and " + "the second (%zu,%zu)."), + m1->size1, m1->size2, + m2->size1, m2->size2); + return NULL; + } + + gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2); + gsl_permutation *p = gsl_permutation_alloc (m1->size1); + int signum; + gsl_linalg_LU_decomp (m1, p, &signum); + for (size_t i = 0; i < m2->size2; i++) + { + gsl_vector bi = gsl_matrix_column (m2, i).vector; + gsl_vector xi = gsl_matrix_column (x, i).vector; + gsl_linalg_LU_solve (m1, p, &bi, &xi); + } + gsl_permutation_free (p); + return x; +} + +static gsl_matrix * +matrix_eval_SQRT (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + { + if (*d < 0) + { + msg (SE, _("Argument to SQRT must be nonnegative.")); + return NULL; + } + *d = sqrt (*d); + } + return m; +} + +static gsl_matrix * +matrix_eval_SSCP (gsl_matrix *m) +{ + gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2); + gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp); + return sscp; +} + +static gsl_matrix * +matrix_eval_SVAL (gsl_matrix *m) +{ + gsl_matrix *tmp_mat = NULL; + if (m->size2 > m->size1) + { + tmp_mat = gsl_matrix_alloc (m->size2, m->size1); + gsl_matrix_transpose_memcpy (tmp_mat, m); + m = tmp_mat; + } + + /* Do SVD. */ + gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2); + gsl_vector *S = gsl_vector_alloc (m->size2); + gsl_vector *work = gsl_vector_alloc (m->size2); + gsl_linalg_SV_decomp (m, V, S, work); + + gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1); + for (size_t i = 0; i < m->size2; i++) + gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i)); + + gsl_matrix_free (V); + gsl_vector_free (S); + gsl_vector_free (work); + gsl_matrix_free (tmp_mat); + + return vals; +} + +static gsl_matrix * +matrix_eval_SWEEP (gsl_matrix *m, double d) +{ + if (d < 1 || d > SIZE_MAX) + { + msg (SE, _("Scalar argument to SWEEP must be integer.")); + return NULL; + } + size_t k = d - 1; + if (k >= MIN (m->size1, m->size2)) + { + msg (SE, _("Scalar argument to SWEEP must be integer less than or " + "equal to the smaller of the matrix argument's rows and " + "columns.")); + return NULL; + } + + double m_kk = gsl_matrix_get (m, k, k); + if (fabs (m_kk) > 1e-19) + { + gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2); + MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a) + { + double m_ij = gsl_matrix_get (m, i, j); + double m_ik = gsl_matrix_get (m, i, k); + double m_kj = gsl_matrix_get (m, k, j); + *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj + : i != k ? -m_ik + : j != k ? m_kj + : 1.0) / m_kk; + } + return a; + } + else + { + for (size_t i = 0; i < m->size1; i++) + { + gsl_matrix_set (m, i, k, 0); + gsl_matrix_set (m, k, i, 0); + } + return m; + } +} + +static double +matrix_eval_TRACE (gsl_matrix *m) +{ + double sum = 0; + size_t n = MIN (m->size1, m->size2); + for (size_t i = 0; i < n; i++) + sum += gsl_matrix_get (m, i, i); + return sum; +} + +static gsl_matrix * +matrix_eval_T (gsl_matrix *m) +{ + return matrix_eval_TRANSPOS (m); +} + +static gsl_matrix * +matrix_eval_TRANSPOS (gsl_matrix *m) +{ + if (m->size1 == m->size2) + { + gsl_matrix_transpose (m); + return m; + } + else + { + gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1); + gsl_matrix_transpose_memcpy (t, m); + return t; + } +} + +static gsl_matrix * +matrix_eval_TRUNC (gsl_matrix *m) +{ + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = trunc (*d); + return m; +} + +static gsl_matrix * +matrix_eval_UNIFORM (double r_, double c_) +{ + if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX) + { + msg (SE, _("Arguments to UNIFORM must be integers.")); + return NULL; + } + size_t r = r_; + size_t c = c_; + if (size_overflow_p (xtimes (r, xmax (c, 1)))) + { + msg (SE, _("Product of arguments to UNIFORM exceeds memory size.")); + return NULL; + } + + gsl_matrix *m = gsl_matrix_alloc (r, c); + MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) + *d = gsl_ran_flat (get_rng (), 0, 1); + return m; +} + +struct matrix_function + { + const char *name; + enum matrix_op op; + size_t min_args, max_args; + }; + +static struct matrix_expr *matrix_parse_expr (struct matrix_state *); + +static const struct matrix_function * +matrix_parse_function_name (const char *token) +{ + static const struct matrix_function functions[] = + { +#define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS }, + MATRIX_FUNCTIONS +#undef F + }; + enum { N_FUNCTIONS = sizeof functions / sizeof *functions }; + + for (size_t i = 0; i < N_FUNCTIONS; i++) + { + if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3)) + return &functions[i]; + } + return NULL; +} + +static struct read_file * +read_file_create (struct matrix_state *s, struct file_handle *fh) +{ + for (size_t i = 0; i < s->n_read_files; i++) + { + struct read_file *rf = s->read_files[i]; + if (rf->file == fh) + { + fh_unref (fh); + return rf; + } + } + + struct read_file *rf = xmalloc (sizeof *rf); + *rf = (struct read_file) { .file = fh }; + + s->read_files = xrealloc (s->read_files, + (s->n_read_files + 1) * sizeof *s->read_files); + s->read_files[s->n_read_files++] = rf; + return rf; +} + +static struct dfm_reader * +read_file_open (struct read_file *rf) +{ + if (!rf->reader) + rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding); + return rf->reader; +} + +static bool +matrix_parse_function (struct matrix_state *s, const char *token, + struct matrix_expr **exprp) +{ + *exprp = NULL; + if (lex_next_token (s->lexer, 1) != T_LPAREN) + return false; + + if (lex_match_id (s->lexer, "EOF")) + { + lex_get (s->lexer); + struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session); + if (!fh) + return true; + + if (!lex_force_match (s->lexer, T_RPAREN)) + { + fh_unref (fh); + return true; + } + + struct read_file *rf = read_file_create (s, fh); + + struct matrix_expr *e = xmalloc (sizeof *e); + *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf }; + *exprp = e; + return true; + } + + const struct matrix_function *f = matrix_parse_function_name (token); + if (!f) + return false; + + lex_get_n (s->lexer, 2); + + struct matrix_expr *e = xmalloc (sizeof *e); + *e = (struct matrix_expr) { .op = f->op, .subs = NULL }; + + size_t allocated_subs = 0; + do + { + struct matrix_expr *sub = matrix_parse_expr (s); + if (!sub) + goto error; + + if (e->n_subs >= allocated_subs) + e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs); + e->subs[e->n_subs++] = sub; + } + while (lex_match (s->lexer, T_COMMA)); + if (!lex_force_match (s->lexer, T_RPAREN)) + goto error; + + *exprp = e; + return true; + +error: + matrix_expr_destroy (e); + return true; +} + +static struct matrix_expr * +matrix_parse_primary (struct matrix_state *s) +{ + if (lex_is_number (s->lexer)) + { + double number = lex_number (s->lexer); + lex_get (s->lexer); + return matrix_expr_create_number (number); + } + else if (lex_is_string (s->lexer)) + { + char string[sizeof (double)]; + buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' '); + lex_get (s->lexer); + + double number; + memcpy (&number, string, sizeof number); + return matrix_expr_create_number (number); + } + else if (lex_match (s->lexer, T_LPAREN)) + { + struct matrix_expr *e = matrix_parse_or_xor (s); + if (!e || !lex_force_match (s->lexer, T_RPAREN)) + { + matrix_expr_destroy (e); + return NULL; + } + return e; + } + else if (lex_match (s->lexer, T_LCURLY)) + { + struct matrix_expr *e = matrix_parse_curly_semi (s); + if (!e || !lex_force_match (s->lexer, T_RCURLY)) + { + matrix_expr_destroy (e); + return NULL; + } + return e; + } + else if (lex_token (s->lexer) == T_ID) + { + struct matrix_expr *retval; + if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval)) + return retval; + + struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer)); + if (!var) + { + lex_error (s->lexer, _("Unknown variable %s."), + lex_tokcstr (s->lexer)); + return NULL; + } + lex_get (s->lexer); + + struct matrix_expr *e = xmalloc (sizeof *e); + *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var }; + return e; + } + else if (lex_token (s->lexer) == T_ALL) + { + struct matrix_expr *retval; + if (matrix_parse_function (s, "ALL", &retval)) + return retval; + } + + lex_error (s->lexer, NULL); + return NULL; +} + +static struct matrix_expr *matrix_parse_postfix (struct matrix_state *); + +static bool +matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp) +{ + if (lex_match (s->lexer, T_COLON)) + { + *indexp = NULL; + return true; + } + else + { + *indexp = matrix_parse_expr (s); + return *indexp != NULL; + } +} + +static struct matrix_expr * +matrix_parse_postfix (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_primary (s); + if (!lhs || !lex_match (s->lexer, T_LPAREN)) + return lhs; + + struct matrix_expr *i0; + if (!matrix_parse_index_expr (s, &i0)) + { + matrix_expr_destroy (lhs); + return NULL; + } + if (lex_match (s->lexer, T_RPAREN)) + return (i0 + ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0) + : matrix_expr_create_1 (MOP_VEC_ALL, lhs)); + else if (lex_match (s->lexer, T_COMMA)) + { + struct matrix_expr *i1; + if (!matrix_parse_index_expr (s, &i1) + || !lex_force_match (s->lexer, T_RPAREN)) + { + matrix_expr_destroy (lhs); + matrix_expr_destroy (i0); + matrix_expr_destroy (i1); + return NULL; + } + return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1) + : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0) + : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1) + : lhs); + } + else + { + lex_error_expecting (s->lexer, "`)'", "`,'"); + return NULL; + } +} + +static struct matrix_expr * +matrix_parse_unary (struct matrix_state *s) +{ + if (lex_match (s->lexer, T_DASH)) + { + struct matrix_expr *lhs = matrix_parse_unary (s); + return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL; + } + else if (lex_match (s->lexer, T_PLUS)) + return matrix_parse_unary (s); + else + return matrix_parse_postfix (s); +} + +static struct matrix_expr * +matrix_parse_seq (struct matrix_state *s) +{ + struct matrix_expr *start = matrix_parse_unary (s); + if (!start || !lex_match (s->lexer, T_COLON)) + return start; + + struct matrix_expr *end = matrix_parse_unary (s); + if (!end) + { + matrix_expr_destroy (start); + return NULL; + } + + if (lex_match (s->lexer, T_COLON)) + { + struct matrix_expr *increment = matrix_parse_unary (s); + if (!increment) + { + matrix_expr_destroy (start); + matrix_expr_destroy (end); + return NULL; + } + return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment); + } + else + return matrix_expr_create_2 (MOP_SEQ, start, end); +} + +static struct matrix_expr * +matrix_parse_exp (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_seq (s); + if (!lhs) + return NULL; + + for (;;) + { + enum matrix_op op; + if (lex_match (s->lexer, T_EXP)) + op = MOP_EXP_MAT; + else if (lex_match_phrase (s->lexer, "&**")) + op = MOP_EXP_ELEMS; + else + return lhs; + + struct matrix_expr *rhs = matrix_parse_seq (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (op, lhs, rhs); + } +} + +static struct matrix_expr * +matrix_parse_mul_div (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_exp (s); + if (!lhs) + return NULL; + + for (;;) + { + enum matrix_op op; + if (lex_match (s->lexer, T_ASTERISK)) + op = MOP_MUL_MAT; + else if (lex_match (s->lexer, T_SLASH)) + op = MOP_DIV_ELEMS; + else if (lex_match_phrase (s->lexer, "&*")) + op = MOP_MUL_ELEMS; + else if (lex_match_phrase (s->lexer, "&/")) + op = MOP_DIV_ELEMS; + else + return lhs; + + struct matrix_expr *rhs = matrix_parse_exp (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (op, lhs, rhs); + } +} + +static struct matrix_expr * +matrix_parse_add_sub (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_mul_div (s); + if (!lhs) + return NULL; + + for (;;) + { + enum matrix_op op; + if (lex_match (s->lexer, T_PLUS)) + op = MOP_ADD_ELEMS; + else if (lex_match (s->lexer, T_DASH)) + op = MOP_SUB_ELEMS; + else if (lex_token (s->lexer) == T_NEG_NUM) + op = MOP_ADD_ELEMS; + else + return lhs; + + struct matrix_expr *rhs = matrix_parse_mul_div (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (op, lhs, rhs); + } +} + +static struct matrix_expr * +matrix_parse_relational (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_add_sub (s); + if (!lhs) + return NULL; + + for (;;) + { + enum matrix_op op; + if (lex_match (s->lexer, T_GT)) + op = MOP_GT; + else if (lex_match (s->lexer, T_GE)) + op = MOP_GE; + else if (lex_match (s->lexer, T_LT)) + op = MOP_LT; + else if (lex_match (s->lexer, T_LE)) + op = MOP_LE; + else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ)) + op = MOP_EQ; + else if (lex_match (s->lexer, T_NE)) + op = MOP_NE; + else + return lhs; + + struct matrix_expr *rhs = matrix_parse_add_sub (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (op, lhs, rhs); + } +} + +static struct matrix_expr * +matrix_parse_not (struct matrix_state *s) +{ + if (lex_match (s->lexer, T_NOT)) + { + struct matrix_expr *lhs = matrix_parse_not (s); + return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL; + } + else + return matrix_parse_relational (s); +} + +static struct matrix_expr * +matrix_parse_and (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_not (s); + if (!lhs) + return NULL; + + while (lex_match (s->lexer, T_AND)) + { + struct matrix_expr *rhs = matrix_parse_not (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs); + } + return lhs; +} + +static struct matrix_expr * +matrix_parse_or_xor (struct matrix_state *s) +{ + struct matrix_expr *lhs = matrix_parse_and (s); + if (!lhs) + return NULL; + + for (;;) + { + enum matrix_op op; + if (lex_match (s->lexer, T_OR)) + op = MOP_OR; + else if (lex_match_id (s->lexer, "XOR")) + op = MOP_XOR; + else + return lhs; + + struct matrix_expr *rhs = matrix_parse_and (s); + if (!rhs) + { + matrix_expr_destroy (lhs); + return NULL; + } + lhs = matrix_expr_create_2 (op, lhs, rhs); + } +} + +static struct matrix_expr * +matrix_parse_expr (struct matrix_state *s) +{ + return matrix_parse_or_xor (s); +} + +/* Expression evaluation. */ + +static double +matrix_op_eval (enum matrix_op op, double a, double b) +{ + switch (op) + { + case MOP_ADD_ELEMS: return a + b; + case MOP_SUB_ELEMS: return a - b; + case MOP_MUL_ELEMS: return a * b; + case MOP_DIV_ELEMS: return a / b; + case MOP_EXP_ELEMS: return pow (a, b); + case MOP_GT: return a > b; + case MOP_GE: return a >= b; + case MOP_LT: return a < b; + case MOP_LE: return a <= b; + case MOP_EQ: return a == b; + case MOP_NE: return a != b; + case MOP_AND: return (a > 0) && (b > 0); + case MOP_OR: return (a > 0) || (b > 0); + case MOP_XOR: return (a > 0) != (b > 0); + +#define F(NAME, PROTOTYPE) case MOP_F_##NAME: + MATRIX_FUNCTIONS +#undef F + case MOP_NEGATE: + case MOP_SEQ: + case MOP_SEQ_BY: + case MOP_MUL_MAT: + case MOP_EXP_MAT: + case MOP_NOT: + case MOP_PASTE_HORZ: + case MOP_PASTE_VERT: + case MOP_EMPTY: + case MOP_VEC_INDEX: + case MOP_VEC_ALL: + case MOP_MAT_INDEX: + case MOP_ROW_INDEX: + case MOP_COL_INDEX: + case MOP_NUMBER: + case MOP_VARIABLE: + case MOP_EOF: + NOT_REACHED (); + } + NOT_REACHED (); +} + +static const char * +matrix_op_name (enum matrix_op op) +{ + switch (op) + { + case MOP_ADD_ELEMS: return "+"; + case MOP_SUB_ELEMS: return "-"; + case MOP_MUL_ELEMS: return "&*"; + case MOP_DIV_ELEMS: return "&/"; + case MOP_EXP_ELEMS: return "&**"; + case MOP_GT: return ">"; + case MOP_GE: return ">="; + case MOP_LT: return "<"; + case MOP_LE: return "<="; + case MOP_EQ: return "="; + case MOP_NE: return "<>"; + case MOP_AND: return "AND"; + case MOP_OR: return "OR"; + case MOP_XOR: return "XOR"; + +#define F(NAME, PROTOTYPE) case MOP_F_##NAME: + MATRIX_FUNCTIONS +#undef F + case MOP_NEGATE: + case MOP_SEQ: + case MOP_SEQ_BY: + case MOP_MUL_MAT: + case MOP_EXP_MAT: + case MOP_NOT: + case MOP_PASTE_HORZ: + case MOP_PASTE_VERT: + case MOP_EMPTY: + case MOP_VEC_INDEX: + case MOP_VEC_ALL: + case MOP_MAT_INDEX: + case MOP_ROW_INDEX: + case MOP_COL_INDEX: + case MOP_NUMBER: + case MOP_VARIABLE: + case MOP_EOF: + NOT_REACHED (); + } + NOT_REACHED (); +} + +static bool +is_scalar (const gsl_matrix *m) +{ + return m->size1 == 1 && m->size2 == 1; +} + +static double +to_scalar (const gsl_matrix *m) +{ + assert (is_scalar (m)); + return gsl_matrix_get (m, 0, 0); +} + +static gsl_matrix * +matrix_expr_evaluate_elementwise (enum matrix_op op, + gsl_matrix *a, gsl_matrix *b) +{ + if (is_scalar (b)) + { + double be = to_scalar (b); + for (size_t r = 0; r < a->size1; r++) + for (size_t c = 0; c < a->size2; c++) + { + double *ae = gsl_matrix_ptr (a, r, c); + *ae = matrix_op_eval (op, *ae, be); + } + return a; + } + else if (is_scalar (a)) + { + double ae = to_scalar (a); + for (size_t r = 0; r < b->size1; r++) + for (size_t c = 0; c < b->size2; c++) + { + double *be = gsl_matrix_ptr (b, r, c); + *be = matrix_op_eval (op, ae, *be); + } + return b; + } + else if (a->size1 == b->size1 && a->size2 == b->size2) + { + for (size_t r = 0; r < a->size1; r++) + for (size_t c = 0; c < a->size2; c++) + { + double *ae = gsl_matrix_ptr (a, r, c); + double be = gsl_matrix_get (b, r, c); + *ae = matrix_op_eval (op, *ae, be); + } + return a; + } + else + { + msg (SE, _("Operands to %s must have the same dimensions or one " + "must be a scalar, not matrices with dimensions (%zu,%zu) " + "and (%zu,%zu)."), + matrix_op_name (op), a->size1, a->size2, b->size1, b->size2); + return NULL; + } +} + +static gsl_matrix * +matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b) +{ + if (is_scalar (a) || is_scalar (b)) + return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b); + + if (a->size2 != b->size1) + { + msg (SE, _("Matrices with dimensions (%zu,%zu) and (%zu,%zu) are " + "not conformable for multiplication."), + a->size1, a->size2, b->size1, b->size2); + return NULL; + } + + gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2); + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c); + return c; +} + +static void +swap_matrix (gsl_matrix **a, gsl_matrix **b) +{ + gsl_matrix *tmp = *a; + *a = *b; + *b = tmp; +} + +static void +mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y, + gsl_matrix **tmp) +{ + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp); + swap_matrix (z, tmp); +} + +static void +square_matrix (gsl_matrix **x, gsl_matrix **tmp) +{ + mul_matrix (x, *x, *x, tmp); +} + +static gsl_matrix * +matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b) +{ + gsl_matrix *x = x_; + if (x->size1 != x->size2) + { + msg (SE, _("Matrix exponentation with ** requires a square matrix on " + "the left-hand size, not one with dimensions (%zu,%zu)."), + x->size1, x->size2); + return NULL; + } + if (!is_scalar (b)) + { + msg (SE, _("Matrix exponentiation with ** requires a scalar on the " + "right-hand side, not a matrix with dimensions (%zu,%zu)."), + b->size1, b->size2); + return NULL; + } + double bf = to_scalar (b); + if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX) + { + msg (SE, _("Exponent %.1f in matrix multiplication is non-integer " + "or outside the valid range."), bf); + return NULL; + } + long int bl = bf; + + gsl_matrix *tmp = gsl_matrix_alloc (x->size1, x->size2); + gsl_matrix *y = gsl_matrix_alloc (x->size1, x->size2); + gsl_matrix_set_identity (y); + if (bl == 0) + return y; + + for (unsigned long int n = labs (bl); n > 1; n /= 2) + if (n & 1) + { + mul_matrix (&y, x, y, &tmp); + square_matrix (&x, &tmp); + } + else + square_matrix (&x, &tmp); + + mul_matrix (&y, x, y, &tmp); + if (bf < 0) + invert_matrix (y); + + if (tmp != x_) + gsl_matrix_free (tmp); + return y; +} + +static gsl_matrix * +matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_, + gsl_matrix *by_) +{ + if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_))) + { + msg (SE, _("All operands of : operator must be scalars.")); + return NULL; + } + + long int start = to_scalar (start_); + long int end = to_scalar (end_); + long int by = by_ ? to_scalar (by_) : 1; + + if (!by) + { + msg (SE, _("The increment operand to : must be nonzero.")); + return NULL; + } + + long int n = (end > start && by > 0 ? (end - start + by) / by + : end < start && by < 0 ? (start - end - by) / -by + : 0); + gsl_matrix *m = gsl_matrix_alloc (1, n); + for (long int i = 0; i < n; i++) + gsl_matrix_set (m, 0, i, start + i * by); + return m; +} + +static gsl_matrix * +matrix_expr_evaluate_not (gsl_matrix *a) +{ + for (size_t r = 0; r < a->size1; r++) + for (size_t c = 0; c < a->size2; c++) + { + double *ae = gsl_matrix_ptr (a, r, c); + *ae = !(*ae > 0); + } + return a; +} + +static gsl_matrix * +matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b) +{ + if (a->size1 != b->size1) + { + if (!a->size1 || !a->size2) + return b; + else if (!b->size1 || !b->size2) + return a; + + msg (SE, _("All columns in a matrix must have the same number of rows, " + "but this tries to paste matrices with %zu and %zu rows."), + a->size1, b->size1); + return NULL; + } + + gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2); + for (size_t y = 0; y < a->size1; y++) + { + for (size_t x = 0; x < a->size2; x++) + gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x)); + for (size_t x = 0; x < b->size2; x++) + gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x)); + } + return c; +} + +static gsl_matrix * +matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b) +{ + if (a->size2 != b->size2) + { + if (!a->size1 || !a->size2) + return b; + else if (!b->size1 || !b->size2) + return a; + + msg (SE, _("All rows in a matrix must have the same number of columns, " + "but this tries to stack matrices with %zu and %zu columns."), + a->size2, b->size2); + return NULL; + } + + gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2); + for (size_t x = 0; x < a->size2; x++) + { + for (size_t y = 0; y < a->size1; y++) + gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x)); + for (size_t y = 0; y < b->size1; y++) + gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x)); + } + return c; +} + +static gsl_vector * +matrix_to_vector (gsl_matrix *m) +{ + assert (m->owner); + gsl_vector v = to_vector (m); + assert (v.block == m->block || !v.block); + assert (!v.owner); + v.owner = 1; + m->owner = 0; + return xmemdup (&v, sizeof v); +} + +enum index_type { + IV_ROW, + IV_COLUMN, + IV_VECTOR +}; + +struct index_vector + { + size_t *indexes; + size_t n; + }; +#define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 } + +static bool +matrix_normalize_index_vector (gsl_matrix *m, size_t size, + enum index_type index_type, size_t other_size, + struct index_vector *iv) +{ + if (m) + { + if (!is_vector (m)) + { + switch (index_type) + { + case IV_VECTOR: + msg (SE, _("Vector index must be scalar or vector, not a " + "matrix with dimensions (%zu,%zu)."), + m->size1, m->size2); + break; + + case IV_ROW: + msg (SE, _("Matrix row index must be scalar or vector, not a " + "matrix with dimensions (%zu,%zu)."), + m->size1, m->size2); + break; + + case IV_COLUMN: + msg (SE, _("Matrix column index must be scalar or vector, not a " + "matrix with dimensions (%zu,%zu)."), + m->size1, m->size2); + break; + } + return false; + } + + gsl_vector v = to_vector (m); + *iv = (struct index_vector) { + .indexes = xnmalloc (v.size, sizeof *iv->indexes), + .n = v.size, + }; + for (size_t i = 0; i < v.size; i++) + { + double index = gsl_vector_get (&v, i); + if (index < 1 || index >= size + 1) + { + switch (index_type) + { + case IV_VECTOR: + msg (SE, _("Index %g is out of range for vector " + "with %zu elements."), index, size); + break; + + case IV_ROW: + msg (SE, _("%g is not a valid row index for a matrix " + "with dimensions (%zu,%zu)."), + index, size, other_size); + break; + + case IV_COLUMN: + msg (SE, _("%g is not a valid column index for a matrix " + "with dimensions (%zu,%zu)."), + index, other_size, size); + break; + } + + free (iv->indexes); + return false; + } + iv->indexes[i] = index - 1; + } + return true; + } + else + { + *iv = (struct index_vector) { + .indexes = xnmalloc (size, sizeof *iv->indexes), + .n = size, + }; + for (size_t i = 0; i < size; i++) + iv->indexes[i] = i; + return true; + } +} + +static gsl_matrix * +matrix_expr_evaluate_vec_all (gsl_matrix *sm) +{ + if (!is_vector (sm)) + { + msg (SE, _("Vector index operator must be applied to vector, " + "not a matrix with dimensions (%zu,%zu)."), + sm->size1, sm->size2); + return NULL; + } + + return sm; +} + +static gsl_matrix * +matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im) +{ + if (!matrix_expr_evaluate_vec_all (sm)) + return NULL; + + gsl_vector sv = to_vector (sm); + struct index_vector iv; + if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv)) + return NULL; + + gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n, + sm->size1 == 1 ? iv.n : 1); + gsl_vector dv = to_vector (dm); + for (size_t dx = 0; dx < iv.n; dx++) + { + size_t sx = iv.indexes[dx]; + gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx)); + } + free (iv.indexes); + + return dm; +} + +static gsl_matrix * +matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0, + gsl_matrix *im1) +{ + struct index_vector iv0; + if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0)) + return NULL; + + struct index_vector iv1; + if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1, + &iv1)) + { + free (iv0.indexes); + return NULL; + } + + gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n); + for (size_t dy = 0; dy < iv0.n; dy++) + { + size_t sy = iv0.indexes[dy]; + + for (size_t dx = 0; dx < iv1.n; dx++) + { + size_t sx = iv1.indexes[dx]; + gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx)); + } + } + free (iv0.indexes); + free (iv1.indexes); + return dm; +} + +#define F(NAME, PROTOTYPE) \ + static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \ + const char *name, gsl_matrix *[], size_t, \ + matrix_proto_##PROTOTYPE *); +MATRIX_FUNCTIONS +#undef F + +static bool +check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index) +{ + if (!is_scalar (subs[index])) + { + msg (SE, _("Function %s argument %zu must be a scalar, " + "but it has dimensions (%zu,%zu)."), + name, index + 1, subs[index]->size1, subs[index]->size2); + return false; + } + return true; +} + +static bool +check_vector_arg (const char *name, gsl_matrix *subs[], size_t index) +{ + if (!is_vector (subs[index])) + { + msg (SE, _("Function %s argument %zu must be a vector, " + "but it has dimensions (%zu,%zu)."), + name, index + 1, subs[index]->size1, subs[index]->size2); + return false; + } + return true; +} + +static bool +to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[]) +{ + for (size_t i = 0; i < n_subs; i++) + { + if (!check_scalar_arg (name, subs, i)) + return false; + d[i] = to_scalar (subs[i]); + } + return true; +} + +static gsl_matrix * +matrix_expr_evaluate_m_d (const char *name, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_d *f) +{ + assert (n_subs == 1); + + double d; + return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL; +} + +static gsl_matrix * +matrix_expr_evaluate_m_dd (const char *name, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_dd *f) +{ + assert (n_subs == 2); + + double d[2]; + return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL; +} + +static gsl_matrix * +matrix_expr_evaluate_m_ddd (const char *name, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_ddd *f) +{ + assert (n_subs == 3); + + double d[3]; + return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL; +} + +static gsl_matrix * +matrix_expr_evaluate_m_m (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_m *f) +{ + assert (n_subs == 1); + return f (subs[0]); +} + +static gsl_matrix * +matrix_expr_evaluate_m_md (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_md *f) +{ + assert (n_subs == 2); + if (!check_scalar_arg (name, subs, 1)) + return NULL; + return f (subs[0], to_scalar (subs[1])); +} + +static gsl_matrix * +matrix_expr_evaluate_m_mdd (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_mdd *f) +{ + assert (n_subs == 3); + if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2)) + return NULL; + return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2])); +} + +static gsl_matrix * +matrix_expr_evaluate_m_mm (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_mm *f) +{ + assert (n_subs == 2); + return f (subs[0], subs[1]); +} + +static gsl_matrix * +matrix_expr_evaluate_m_v (const char *name, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_v *f) +{ + assert (n_subs == 1); + if (!check_vector_arg (name, subs, 0)) + return NULL; + gsl_vector v = to_vector (subs[0]); + return f (&v); +} + +static gsl_matrix * +matrix_expr_evaluate_d_m (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_d_m *f) +{ + assert (n_subs == 1); + + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, f (subs[0])); + return m; +} + +static gsl_matrix * +matrix_expr_evaluate_m_any (const char *name UNUSED, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_any *f) +{ + return f (subs, n_subs); +} + +static gsl_matrix * +matrix_expr_evaluate_IDENT (const char *name, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_IDENT *f) +{ + assert (n_subs <= 2); + + double d[2]; + if (!to_scalar_args (name, subs, n_subs, d)) + return NULL; + + return f (d[0], d[n_subs - 1]); +} + +static gsl_matrix * +matrix_expr_evaluate (const struct matrix_expr *e) +{ + if (e->op == MOP_NUMBER) + { + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, e->number); + return m; + } + else if (e->op == MOP_VARIABLE) + { + const gsl_matrix *src = e->variable->value; + if (!src) + { + msg (SE, _("Uninitialized variable %s used in expression."), + e->variable->name); + return NULL; + } + + gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2); + gsl_matrix_memcpy (dst, src); + return dst; + } + else if (e->op == MOP_EOF) + { + struct dfm_reader *reader = read_file_open (e->eof); + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader)); + return m; + } + + enum { N_LOCAL = 3 }; + gsl_matrix *local_subs[N_LOCAL]; + gsl_matrix **subs = (e->n_subs < N_LOCAL + ? local_subs + : xmalloc (e->n_subs * sizeof *subs)); + + for (size_t i = 0; i < e->n_subs; i++) + { + subs[i] = matrix_expr_evaluate (e->subs[i]); + if (!subs[i]) + { + for (size_t j = 0; j < i; j++) + gsl_matrix_free (subs[j]); + if (subs != local_subs) + free (subs); + return NULL; + } + } + + gsl_matrix *result = NULL; + switch (e->op) + { +#define F(NAME, PROTOTYPE) \ + case MOP_F_##NAME: \ + result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \ + subs, e->n_subs, \ + matrix_eval_##NAME); \ + break; + MATRIX_FUNCTIONS +#undef F + + case MOP_NEGATE: + gsl_matrix_scale (subs[0], -1.0); + result = subs[0]; + break; + + case MOP_ADD_ELEMS: + case MOP_SUB_ELEMS: + case MOP_MUL_ELEMS: + case MOP_DIV_ELEMS: + case MOP_EXP_ELEMS: + case MOP_GT: + case MOP_GE: + case MOP_LT: + case MOP_LE: + case MOP_EQ: + case MOP_NE: + case MOP_AND: + case MOP_OR: + case MOP_XOR: + result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]); + break; + + case MOP_NOT: + result = matrix_expr_evaluate_not (subs[0]); + break; + + case MOP_SEQ: + result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL); + break; + + case MOP_SEQ_BY: + result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]); + break; + + case MOP_MUL_MAT: + result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]); + break; + + case MOP_EXP_MAT: + result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]); + break; + + case MOP_PASTE_HORZ: + result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]); + break; + + case MOP_PASTE_VERT: + result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]); + break; + + case MOP_EMPTY: + result = gsl_matrix_alloc (0, 0); + break; + + case MOP_VEC_INDEX: + result = matrix_expr_evaluate_vec_index (subs[0], subs[1]); + break; + + case MOP_VEC_ALL: + result = matrix_expr_evaluate_vec_all (subs[0]); + break; + + case MOP_MAT_INDEX: + result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]); + break; + + case MOP_ROW_INDEX: + result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL); + break; + + case MOP_COL_INDEX: + result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]); + break; + + case MOP_NUMBER: + case MOP_VARIABLE: + case MOP_EOF: + NOT_REACHED (); + } + + for (size_t i = 0; i < e->n_subs; i++) + if (subs[i] != result) + gsl_matrix_free (subs[i]); + if (subs != local_subs) + free (subs); + return result; +} + +static bool +matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context, + double *d) +{ + gsl_matrix *m = matrix_expr_evaluate (e); + if (!m) + return false; + + if (!is_scalar (m)) + { + msg (SE, _("Expression for %s must evaluate to scalar, " + "not a matrix with dimensions (%zu,%zu)."), + context, m->size1, m->size2); + gsl_matrix_free (m); + return false; + } + + *d = to_scalar (m); + gsl_matrix_free (m); + return true; +} + +static bool +matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context, + long int *integer) +{ + double d; + if (!matrix_expr_evaluate_scalar (e, context, &d)) + return false; + + d = trunc (d); + if (d < LONG_MIN || d > LONG_MAX) + { + msg (SE, _("Expression for %s is outside the integer range."), context); + return false; + } + *integer = d; + return true; +} + +struct matrix_lvalue + { + struct matrix_var *var; + struct matrix_expr *indexes[2]; + size_t n_indexes; + }; + +static void +matrix_lvalue_destroy (struct matrix_lvalue *lvalue) +{ + if (lvalue) + for (size_t i = 0; i < lvalue->n_indexes; i++) + matrix_expr_destroy (lvalue->indexes[i]); +} + +static struct matrix_lvalue * +matrix_lvalue_parse (struct matrix_state *s) +{ + struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue); + if (!lex_force_id (s->lexer)) + goto error; + lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer)); + if (lex_next_token (s->lexer, 1) == T_LPAREN) + { + if (!lvalue->var) + { + msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer)); + return NULL; + } + + lex_get_n (s->lexer, 2); + + if (!matrix_parse_index_expr (s, &lvalue->indexes[0])) + goto error; + lvalue->n_indexes++; + + if (lex_match (s->lexer, T_COMMA)) + { + if (!matrix_parse_index_expr (s, &lvalue->indexes[1])) + goto error; + lvalue->n_indexes++; + } + if (!lex_force_match (s->lexer, T_RPAREN)) + goto error; + } + else + { + if (!lvalue->var) + lvalue->var = matrix_var_create (s, lex_tokss (s->lexer)); + lex_get (s->lexer); + } + return lvalue; + +error: + matrix_lvalue_destroy (lvalue); + return NULL; +} + +static bool +matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size, + enum index_type index_type, size_t other_size, + struct index_vector *iv) +{ + gsl_matrix *m; + if (e) + { + m = matrix_expr_evaluate (e); + if (!m) + return false; + } + else + m = NULL; + + return matrix_normalize_index_vector (m, size, index_type, other_size, iv); +} + +static bool +matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue, + struct index_vector *iv, gsl_matrix *sm) +{ + gsl_vector dv = to_vector (lvalue->var->value); + + /* Convert source matrix 'sm' to source vector 'sv'. */ + if (!is_vector (sm)) + { + msg (SE, _("Can't assign matrix with dimensions (%zu,%zu) to subvector."), + sm->size1, sm->size2); + return false; + } + gsl_vector sv = to_vector (sm); + if (iv->n != sv.size) + { + msg (SE, _("Can't assign vector with %zu elements " + "to subvector with %zu."), sv.size, iv->n); + return false; + } + + /* Assign elements. */ + for (size_t x = 0; x < iv->n; x++) + gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x)); + return true; +} + +static bool +matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue, + struct index_vector *iv0, + struct index_vector *iv1, gsl_matrix *sm) +{ + gsl_matrix *dm = lvalue->var->value; + + /* Convert source matrix 'sm' to source vector 'sv'. */ + if (iv0->n != sm->size1) + { + msg (SE, _("Row index vector for assignment to %s has %zu elements " + "but source matrix has %zu rows."), + lvalue->var->name, iv0->n, sm->size1); + return false; + } + if (iv1->n != sm->size2) + { + msg (SE, _("Column index vector for assignment to %s has %zu elements " + "but source matrix has %zu columns."), + lvalue->var->name, iv1->n, sm->size2); + return false; + } + + /* Assign elements. */ + for (size_t y = 0; y < iv0->n; y++) + { + size_t f0 = iv0->indexes[y]; + for (size_t x = 0; x < iv1->n; x++) + { + size_t f1 = iv1->indexes[x]; + gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x)); + } + } + return true; +} + +/* Takes ownership of M. */ +static bool +matrix_lvalue_assign (struct matrix_lvalue *lvalue, + struct index_vector *iv0, struct index_vector *iv1, + gsl_matrix *sm) +{ + if (!lvalue->n_indexes) + { + gsl_matrix_free (lvalue->var->value); + lvalue->var->value = sm; + return true; + } + else + { + bool ok = (lvalue->n_indexes == 1 + ? matrix_lvalue_assign_vector (lvalue, iv0, sm) + : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm)); + free (iv0->indexes); + free (iv1->indexes); + gsl_matrix_free (sm); + return ok; + } +} + +static bool +matrix_lvalue_evaluate (struct matrix_lvalue *lvalue, + struct index_vector *iv0, + struct index_vector *iv1) +{ + *iv0 = INDEX_VECTOR_INIT; + *iv1 = INDEX_VECTOR_INIT; + if (!lvalue->n_indexes) + return true; + + /* Validate destination matrix exists and has the right shape. */ + gsl_matrix *dm = lvalue->var->value; + if (!dm) + { + msg (SE, _("Undefined variable %s."), lvalue->var->name); + return false; + } + else if (dm->size1 == 0 || dm->size2 == 0) + { + msg (SE, _("Cannot index matrix with dimensions (%zu,%zu)."), + dm->size1, dm->size2); + return false; + } + else if (lvalue->n_indexes == 1) + { + if (!is_vector (dm)) + { + msg (SE, _("Can't use vector indexing on matrix %s with " + "dimensions (%zu,%zu)."), + lvalue->var->name, dm->size1, dm->size2); + return false; + } + return matrix_lvalue_evaluate_vector (lvalue->indexes[0], + MAX (dm->size1, dm->size2), + IV_VECTOR, 0, iv0); + } + else + { + assert (lvalue->n_indexes == 2); + if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1, + IV_ROW, dm->size2, iv0)) + return false; + + if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2, + IV_COLUMN, dm->size1, iv1)) + { + free (iv0->indexes); + return false; + } + return true; + } +} + +/* Takes ownership of M. */ +static bool +matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm) +{ + struct index_vector iv0, iv1; + return (matrix_lvalue_evaluate (lvalue, &iv0, &iv1) + && matrix_lvalue_assign (lvalue, &iv0, &iv1, sm)); +} + + +/* Matrix command. */ + +struct matrix_cmds + { + struct matrix_cmd **commands; + size_t n; + }; + +struct matrix_cmd + { + enum matrix_cmd_type + { + MCMD_COMPUTE, + MCMD_PRINT, + MCMD_DO_IF, + MCMD_LOOP, + MCMD_BREAK, + MCMD_DISPLAY, + MCMD_RELEASE, + MCMD_SAVE, + MCMD_READ, + MCMD_WRITE, + MCMD_GET, + MCMD_MSAVE, + MCMD_MGET, + MCMD_EIGEN, + MCMD_SETDIAG, + MCMD_SVD, + } + type; + + union + { + struct compute_command + { + struct matrix_lvalue *lvalue; + struct matrix_expr *rvalue; + } + compute; + + struct print_command + { + struct matrix_expr *expression; + bool use_default_format; + struct fmt_spec format; + char *title; + int space; /* -1 means NEWPAGE. */ + + struct print_labels + { + struct string_array literals; /* CLABELS/RLABELS. */ + struct matrix_expr *expr; /* CNAMES/RNAMES. */ + } + *rlabels, *clabels; + } + print; + + struct do_if_command + { + struct do_if_clause + { + struct matrix_expr *condition; + struct matrix_cmds commands; + } + *clauses; + size_t n_clauses; + } + do_if; + + struct loop_command + { + /* Index. */ + struct matrix_var *var; + struct matrix_expr *start; + struct matrix_expr *end; + struct matrix_expr *increment; + + /* Loop conditions. */ + struct matrix_expr *top_condition; + struct matrix_expr *bottom_condition; + + /* Commands. */ + struct matrix_cmds commands; + } + loop; + + struct save_command + { + struct matrix_expr *expression; + struct file_handle *outfile; + struct string_array *variables; + struct matrix_expr *names; + struct stringi_set strings; + } + save; + + struct read_command + { + struct read_file *rf; + struct matrix_lvalue *dst; + struct matrix_expr *size; + int c1, c2; + enum fmt_type format; + int w; + int d; + bool symmetric; + bool reread; + } + read; + + struct write_command + { + struct matrix_expr *expression; + struct file_handle *outfile; + char *encoding; + int c1, c2; + enum fmt_type format; + int w; + int d; + bool triangular; + bool hold; + } + write; + + struct display_command + { + struct matrix_state *state; + bool dictionary; + bool status; + } + display; + + struct release_command + { + struct matrix_var **vars; + size_t n_vars; + } + release; + + struct get_command + { + struct matrix_lvalue *dst; + struct file_handle *file; + char *encoding; + struct string_array variables; + struct matrix_var *names; + + /* Treatment of missing values. */ + struct + { + enum + { + MGET_ERROR, /* Flag error on command. */ + MGET_ACCEPT, /* Accept without change, user-missing only. */ + MGET_OMIT, /* Drop this case. */ + MGET_RECODE /* Recode to 'substitute'. */ + } + treatment; + double substitute; /* MGET_RECODE only. */ + } + user, system; + } + get; + + struct msave_command + { + struct msave_common *common; + char *varname_; + struct matrix_expr *expr; + const char *rowtype; + struct matrix_expr *factors; + struct matrix_expr *splits; + } + msave; + + struct mget_command + { + struct matrix_state *state; + struct file_handle *file; + struct stringi_set rowtypes; + } + mget; + + struct eigen_command + { + struct matrix_expr *expr; + struct matrix_var *evec; + struct matrix_var *eval; + } + eigen; + + struct setdiag_command + { + struct matrix_var *dst; + struct matrix_expr *expr; + } + setdiag; + + struct svd_command + { + struct matrix_expr *expr; + struct matrix_var *u; + struct matrix_var *s; + struct matrix_var *v; + } + svd; + }; + }; + +static struct matrix_cmd *matrix_parse_command (struct matrix_state *); +static bool matrix_cmd_execute (struct matrix_cmd *); + +static void +matrix_cmd_destroy (struct matrix_cmd *cmd) +{ + /* XXX */ + free (cmd); +} + +static struct matrix_cmd * +matrix_parse_compute (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_COMPUTE, + .compute = { .lvalue = NULL }, + }; + + struct compute_command *compute = &cmd->compute; + compute->lvalue = matrix_lvalue_parse (s); + if (!compute->lvalue) + goto error; + + if (!lex_force_match (s->lexer, T_EQUALS)) + goto error; + + compute->rvalue = matrix_parse_expr (s); + if (!compute->rvalue) + goto error; + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_compute (struct compute_command *compute) +{ + gsl_matrix *value = matrix_expr_evaluate (compute->rvalue); + if (!value) + return; + + matrix_lvalue_evaluate_and_assign (compute->lvalue, value); +} + +static void +print_labels_destroy (struct print_labels *labels) +{ + if (labels) + { + string_array_destroy (&labels->literals); + matrix_expr_destroy (labels->expr); + free (labels); + } +} + +static void +parse_literal_print_labels (struct matrix_state *s, + struct print_labels **labelsp) +{ + lex_match (s->lexer, T_EQUALS); + print_labels_destroy (*labelsp); + *labelsp = xzalloc (sizeof **labelsp); + while (lex_token (s->lexer) != T_SLASH + && lex_token (s->lexer) != T_ENDCMD + && lex_token (s->lexer) != T_STOP) + { + struct string label = DS_EMPTY_INITIALIZER; + while (lex_token (s->lexer) != T_COMMA + && lex_token (s->lexer) != T_SLASH + && lex_token (s->lexer) != T_ENDCMD + && lex_token (s->lexer) != T_STOP) + { + if (!ds_is_empty (&label)) + ds_put_byte (&label, ' '); + + if (lex_token (s->lexer) == T_STRING) + ds_put_cstr (&label, lex_tokcstr (s->lexer)); + else + { + char *rep = lex_next_representation (s->lexer, 0, 0); + ds_put_cstr (&label, rep); + free (rep); + } + lex_get (s->lexer); + } + string_array_append_nocopy (&(*labelsp)->literals, + ds_steal_cstr (&label)); + + if (!lex_match (s->lexer, T_COMMA)) + break; + } +} + +static bool +parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp) +{ + lex_match (s->lexer, T_EQUALS); + struct matrix_expr *e = matrix_parse_exp (s); + if (!e) + return false; + + print_labels_destroy (*labelsp); + *labelsp = xzalloc (sizeof **labelsp); + (*labelsp)->expr = e; + return true; +} + +static struct matrix_cmd * +matrix_parse_print (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_PRINT, + .print = { + .use_default_format = true, + } + }; + + if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD) + { + size_t depth = 0; + for (size_t i = 0; ; i++) + { + enum token_type t = lex_next_token (s->lexer, i); + if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY) + depth++; + else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth) + depth--; + else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP) + { + if (i > 0) + cmd->print.title = lex_next_representation (s->lexer, 0, i - 1); + break; + } + } + + cmd->print.expression = matrix_parse_exp (s); + if (!cmd->print.expression) + goto error; + } + + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "FORMAT")) + { + lex_match (s->lexer, T_EQUALS); + if (!parse_format_specifier (s->lexer, &cmd->print.format)) + goto error; + cmd->print.use_default_format = false; + } + else if (lex_match_id (s->lexer, "TITLE")) + { + lex_match (s->lexer, T_EQUALS); + if (!lex_force_string (s->lexer)) + goto error; + free (cmd->print.title); + cmd->print.title = ss_xstrdup (lex_tokss (s->lexer)); + lex_get (s->lexer); + } + else if (lex_match_id (s->lexer, "SPACE")) + { + lex_match (s->lexer, T_EQUALS); + if (lex_match_id (s->lexer, "NEWPAGE")) + cmd->print.space = -1; + else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX)) + { + cmd->print.space = lex_integer (s->lexer); + lex_get (s->lexer); + } + else + goto error; + } + else if (lex_match_id (s->lexer, "RLABELS")) + parse_literal_print_labels (s, &cmd->print.rlabels); + else if (lex_match_id (s->lexer, "CLABELS")) + parse_literal_print_labels (s, &cmd->print.clabels); + else if (lex_match_id (s->lexer, "RNAMES")) + { + if (!parse_expr_print_labels (s, &cmd->print.rlabels)) + goto error; + } + else if (lex_match_id (s->lexer, "CNAMES")) + { + if (!parse_expr_print_labels (s, &cmd->print.clabels)) + goto error; + } + else + { + lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE", + "RLABELS", "CLABELS", "RNAMES", "CNAMES"); + goto error; + } + + } + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static bool +matrix_is_integer (const gsl_matrix *m) +{ + for (size_t y = 0; y < m->size1; y++) + for (size_t x = 0; x < m->size2; x++) + { + double d = gsl_matrix_get (m, y, x); + if (d != floor (d)) + return false; + } + return true; +} + +static double +matrix_max_magnitude (const gsl_matrix *m) +{ + double max = 0.0; + for (size_t y = 0; y < m->size1; y++) + for (size_t x = 0; x < m->size2; x++) + { + double d = fabs (gsl_matrix_get (m, y, x)); + if (d > max) + max = d; + } + return max; +} + +static bool +format_fits (struct fmt_spec format, double x) +{ + char *s = data_out (&(union value) { .f = x }, NULL, + &format, settings_get_fmt_settings ()); + bool fits = *s != '*' && !strchr (s, 'E'); + free (s); + return fits; +} + +static struct fmt_spec +default_format (const gsl_matrix *m, int *log_scale) +{ + *log_scale = 0; + + double max = matrix_max_magnitude (m); + + if (matrix_is_integer (m)) + for (int w = 1; w <= 12; w++) + { + struct fmt_spec format = { .type = FMT_F, .w = w }; + if (format_fits (format, -max)) + return format; + } + + if (max >= 1e9 || max <= 1e-4) + { + char s[64]; + snprintf (s, sizeof s, "%.1e", max); + + const char *e = strchr (s, 'e'); + if (e) + *log_scale = atoi (e + 1); + } + + return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 }; +} + +static char * +trimmed_string (double d) +{ + struct substring s = ss_buffer ((char *) &d, sizeof d); + ss_rtrim (&s, ss_cstr (" ")); + return ss_xstrdup (s); +} + +static struct string_array * +print_labels_get (const struct print_labels *labels, size_t n, + const char *prefix, bool truncate) +{ + if (!labels) + return NULL; + + struct string_array *out = xzalloc (sizeof *out); + if (labels->literals.n) + string_array_clone (out, &labels->literals); + else if (labels->expr) + { + gsl_matrix *m = matrix_expr_evaluate (labels->expr); + if (m && is_vector (m)) + { + gsl_vector v = to_vector (m); + for (size_t i = 0; i < v.size; i++) + string_array_append_nocopy (out, trimmed_string ( + gsl_vector_get (&v, i))); + } + gsl_matrix_free (m); + } + + while (out->n < n) + { + if (labels->expr) + string_array_append_nocopy (out, xasprintf ("%s%zu", prefix, + out->n + 1)); + else + string_array_append (out, ""); + } + while (out->n > n) + string_array_delete (out, out->n - 1); + + if (truncate) + for (size_t i = 0; i < out->n; i++) + { + char *s = out->strings[i]; + s[strnlen (s, 8)] = '\0'; + } + + return out; +} + +static void +matrix_cmd_print_space (int space) +{ + if (space < 0) + output_item_submit (page_break_item_create ()); + for (int i = 0; i < space; i++) + output_log ("%s", ""); +} + +static void +matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m, + struct fmt_spec format, int log_scale) +{ + matrix_cmd_print_space (print->space); + if (print->title) + output_log ("%s", print->title); + if (log_scale != 0) + output_log (" 10 ** %d X", log_scale); + + struct string_array *clabels = print_labels_get (print->clabels, + m->size2, "col", true); + if (clabels && format.w < 8) + format.w = 8; + struct string_array *rlabels = print_labels_get (print->rlabels, + m->size1, "row", true); + + if (clabels) + { + struct string line = DS_EMPTY_INITIALIZER; + if (rlabels) + ds_put_byte_multiple (&line, ' ', 8); + for (size_t x = 0; x < m->size2; x++) + ds_put_format (&line, " %*s", format.w, clabels->strings[x]); + output_log_nocopy (ds_steal_cstr (&line)); + } + + double scale = pow (10.0, log_scale); + bool numeric = fmt_is_numeric (format.type); + for (size_t y = 0; y < m->size1; y++) + { + struct string line = DS_EMPTY_INITIALIZER; + if (rlabels) + ds_put_format (&line, "%-8s", rlabels->strings[y]); + + for (size_t x = 0; x < m->size2; x++) + { + double f = gsl_matrix_get (m, y, x); + char *s = (numeric + ? data_out (&(union value) { .f = f / scale}, NULL, + &format, settings_get_fmt_settings ()) + : trimmed_string (f)); + ds_put_format (&line, " %s", s); + free (s); + } + output_log_nocopy (ds_steal_cstr (&line)); + } + + string_array_destroy (rlabels); + string_array_destroy (clabels); +} + +static void +create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis, + const struct print_labels *print_labels, size_t n, + const char *name, const char *prefix) +{ + struct string_array *labels = print_labels_get (print_labels, n, prefix, + false); + struct pivot_dimension *d = pivot_dimension_create (table, axis, name); + for (size_t i = 0; i < n; i++) + pivot_category_create_leaf ( + d->root, (labels + ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX) + : pivot_value_new_integer (i + 1))); + if (!labels) + d->hide_all_labels = true; + string_array_destroy (labels); +} + +static void +matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m, + struct fmt_spec format, int log_scale) +{ + struct pivot_table *table = pivot_table_create__ ( + pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX), + "Matrix Print"); + + create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1, + N_("Rows"), "row"); + create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2, + N_("Columns"), "col"); + + struct pivot_footnote *footnote = NULL; + if (log_scale != 0) + { + char *s = xasprintf ("× 10**%d", log_scale); + footnote = pivot_table_create_footnote ( + table, pivot_value_new_user_text_nocopy (s)); + } + + double scale = pow (10.0, log_scale); + bool numeric = fmt_is_numeric (format.type); + for (size_t y = 0; y < m->size1; y++) + for (size_t x = 0; x < m->size2; x++) + { + double f = gsl_matrix_get (m, y, x); + struct pivot_value *v; + if (numeric) + { + v = pivot_value_new_number (f / scale); + v->numeric.format = format; + } + else + v = pivot_value_new_user_text_nocopy (trimmed_string (f)); + if (footnote) + pivot_value_add_footnote (v, footnote); + pivot_table_put2 (table, y, x, v); + } + + pivot_table_submit (table); +} + +static void +matrix_cmd_execute_print (const struct print_command *print) +{ + if (print->expression) + { + gsl_matrix *m = matrix_expr_evaluate (print->expression); + if (!m) + return; + + int log_scale = 0; + struct fmt_spec format = (print->use_default_format + ? default_format (m, &log_scale) + : print->format); + + if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT) + matrix_cmd_print_text (print, m, format, log_scale); + else + matrix_cmd_print_tables (print, m, format, log_scale); + + gsl_matrix_free (m); + } + else + { + matrix_cmd_print_space (print->space); + if (print->title) + output_log ("%s", print->title); + } +} + +/* DO IF. */ + +static bool +matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c, + const char *command_name, + const char *stop1, const char *stop2) +{ + lex_end_of_command (s->lexer); + lex_discard_rest_of_command (s->lexer); + + size_t allocated = 0; + for (;;) + { + while (lex_token (s->lexer) == T_ENDCMD) + lex_get (s->lexer); + + if (lex_at_phrase (s->lexer, stop1) + || (stop2 && lex_at_phrase (s->lexer, stop2))) + return true; + + if (lex_at_phrase (s->lexer, "END MATRIX")) + { + msg (SE, _("Premature END MATRIX within %s."), command_name); + return false; + } + + struct matrix_cmd *cmd = matrix_parse_command (s); + if (!cmd) + return false; + + if (c->n >= allocated) + c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands); + c->commands[c->n++] = cmd; + } +} + +static bool +matrix_parse_do_if_clause (struct matrix_state *s, + struct do_if_command *ifc, + bool parse_condition, + size_t *allocated_clauses) +{ + if (ifc->n_clauses >= *allocated_clauses) + ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses, + sizeof *ifc->clauses); + struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++]; + *c = (struct do_if_clause) { .condition = NULL }; + + if (parse_condition) + { + c->condition = matrix_parse_expr (s); + if (!c->condition) + return false; + } + + return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF"); +} + +static struct matrix_cmd * +matrix_parse_do_if (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_DO_IF, + .do_if = { .n_clauses = 0 } + }; + + size_t allocated_clauses = 0; + do + { + if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses)) + goto error; + } + while (lex_match_phrase (s->lexer, "ELSE IF")); + + if (lex_match_id (s->lexer, "ELSE") + && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses)) + goto error; + + if (!lex_match_phrase (s->lexer, "END IF")) + NOT_REACHED (); + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static bool +matrix_cmd_execute_do_if (struct do_if_command *cmd) +{ + for (size_t i = 0; i < cmd->n_clauses; i++) + { + struct do_if_clause *c = &cmd->clauses[i]; + if (c->condition) + { + double d; + if (!matrix_expr_evaluate_scalar (c->condition, + i ? "ELSE IF" : "DO IF", + &d) || d <= 0) + continue; + } + + for (size_t j = 0; j < c->commands.n; j++) + if (!matrix_cmd_execute (c->commands.commands[j])) + return false; + return true; + } + return true; +} + +static struct matrix_cmd * +matrix_parse_loop (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } }; + + struct loop_command *loop = &cmd->loop; + if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS) + { + struct substring name = lex_tokss (s->lexer); + loop->var = matrix_var_lookup (s, name); + if (!loop->var) + loop->var = matrix_var_create (s, name); + + lex_get (s->lexer); + lex_get (s->lexer); + + loop->start = matrix_parse_expr (s); + if (!loop->start || !lex_force_match (s->lexer, T_TO)) + goto error; + + loop->end = matrix_parse_expr (s); + if (!loop->end) + goto error; + + if (lex_match (s->lexer, T_BY)) + { + loop->increment = matrix_parse_expr (s); + if (!loop->increment) + goto error; + } + } + + if (lex_match_id (s->lexer, "IF")) + { + loop->top_condition = matrix_parse_expr (s); + if (!loop->top_condition) + goto error; + } + + bool was_in_loop = s->in_loop; + s->in_loop = true; + bool ok = matrix_parse_commands (s, &loop->commands, "LOOP", + "END LOOP", NULL); + s->in_loop = was_in_loop; + if (!ok) + goto error; + + if (!lex_match_phrase (s->lexer, "END LOOP")) + NOT_REACHED (); + + if (lex_match_id (s->lexer, "IF")) + { + loop->bottom_condition = matrix_parse_expr (s); + if (!loop->bottom_condition) + goto error; + } + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_loop (struct loop_command *cmd) +{ + long int value = 0; + long int end = 0; + long int increment = 1; + if (cmd->var) + { + if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value) + || !matrix_expr_evaluate_integer (cmd->end, "TO", &end) + || (cmd->increment + && !matrix_expr_evaluate_integer (cmd->increment, "BY", + &increment))) + return; + + if (increment > 0 ? value > end + : increment < 0 ? value < end + : true) + return; + } + + int mxloops = settings_get_mxloops (); + for (int i = 0; i < mxloops; i++) + { + if (cmd->var) + { + if (cmd->var->value + && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1)) + { + gsl_matrix_free (cmd->var->value); + cmd->var->value = NULL; + } + if (!cmd->var->value) + cmd->var->value = gsl_matrix_alloc (1, 1); + + gsl_matrix_set (cmd->var->value, 0, 0, value); + } + + if (cmd->top_condition) + { + double d; + if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF", + &d) || d <= 0) + return; + } + + for (size_t j = 0; j < cmd->commands.n; j++) + if (!matrix_cmd_execute (cmd->commands.commands[j])) + return; + + if (cmd->bottom_condition) + { + double d; + if (!matrix_expr_evaluate_scalar (cmd->bottom_condition, + "END LOOP IF", &d) || d > 0) + return; + } + + if (cmd->var) + { + value += increment; + if (increment > 0 ? value > end : value < end) + return; + } + } +} + +static struct matrix_cmd * +matrix_parse_break (struct matrix_state *s) +{ + if (!s->in_loop) + { + msg (SE, _("BREAK not inside LOOP.")); + return NULL; + } + + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { .type = MCMD_BREAK }; + return cmd; +} + +static struct matrix_cmd * +matrix_parse_display (struct matrix_state *s) +{ + bool dictionary = false; + bool status = false; + while (lex_token (s->lexer) == T_ID) + { + if (lex_match_id (s->lexer, "DICTIONARY")) + dictionary = true; + else if (lex_match_id (s->lexer, "STATUS")) + status = true; + else + { + lex_error_expecting (s->lexer, "DICTIONARY", "STATUS"); + return NULL; + } + } + if (!dictionary && !status) + dictionary = status = true; + + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_DISPLAY, + .display = { + .state = s, + .dictionary = dictionary, + .status = status, + } + }; + return cmd; +} + +static int +compare_matrix_var_pointers (const void *a_, const void *b_) +{ + const struct matrix_var *const *ap = a_; + const struct matrix_var *const *bp = b_; + const struct matrix_var *a = *ap; + const struct matrix_var *b = *bp; + return strcmp (a->name, b->name); +} + +static void +matrix_cmd_execute_display (struct display_command *cmd) +{ + const struct matrix_state *s = cmd->state; + + struct pivot_table *table = pivot_table_create (N_("Matrix Variables")); + pivot_dimension_create ( + table, PIVOT_AXIS_COLUMN, N_("Property"), + N_("Rows"), N_("Columns"), N_("Size (kB)")); + + const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars); + size_t n_vars = 0; + + const struct matrix_var *var; + HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars) + if (var->value) + vars[n_vars++] = var; + qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers); + + struct pivot_dimension *rows = pivot_dimension_create ( + table, PIVOT_AXIS_ROW, N_("Variable")); + for (size_t i = 0; i < n_vars; i++) + { + const struct matrix_var *var = vars[i]; + pivot_category_create_leaf ( + rows->root, pivot_value_new_user_text (var->name, SIZE_MAX)); + + size_t r = var->value->size1; + size_t c = var->value->size2; + double values[] = { r, c, r * c * sizeof (double) / 1024 }; + for (size_t j = 0; j < sizeof values / sizeof *values; j++) + pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j])); + } + pivot_table_submit (table); +} + +static struct matrix_cmd * +matrix_parse_release (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE }; + + size_t allocated_vars = 0; + while (lex_token (s->lexer) == T_ID) + { + struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer)); + if (var) + { + if (cmd->release.n_vars >= allocated_vars) + cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars, + sizeof *cmd->release.vars); + cmd->release.vars[cmd->release.n_vars++] = var; + } + else + lex_error (s->lexer, _("Variable name expected.")); + lex_get (s->lexer); + + if (!lex_match (s->lexer, T_COMMA)) + break; + } + + return cmd; +} + +static void +matrix_cmd_execute_release (struct release_command *cmd) +{ + for (size_t i = 0; i < cmd->n_vars; i++) + { + struct matrix_var *var = cmd->vars[i]; + gsl_matrix_free (var->value); + var->value = NULL; + } +} + +static struct matrix_cmd * +matrix_parse_save (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_SAVE, + .save = { + .strings = STRINGI_SET_INITIALIZER (cmd->save.strings) + } + }; + + struct save_command *save = &cmd->save; + save->expression = matrix_parse_exp (s); + if (!save->expression) + goto error; + + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "OUTFILE")) + { + lex_match (s->lexer, T_EQUALS); + + fh_unref (save->outfile); + save->outfile = (lex_match (s->lexer, T_ASTERISK) + ? fh_inline_file () + : fh_parse (s->lexer, FH_REF_FILE, s->session)); + if (!save->outfile) + goto error; + } + else if (lex_match_id (s->lexer, "VARIABLES")) + { + lex_match (s->lexer, T_EQUALS); + + char **names; + size_t n; + struct dictionary *d = dict_create (get_default_encoding ()); + bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n, + PV_NO_SCRATCH | PV_NO_DUPLICATE); + dict_unref (d); + if (!ok) + goto error; + + string_array_destroy (save->variables); + if (!save->variables) + save->variables = xmalloc (sizeof *save->variables); + *save->variables = (struct string_array) { + .strings = names, + .n = n, + .allocated = n, + }; + } + else if (lex_match_id (s->lexer, "NAMES")) + { + lex_match (s->lexer, T_EQUALS); + matrix_expr_destroy (save->names); + save->names = matrix_parse_exp (s); + if (!save->names) + goto error; + } + else if (lex_match_id (s->lexer, "STRINGS")) + { + lex_match (s->lexer, T_EQUALS); + while (lex_token (s->lexer) == T_ID) + { + stringi_set_insert (&save->strings, lex_tokcstr (s->lexer)); + lex_get (s->lexer); + if (!lex_match (s->lexer, T_COMMA)) + break; + } + } + else + { + lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES", + "STRINGS"); + goto error; + } + } + + if (!save->outfile) + { + if (s->prev_save_outfile) + save->outfile = fh_ref (s->prev_save_outfile); + else + { + lex_sbc_missing ("OUTFILE"); + goto error; + } + } + fh_unref (s->prev_save_outfile); + s->prev_save_outfile = fh_ref (save->outfile); + + if (save->variables && save->names) + { + msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES.")); + matrix_expr_destroy (save->names); + save->names = NULL; + } + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_save (const struct save_command *save) +{ + assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */ + + gsl_matrix *m = matrix_expr_evaluate (save->expression); + if (!m) + return; + + bool ok = true; + struct dictionary *dict = dict_create (get_default_encoding ()); + struct string_array names = { .n = 0 }; + if (save->variables) + string_array_clone (&names, save->variables); + else if (save->names) + { + gsl_matrix *nm = matrix_expr_evaluate (save->names); + if (nm && is_vector (nm)) + { + gsl_vector nv = to_vector (nm); + for (size_t i = 0; i < nv.size; i++) + { + char *name = trimmed_string (gsl_vector_get (&nv, i)); + if (dict_id_is_valid (dict, name, true)) + string_array_append_nocopy (&names, name); + else + ok = false; + } + } + } + + struct stringi_set strings; + stringi_set_clone (&strings, &save->strings); + + for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++) + { + char tmp_name[64]; + const char *name; + if (i < names.n) + name = names.strings[i]; + else + { + snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1); + name = tmp_name; + } + + int width = stringi_set_delete (&strings, name) ? 8 : 0; + struct variable *var = dict_create_var (dict, name, width); + if (!var) + { + msg (SE, _("Duplicate variable name %s in SAVE statement."), + name); + ok = false; + } + } + + if (!stringi_set_is_empty (&strings)) + { + const char *example = stringi_set_node_get_string ( + stringi_set_first (&strings)); + msg (SE, ngettext ("STRINGS specified a variable %s, but no variable " + "with that name was found on SAVE.", + "STRINGS specified %2$zu variables, including %1$s, " + "whose names were not found on SAVE.", + stringi_set_count (&strings)), + example, stringi_set_count (&strings)); + ok = false; + } + stringi_set_destroy (&strings); + string_array_destroy (&names); + + if (!ok) + { + gsl_matrix_free (m); + dict_unref (dict); + return; + } + + struct casewriter *writer = any_writer_open (save->outfile, dict); + if (!writer) + { + gsl_matrix_free (m); + dict_unref (dict); + return; + } + + const struct caseproto *proto = dict_get_proto (dict); + for (size_t y = 0; y < m->size1; y++) + { + struct ccase *c = case_create (proto); + for (size_t x = 0; x < m->size2; x++) + { + double d = gsl_matrix_get (m, y, x); + int width = caseproto_get_width (proto, x); + union value *value = case_data_rw_idx (c, x); + if (width == 0) + value->f = d; + else + memcpy (value->s, &d, width); + } + casewriter_write (writer, c); + } + casewriter_destroy (writer); + + gsl_matrix_free (m); + dict_unref (dict); +} + +static struct matrix_cmd * +matrix_parse_read (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_READ, + .read = { .format = FMT_F }, + }; + + struct file_handle *fh = NULL; + char *encoding = NULL; + struct read_command *read = &cmd->read; + read->dst = matrix_lvalue_parse (s); + if (!read->dst) + goto error; + + int by = 0; + int repetitions = 0; + int record_width = 0; + bool seen_format = false; + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "FILE")) + { + lex_match (s->lexer, T_EQUALS); + + fh_unref (fh); + fh = fh_parse (s->lexer, FH_REF_FILE, NULL); + if (!fh) + goto error; + } + else if (lex_match_id (s->lexer, "ENCODING")) + { + lex_match (s->lexer, T_EQUALS); + if (!lex_force_string (s->lexer)) + goto error; + + free (encoding); + encoding = ss_xstrdup (lex_tokss (s->lexer)); + + lex_get (s->lexer); + } + else if (lex_match_id (s->lexer, "FIELD")) + { + lex_match (s->lexer, T_EQUALS); + + if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX)) + goto error; + read->c1 = lex_integer (s->lexer); + lex_get (s->lexer); + if (!lex_force_match (s->lexer, T_TO) + || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX)) + goto error; + read->c2 = lex_integer (s->lexer) + 1; + lex_get (s->lexer); + + record_width = read->c2 - read->c1; + if (lex_match (s->lexer, T_BY)) + { + if (!lex_force_int_range (s->lexer, "BY", 1, + read->c2 - read->c1)) + goto error; + by = lex_integer (s->lexer); + lex_get (s->lexer); + + if (record_width % by) + { + msg (SE, _("BY %d does not evenly divide record width %d."), + by, record_width); + goto error; + } + } + else + by = 0; + } + else if (lex_match_id (s->lexer, "SIZE")) + { + lex_match (s->lexer, T_EQUALS); + read->size = matrix_parse_exp (s); + if (!read->size) + goto error; + } + else if (lex_match_id (s->lexer, "MODE")) + { + lex_match (s->lexer, T_EQUALS); + if (lex_match_id (s->lexer, "RECTANGULAR")) + read->symmetric = false; + else if (lex_match_id (s->lexer, "SYMMETRIC")) + read->symmetric = true; + else + { + lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC"); + goto error; + } + } + else if (lex_match_id (s->lexer, "REREAD")) + read->reread = true; + else if (lex_match_id (s->lexer, "FORMAT")) + { + if (seen_format) + { + lex_sbc_only_once ("FORMAT"); + goto error; + } + seen_format = true; + + lex_match (s->lexer, T_EQUALS); + + if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer)) + goto error; + + const char *p = lex_tokcstr (s->lexer); + if (c_isdigit (p[0])) + { + repetitions = atoi (p); + p += strspn (p, "0123456789"); + if (!fmt_from_name (p, &read->format)) + { + lex_error (s->lexer, _("Unknown format %s."), p); + goto error; + } + lex_get (s->lexer); + } + else if (!fmt_from_name (p, &read->format)) + { + struct fmt_spec format; + if (!parse_format_specifier (s->lexer, &format)) + goto error; + read->format = format.type; + read->w = format.w; + read->d = format.d; + } + } + else + { + lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE", + "REREAD", "FORMAT"); + goto error; + } + } + + if (!read->c1) + { + lex_sbc_missing ("FIELD"); + goto error; + } + + if (!read->dst->n_indexes && !read->size) + { + msg (SE, _("SIZE is required for reading data into a full matrix " + "(as opposed to a submatrix).")); + goto error; + } + + if (!fh) + { + if (s->prev_read_file) + fh = fh_ref (s->prev_read_file); + else + { + lex_sbc_missing ("FILE"); + goto error; + } + } + fh_unref (s->prev_read_file); + s->prev_read_file = fh_ref (fh); + + read->rf = read_file_create (s, fh); + if (encoding) + { + free (read->rf->encoding); + read->rf->encoding = encoding; + encoding = NULL; + } + + /* Field width may be specified in multiple ways: + + 1. BY on FIELD. + 2. The format on FORMAT. + 3. The repetition factor on FORMAT. + + (2) and (3) are mutually exclusive. + + If more than one of these is present, they must agree. If none of them is + present, then free-field format is used. + */ + if (repetitions > record_width) + { + msg (SE, _("%d repetitions cannot fit in record width %d."), + repetitions, record_width); + goto error; + } + int w = (repetitions ? record_width / repetitions + : read->w ? read->w + : by); + if (by && w != by) + { + if (repetitions) + msg (SE, _("FORMAT specifies %d repetitions with record width %d, " + "which implies field width %d, " + "but BY specifies field width %d."), + repetitions, record_width, w, by); + else + msg (SE, _("FORMAT specifies field width %d but BY specifies %d."), + w, by); + goto error; + } + read->w = w; + return cmd; + +error: + matrix_cmd_destroy (cmd); + free (encoding); + return NULL; +} + +static void +matrix_read_set_field (struct read_command *read, struct dfm_reader *reader, + gsl_matrix *m, struct substring p, size_t y, size_t x) +{ + const char *input_encoding = dfm_reader_get_encoding (reader); + union value v; + char *error = data_in (p, input_encoding, read->format, + settings_get_fmt_settings (), &v, 0, NULL); + /* XXX report error if value is missing */ + if (error) + msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error); + else + { + gsl_matrix_set (m, y, x, v.f); + if (read->symmetric && x != y) + gsl_matrix_set (m, x, y, v.f); + } +} + +static bool +matrix_read_line (struct read_command *read, struct dfm_reader *reader, + struct substring *line) +{ + if (dfm_eof (reader)) + { + msg (SE, _("Unexpected end of file reading matrix data.")); + return false; + } + dfm_expand_tabs (reader); + *line = ss_substr (dfm_get_record (reader), + read->c1 - 1, read->c2 - read->c1); + return true; +} + +static void +matrix_read (struct read_command *read, struct dfm_reader *reader, + gsl_matrix *m) +{ + for (size_t y = 0; y < m->size1; y++) + { + size_t nx = read->symmetric ? y + 1 : m->size2; + + struct substring line = ss_empty (); + for (size_t x = 0; x < nx; x++) + { + struct substring p; + if (!read->w) + { + for (;;) + { + ss_ltrim (&line, ss_cstr (" ,")); + if (!ss_is_empty (line)) + break; + if (!matrix_read_line (read, reader, &line)) + return; + dfm_forward_record (reader); + } + + ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p); + } + else + { + if (!matrix_read_line (read, reader, &line)) + return; + size_t fields_per_line = (read->c2 - read->c1) / read->w; + int f = x % fields_per_line; + if (f == fields_per_line - 1) + dfm_forward_record (reader); + + p = ss_substr (line, read->w * f, read->w); + } + + matrix_read_set_field (read, reader, m, p, y, x); + } + + if (read->w) + dfm_forward_record (reader); + else + { + ss_ltrim (&line, ss_cstr (" ,")); + if (!ss_is_empty (line)) + msg (SW, _("Trailing garbage on line \"%.*s\""), + (int) line.length, line.string); + } + } +} + +static void +matrix_cmd_execute_read (struct read_command *read) +{ + struct index_vector iv0, iv1; + if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1)) + return; + + size_t size[2] = { SIZE_MAX, SIZE_MAX }; + if (read->size) + { + gsl_matrix *m = matrix_expr_evaluate (read->size); + if (!m) + return; + + if (!is_vector (m)) + { + msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector")); + gsl_matrix_free (m); + free (iv0.indexes); + free (iv1.indexes); + return; + } + + gsl_vector v = to_vector (m); + double d[2]; + if (v.size == 1) + { + d[0] = gsl_vector_get (&v, 0); + d[1] = 1; + } + else if (v.size == 2) + { + d[0] = gsl_vector_get (&v, 0); + d[1] = gsl_vector_get (&v, 1); + } + else + { + msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector")); + gsl_matrix_free (m); + free (iv0.indexes); + free (iv1.indexes); + return; + } + gsl_matrix_free (m); + + if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX) + { + msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]); + free (iv0.indexes); + free (iv1.indexes); + return; + } + + size[0] = d[0]; + size[1] = d[1]; + } + + if (read->dst->n_indexes) + { + size_t submatrix_size[2]; + if (read->dst->n_indexes == 2) + { + submatrix_size[0] = iv0.n; + submatrix_size[1] = iv1.n; + } + else if (read->dst->var->value->size1 == 1) + { + submatrix_size[0] = 1; + submatrix_size[1] = iv0.n; + } + else + { + submatrix_size[0] = iv0.n; + submatrix_size[1] = 1; + } + + if (read->size) + { + if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1]) + { + msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions " + "(%zu,%zu)."), + size[0], size[1], + submatrix_size[0], submatrix_size[1]); + free (iv0.indexes); + free (iv1.indexes); + return; + } + } + else + { + size[0] = submatrix_size[0]; + size[1] = submatrix_size[1]; + } + } + + struct dfm_reader *reader = read_file_open (read->rf); + if (read->reread) + dfm_reread_record (reader, 1); + + if (read->symmetric && size[0] != size[1]) + { + msg (SE, _("Cannot read matrix with non-square dimensions (%zu,%zu) " + "using READ with MODE=SYMMETRIC."), + size[0], size[1]); + free (iv0.indexes); + free (iv1.indexes); + return; + } + gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]); + matrix_read (read, reader, tmp); + matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp); +} + +static struct matrix_cmd * +matrix_parse_write (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_WRITE, + .write = { .format = FMT_F }, + }; + + struct write_command *write = &cmd->write; + write->expression = matrix_parse_exp (s); + if (!write->expression) + goto error; + + int by = 0; + int repetitions = 0; + int record_width = 0; + bool seen_format = false; + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "OUTFILE")) + { + lex_match (s->lexer, T_EQUALS); + + fh_unref (write->outfile); + write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL); + if (!write->outfile) + goto error; + } + else if (lex_match_id (s->lexer, "ENCODING")) + { + lex_match (s->lexer, T_EQUALS); + if (!lex_force_string (s->lexer)) + goto error; + + free (write->encoding); + write->encoding = ss_xstrdup (lex_tokss (s->lexer)); + + lex_get (s->lexer); + } + else if (lex_match_id (s->lexer, "FIELD")) + { + lex_match (s->lexer, T_EQUALS); + + if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX)) + goto error; + write->c1 = lex_integer (s->lexer); + lex_get (s->lexer); + if (!lex_force_match (s->lexer, T_TO) + || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX)) + goto error; + write->c2 = lex_integer (s->lexer) + 1; + lex_get (s->lexer); + + record_width = write->c2 - write->c1; + if (lex_match (s->lexer, T_BY)) + { + if (!lex_force_int_range (s->lexer, "BY", 1, + write->c2 - write->c1)) + goto error; + by = lex_integer (s->lexer); + lex_get (s->lexer); + + if (record_width % by) + { + msg (SE, _("BY %d does not evenly divide record width %d."), + by, record_width); + goto error; + } + } + else + by = 0; + } + else if (lex_match_id (s->lexer, "MODE")) + { + lex_match (s->lexer, T_EQUALS); + if (lex_match_id (s->lexer, "RECTANGULAR")) + write->triangular = false; + else if (lex_match_id (s->lexer, "TRIANGULAR")) + write->triangular = true; + else + { + lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR"); + goto error; + } + } + else if (lex_match_id (s->lexer, "HOLD")) + write->hold = true; + else if (lex_match_id (s->lexer, "FORMAT")) + { + if (seen_format) + { + lex_sbc_only_once ("FORMAT"); + goto error; + } + seen_format = true; + + lex_match (s->lexer, T_EQUALS); + + if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer)) + goto error; + + const char *p = lex_tokcstr (s->lexer); + if (c_isdigit (p[0])) + { + repetitions = atoi (p); + p += strspn (p, "0123456789"); + if (!fmt_from_name (p, &write->format)) + { + lex_error (s->lexer, _("Unknown format %s."), p); + goto error; + } + lex_get (s->lexer); + } + else if (!fmt_from_name (p, &write->format)) + { + struct fmt_spec format; + if (!parse_format_specifier (s->lexer, &format)) + goto error; + write->format = format.type; + write->w = format.w; + write->d = format.d; + } + } + else + { + lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE", + "HOLD", "FORMAT"); + goto error; + } + } + + if (!write->c1) + { + lex_sbc_missing ("FIELD"); + goto error; + } + + if (!write->outfile) + { + if (s->prev_write_outfile) + write->outfile = fh_ref (s->prev_write_outfile); + else + { + lex_sbc_missing ("OUTFILE"); + goto error; + } + } + fh_unref (s->prev_write_outfile); + s->prev_write_outfile = fh_ref (write->outfile); + + /* Field width may be specified in multiple ways: + + 1. BY on FIELD. + 2. The format on FORMAT. + 3. The repetition factor on FORMAT. + + (2) and (3) are mutually exclusive. + + If more than one of these is present, they must agree. If none of them is + present, then free-field format is used. + */ + if (repetitions > record_width) + { + msg (SE, _("%d repetitions cannot fit in record width %d."), + repetitions, record_width); + goto error; + } + int w = (repetitions ? record_width / repetitions + : write->w ? write->w + : by); + if (by && w != by) + { + if (repetitions) + msg (SE, _("FORMAT specifies %d repetitions with record width %d, " + "which implies field width %d, " + "but BY specifies field width %d."), + repetitions, record_width, w, by); + else + msg (SE, _("FORMAT specifies field width %d but BY specifies %d."), + w, by); + goto error; + } + write->w = w; + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_write (struct write_command *write) +{ + gsl_matrix *m = matrix_expr_evaluate (write->expression); + if (!m) + return; + + if (write->triangular && m->size1 != m->size2) + { + msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but " + "the matrix to be written has dimensions (%zu,%zu)."), + m->size1, m->size2); + gsl_matrix_free (m); + return; + } + + struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding); + if (!writer) + { + gsl_matrix_free (m); + return; + } + + const struct fmt_settings *settings = settings_get_fmt_settings (); + struct fmt_spec format = { + .type = write->format, + .w = write->w ? write->w : 40, + .d = write->d + }; + struct u8_line line = U8_LINE_INITIALIZER; + for (size_t y = 0; y < m->size1; y++) + { + size_t nx = write->triangular ? y + 1 : m->size2; + int x0 = write->c1; + for (size_t x = 0; x < nx; x++) + { + /* XXX string values */ + union value v = { .f = gsl_matrix_get (m, y, x) }; + char *s = (write->w + ? data_out (&v, NULL, &format, settings) + : data_out_stretchy (&v, NULL, &format, settings, NULL)); + size_t len = strlen (s); + int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8); + if (width + x0 > write->c2) + { + dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length); + u8_line_clear (&line); + x0 = write->c1; + } + u8_line_put (&line, x0, x0 + width, s, len); + free (s); + + x0 += write->w ? write->w : width + 1; + } + + dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length); + u8_line_clear (&line); + } + u8_line_destroy (&line); + dfm_close_writer (writer); + + gsl_matrix_free (m); +} + +static struct matrix_cmd * +matrix_parse_get (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_GET, + .get = { + .user = { .treatment = MGET_ERROR }, + .system = { .treatment = MGET_ERROR }, + } + }; + + struct get_command *get = &cmd->get; + get->dst = matrix_lvalue_parse (s); + if (!get->dst) + goto error; + + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "FILE")) + { + if (get->variables.n) + { + lex_error (s->lexer, _("FILE must precede VARIABLES")); + goto error; + } + lex_match (s->lexer, T_EQUALS); + + fh_unref (get->file); + get->file = fh_parse (s->lexer, FH_REF_FILE, s->session); + if (!get->file) + goto error; + } + else if (lex_match_id (s->lexer, "ENCODING")) + { + if (get->variables.n) + { + lex_error (s->lexer, _("ENCODING must precede VARIABLES")); + goto error; + } + lex_match (s->lexer, T_EQUALS); + if (!lex_force_string (s->lexer)) + goto error; + + free (get->encoding); + get->encoding = ss_xstrdup (lex_tokss (s->lexer)); + + lex_get (s->lexer); + } + else if (lex_match_id (s->lexer, "VARIABLES")) + { + lex_match (s->lexer, T_EQUALS); + + struct dictionary *dict = NULL; + if (!get->file) + { + dict = dataset_dict (s->dataset); + if (dict_get_var_cnt (dict) == 0) + { + lex_error (s->lexer, _("GET cannot read empty active file.")); + goto error; + } + } + else + { + struct casereader *reader = any_reader_open_and_decode ( + get->file, get->encoding, &dict, NULL); + if (!reader) + goto error; + casereader_destroy (reader); + } + + struct variable **vars; + size_t n_vars; + bool ok = parse_variables (s->lexer, dict, &vars, &n_vars, + PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH); + if (!ok) + { + dict_unref (dict); + goto error; + } + + string_array_clear (&get->variables); + for (size_t i = 0; i < n_vars; i++) + string_array_append (&get->variables, var_get_name (vars[i])); + free (vars); + dict_unref (dict); + } + else if (lex_match_id (s->lexer, "NAMES")) + { + lex_match (s->lexer, T_EQUALS); + if (!lex_force_id (s->lexer)) + goto error; + + struct substring name = lex_tokss (s->lexer); + get->names = matrix_var_lookup (s, name); + if (!get->names) + get->names = matrix_var_create (s, name); + lex_get (s->lexer); + } + else if (lex_match_id (s->lexer, "MISSING")) + { + lex_match (s->lexer, T_EQUALS); + if (lex_match_id (s->lexer, "ACCEPT")) + get->user.treatment = MGET_ACCEPT; + else if (lex_match_id (s->lexer, "OMIT")) + get->user.treatment = MGET_OMIT; + else if (lex_is_number (s->lexer)) + { + get->user.treatment = MGET_RECODE; + get->user.substitute = lex_number (s->lexer); + lex_get (s->lexer); + } + else + { + lex_error (s->lexer, NULL); + goto error; + } + } + else if (lex_match_id (s->lexer, "SYSMIS")) + { + lex_match (s->lexer, T_EQUALS); + if (lex_match_id (s->lexer, "OMIT")) + get->user.treatment = MGET_OMIT; + else if (lex_is_number (s->lexer)) + { + get->user.treatment = MGET_RECODE; + get->user.substitute = lex_number (s->lexer); + lex_get (s->lexer); + } + else + { + lex_error (s->lexer, NULL); + goto error; + } + } + else + { + lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES", + "MISSING", "SYSMIS"); + goto error; + } + } + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_get (struct get_command *get) +{ + assert (get->file); /* XXX */ + + struct dictionary *dict; + struct casereader *reader = any_reader_open_and_decode ( + get->file, get->encoding, &dict, NULL); + if (!reader) + return; + + const struct variable **vars = xnmalloc ( + get->variables.n ? get->variables.n : dict_get_var_cnt (dict), + sizeof *vars); + size_t n_vars = 0; + + if (get->variables.n) + { + for (size_t i = 0; i < get->variables.n; i++) + { + const char *name = get->variables.strings[i]; + const struct variable *var = dict_lookup_var (dict, name); + if (!var) + { + msg (SE, _("GET: Data file does not contain variable %s."), + name); + dict_unref (dict); + free (vars); + return; + } + if (!var_is_numeric (var)) + { + msg (SE, _("GET: Variable %s is not numeric."), name); + dict_unref (dict); + free (vars); + return; + } + vars[n_vars++] = var; + } + } + else + { + for (size_t i = 0; i < dict_get_var_cnt (dict); i++) + { + const struct variable *var = dict_get_var (dict, i); + if (!var_is_numeric (var)) + { + msg (SE, _("GET: Variable %s is not numeric."), + var_get_name (var)); + dict_unref (dict); + free (vars); + return; + } + vars[n_vars++] = var; + } + } + + size_t n_rows = 0; + gsl_matrix *m = gsl_matrix_alloc (4, n_vars); + long long int casenum = 1; + bool error = false; + for (struct ccase *c = casereader_read (reader); c; + c = casereader_read (reader), casenum++) + { + if (n_rows >= m->size1) + { + gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars); + for (size_t y = 0; y < n_rows; y++) + for (size_t x = 0; x < n_vars; x++) + gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x)); + gsl_matrix_free (m); + m = p; + } + + bool keep = true; + for (size_t x = 0; x < n_vars; x++) + { + const struct variable *var = vars[x]; + double d = case_num (c, var); + if (d == SYSMIS) + { + if (get->system.treatment == MGET_RECODE) + d = get->system.substitute; + else if (get->system.treatment == MGET_OMIT) + keep = false; + else + { + msg (SE, _("GET: Variable %s in case %lld " + "is system-missing."), + var_get_name (var), casenum); + error = true; + } + } + else if (var_is_num_missing (var, d, MV_USER)) + { + if (get->user.treatment == MGET_RECODE) + d = get->user.substitute; + else if (get->user.treatment == MGET_OMIT) + keep = false; + else if (get->user.treatment != MGET_ACCEPT) + { + msg (SE, _("GET: Variable %s in case %lld has user-missing " + "value %g."), + var_get_name (var), casenum, d); + error = true; + } + } + gsl_matrix_set (m, n_rows, x, d); + } + case_unref (c); + if (error) + break; + if (keep) + n_rows++; + } + casereader_destroy (reader); + if (!error) + { + m->size1 = n_rows; + matrix_lvalue_evaluate_and_assign (get->dst, m); + } + else + gsl_matrix_free (m); + dict_unref (dict); + free (vars); +} + +static const char * +match_rowtype (struct lexer *lexer) +{ + static const char *rowtypes[] = { + "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT" + }; + size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes; + + for (size_t i = 0; i < n_rowtypes; i++) + if (lex_match_id (lexer, rowtypes[i])) + return rowtypes[i]; + + lex_error_expecting_array (lexer, rowtypes, n_rowtypes); + return NULL; +} + +static bool +parse_var_names (struct lexer *lexer, struct string_array *sa) +{ + lex_match (lexer, T_EQUALS); + + string_array_clear (sa); + + struct dictionary *dict = dict_create (get_default_encoding ()); + dict_create_var_assert (dict, "ROWTYPE_", 8); + dict_create_var_assert (dict, "VARNAME_", 8); + char **names; + size_t n_names; + bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names, + PV_NO_DUPLICATE | PV_NO_SCRATCH); + dict_unref (dict); + + if (ok) + { + string_array_clear (sa); + sa->strings = names; + sa->n = sa->allocated = n_names; + } + return ok; +} + +static void +msave_common_uninit (struct msave_common *common) +{ + if (common) + { + fh_unref (common->outfile); + string_array_destroy (&common->variables); + string_array_destroy (&common->fnames); + string_array_destroy (&common->snames); + } +} + +static bool +compare_variables (const char *keyword, + const struct string_array *new, + const struct string_array *old) +{ + if (new->n) + { + if (!old->n) + { + msg (SE, _("%s may only be specified on MSAVE if it was specified " + "on the first MSAVE within MATRIX."), keyword); + return false; + } + else if (!string_array_equal_case (old, new)) + { + msg (SE, _("%s must specify the same variables each time within " + "a given MATRIX."), keyword); + return false; + } + } + return true; +} + +static struct matrix_cmd * +matrix_parse_msave (struct matrix_state *s) +{ + struct msave_common common = { .outfile = NULL }; + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } }; + + struct msave_command *msave = &cmd->msave; + if (lex_token (s->lexer) == T_ID) + msave->varname_ = ss_xstrdup (lex_tokss (s->lexer)); + msave->expr = matrix_parse_exp (s); + if (!msave->expr) + return NULL; + + while (lex_match (s->lexer, T_SLASH)) + { + if (lex_match_id (s->lexer, "TYPE")) + { + lex_match (s->lexer, T_EQUALS); + + msave->rowtype = match_rowtype (s->lexer); + if (!msave->rowtype) + goto error; + } + else if (lex_match_id (s->lexer, "OUTFILE")) + { + lex_match (s->lexer, T_EQUALS); + + fh_unref (common.outfile); + common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL); + if (!common.outfile) + goto error; + } + else if (lex_match_id (s->lexer, "VARIABLES")) + { + if (!parse_var_names (s->lexer, &common.variables)) + goto error; + } + else if (lex_match_id (s->lexer, "FNAMES")) + { + if (!parse_var_names (s->lexer, &common.fnames)) + goto error; + } + else if (lex_match_id (s->lexer, "SNAMES")) + { + if (!parse_var_names (s->lexer, &common.snames)) + goto error; + } + else if (lex_match_id (s->lexer, "SPLIT")) + { + lex_match (s->lexer, T_EQUALS); + + matrix_expr_destroy (msave->splits); + msave->splits = matrix_parse_exp (s); + if (!msave->splits) + goto error; + } + else if (lex_match_id (s->lexer, "FACTOR")) + { + lex_match (s->lexer, T_EQUALS); + + matrix_expr_destroy (msave->factors); + msave->factors = matrix_parse_exp (s); + if (!msave->factors) + goto error; + } + else + { + lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES", + "FNAMES", "SNAMES", "SPLIT", "FACTOR"); + goto error; + } + } + if (!msave->rowtype) + { + lex_sbc_missing ("TYPE"); + goto error; + } + common.has_splits = msave->splits || common.snames.n; + common.has_factors = msave->factors || common.fnames.n; + + struct msave_common *c = s->common ? s->common : &common; + if (c->fnames.n && !msave->factors) + { + msg (SE, _("FNAMES requires FACTOR.")); + goto error; + } + if (c->snames.n && !msave->splits) + { + msg (SE, _("SNAMES requires SPLIT.")); + goto error; + } + if (c->has_factors && !common.has_factors) + { + msg (SE, _("%s is required because it was present on the first " + "MSAVE in this MATRIX command."), "FACTOR"); + goto error; + } + if (c->has_splits && !common.has_splits) + { + msg (SE, _("%s is required because it was present on the first " + "MSAVE in this MATRIX command."), "SPLIT"); + goto error; + } + + if (!s->common) + { + if (!common.outfile) + { + lex_sbc_missing ("OUTFILE"); + goto error; + } + s->common = xmemdup (&common, sizeof common); + } + else + { + if (common.outfile) + { + bool same = common.outfile == s->common->outfile; + fh_unref (common.outfile); + if (!same) + { + msg (SE, _("OUTFILE must name the same file on each MSAVE " + "within a single MATRIX command.")); + goto error; + } + } + if (!compare_variables ("VARIABLES", + &common.variables, &s->common->variables) + || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames) + || !compare_variables ("SNAMES", &common.snames, &s->common->snames)) + goto error; + msave_common_uninit (&common); + } + msave->common = s->common; + if (!msave->varname_) + msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames); + return cmd; + +error: + msave_common_uninit (&common); + matrix_cmd_destroy (cmd); + return NULL; +} + +static gsl_vector * +matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name) +{ + gsl_matrix *m = matrix_expr_evaluate (e); + if (!m) + return NULL; + + if (!is_vector (m)) + { + msg (SE, _("%s expression must evaluate to vector, not a matrix with " + "dimensions (%zu,%zu)."), + name, m->size1, m->size2); + gsl_matrix_free (m); + return NULL; + } + + return matrix_to_vector (m); +} + +static const char * +msave_add_vars (struct dictionary *d, const struct string_array *vars) +{ + for (size_t i = 0; i < vars->n; i++) + if (!dict_create_var (d, vars->strings[i], 0)) + return vars->strings[i]; + return NULL; +} + +static struct dictionary * +msave_create_dict (const struct msave_common *common) +{ + struct dictionary *dict = dict_create (get_default_encoding ()); + + const char *dup_split = msave_add_vars (dict, &common->snames); + if (dup_split) + { + msg (SE, _("Duplicate SPLIT variable name %s."), dup_split); + goto error; + } + + dict_create_var_assert (dict, "ROWTYPE_", 8); + + const char *dup_factor = msave_add_vars (dict, &common->fnames); + if (dup_factor) + { + msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor); + goto error; + } + + dict_create_var_assert (dict, "VARNAME_", 8); + + const char *dup_var = msave_add_vars (dict, &common->variables); + if (dup_var) + { + msg (SE, _("Duplicate or invalid variable name %s."), dup_var); + goto error; + } + + return dict; + +error: + dict_unref (dict); + return NULL; +} + +static void +matrix_cmd_execute_msave (struct msave_command *msave) +{ + struct msave_common *common = msave->common; + gsl_matrix *m = NULL; + gsl_vector *factors = NULL; + gsl_vector *splits = NULL; + + m = matrix_expr_evaluate (msave->expr); + if (!m) + goto error; + + if (!common->variables.n) + for (size_t i = 0; i < m->size2; i++) + string_array_append_nocopy (&common->variables, + xasprintf ("COL%zu", i + 1)); + + if (m->size2 != common->variables.n) + { + msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."), + m->size2, common->variables.n); + goto error; + } + + if (msave->factors) + { + factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR"); + if (!factors) + goto error; + + if (!common->fnames.n) + for (size_t i = 0; i < factors->size; i++) + string_array_append_nocopy (&common->fnames, + xasprintf ("FAC%zu", i + 1)); + + if (factors->size != common->fnames.n) + { + msg (SE, _("There are %zu factor variables, " + "but %zu split values were supplied."), + common->fnames.n, factors->size); + goto error; + } + } + + if (msave->splits) + { + splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT"); + if (!splits) + goto error; + + if (!common->fnames.n) + for (size_t i = 0; i < splits->size; i++) + string_array_append_nocopy (&common->fnames, + xasprintf ("SPL%zu", i + 1)); + + if (splits->size != common->fnames.n) + { + msg (SE, _("There are %zu split variables, " + "but %zu split values were supplied."), + common->fnames.n, splits->size); + goto error; + } + } + + if (!common->writer) + { + struct dictionary *dict = msave_create_dict (common); + if (!dict) + goto error; + + common->writer = any_writer_open (common->outfile, dict); + if (!common->writer) + { + dict_unref (dict); + goto error; + } + + common->dict = dict; + } + + for (size_t y = 0; y < m->size1; y++) + { + struct ccase *c = case_create (dict_get_proto (common->dict)); + size_t idx = 0; + + /* Split variables */ + if (splits) + for (size_t i = 0; i < splits->size; i++) + case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i); + + /* ROWTYPE_. */ + buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8, + msave->rowtype, ' '); + + /* Factors. */ + if (factors) + for (size_t i = 0; i < factors->size; i++) + *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i); + + /* VARNAME_. */ + buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8, + msave->varname_, ' '); + + /* Continuous variables. */ + for (size_t x = 0; x < m->size2; x++) + case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x); + casewriter_write (common->writer, c); + } + +error: + gsl_matrix_free (m); + gsl_vector_free (factors); + gsl_vector_free (splits); +} + +static struct matrix_cmd * +matrix_parse_mget (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } }; + + struct mget_command *mget = &cmd->mget; + + for (;;) + { + if (lex_match_id (s->lexer, "FILE")) + { + lex_match (s->lexer, T_EQUALS); + + fh_unref (mget->file); + mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session); + if (!mget->file) + goto error; + } + else if (lex_match_id (s->lexer, "TYPE")) + { + lex_match (s->lexer, T_EQUALS); + while (lex_token (s->lexer) != T_SLASH + && lex_token (s->lexer) != T_ENDCMD) + { + const char *rowtype = match_rowtype (s->lexer); + if (!rowtype) + goto error; + + stringi_set_insert (&mget->rowtypes, rowtype); + } + } + else + { + lex_error_expecting (s->lexer, "FILE", "TYPE"); + goto error; + } + if (lex_token (s->lexer) == T_ENDCMD) + break; + + if (!lex_force_match (s->lexer, T_SLASH)) + goto error; + } + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static const struct variable * +get_a8_var (const struct dictionary *d, const char *name) +{ + const struct variable *v = dict_lookup_var (d, name); + if (!v) + { + msg (SE, _("Matrix data file lacks %s variable."), name); + return NULL; + } + if (var_get_width (v) != 8) + { + msg (SE, _("%s variable in matrix data file must be " + "8-byte string, but it has width %d."), + name, var_get_width (v)); + return NULL; + } + return v; +} + +static bool +is_rowtype (const union value *v, const char *rowtype) +{ + struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8); + ss_rtrim (&vs, ss_cstr (" ")); + return ss_equals_case (vs, ss_cstr (rowtype)); +} + +static void +matrix_mget_commit_var (struct ccase **rows, size_t n_rows, + const struct dictionary *d, + const struct variable *rowtype_var, + struct matrix_state *s, size_t si, size_t fi, + size_t cs, size_t cn) +{ + if (!n_rows) + return; + + const union value *rowtype_ = case_data (rows[0], rowtype_var); + const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV" + : is_rowtype (rowtype_, "CORR") ? "CR" + : is_rowtype (rowtype_, "MEAN") ? "MN" + : is_rowtype (rowtype_, "STDDEV") ? "SD" + : is_rowtype (rowtype_, "N") ? "NC" + : "CN"); + + struct string name = DS_EMPTY_INITIALIZER; + ds_put_cstr (&name, name_prefix); + if (fi > 0) + ds_put_format (&name, "F%zu", fi); + if (si > 0) + ds_put_format (&name, "S%zu", si); + + struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name)); + if (!mv) + mv = matrix_var_create (s, ds_ss (&name)); + else if (mv->value) + { + msg (SW, _("Matrix data file contains variable with existing name %s."), + ds_cstr (&name)); + goto exit; + } + + gsl_matrix *m = gsl_matrix_alloc (n_rows, cn); + size_t n_missing = 0; + for (size_t y = 0; y < n_rows; y++) + { + for (size_t x = 0; x < cn; x++) + { + struct variable *var = dict_get_var (d, cs + x); + double value = case_num (rows[y], var); + if (var_is_num_missing (var, value, MV_ANY)) + { + n_missing++; + value = 0.0; + } + gsl_matrix_set (m, y, x, value); + } + } + + if (n_missing) + msg (SE, ngettext ("Matrix data file variable %s contains a missing " + "value, which was treated as zero.", + "Matrix data file variable %s contains %zu missing " + "values, which were treated as zero.", n_missing), + ds_cstr (&name), n_missing); + mv->value = m; + +exit: + ds_destroy (&name); + for (size_t y = 0; y < n_rows; y++) + case_unref (rows[y]); +} + +static bool +var_changed (const struct ccase *ca, const struct ccase *cb, + const struct variable *var) +{ + return !value_equal (case_data (ca, var), case_data (cb, var), + var_get_width (var)); +} + +static bool +vars_changed (const struct ccase *ca, const struct ccase *cb, + const struct dictionary *d, + size_t first_var, size_t n_vars) +{ + for (size_t i = 0; i < n_vars; i++) + { + const struct variable *v = dict_get_var (d, first_var + i); + if (var_changed (ca, cb, v)) + return true; + } + return false; +} + +static void +matrix_cmd_execute_mget (struct mget_command *mget) +{ + struct dictionary *d; + struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8", + &d, NULL); + if (!r) + return; + + const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_"); + const struct variable *varname_ = get_a8_var (d, "VARNAME_"); + if (!rowtype_ || !varname_) + goto exit; + + if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_)) + { + msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file.")); + goto exit; + } + if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d)) + { + msg (SE, _("Matrix data file contains no continuous variables.")); + goto exit; + } + + for (size_t i = 0; i < dict_get_var_cnt (d); i++) + { + const struct variable *v = dict_get_var (d, i); + if (v != rowtype_ && v != varname_ && var_get_width (v) != 0) + { + msg (SE, + _("Matrix data file contains unexpected string variable %s."), + var_get_name (v)); + goto exit; + } + } + + /* SPLIT variables. */ + size_t ss = 0; + size_t sn = var_get_dict_index (rowtype_); + struct ccase *sc = NULL; + size_t si = 0; + + /* FACTOR variables. */ + size_t fs = var_get_dict_index (rowtype_) + 1; + size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1; + struct ccase *fc = NULL; + size_t fi = 0; + + /* Continuous variables. */ + size_t cs = var_get_dict_index (varname_) + 1; + size_t cn = dict_get_var_cnt (d) - cs; + struct ccase *cc = NULL; + + /* Matrix. */ + struct ccase **rows = NULL; + size_t allocated_rows = 0; + size_t n_rows = 0; + + struct ccase *c; + while ((c = casereader_read (r)) != NULL) + { + bool sd = vars_changed (sc, c, d, ss, sn); + bool fd = sd || vars_changed (fc, c, d, fs, fn); + bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_); + if (sd) + { + si++; + case_unref (sc); + sc = case_ref (c); + } + if (fd) + { + fi++; + case_unref (fc); + fc = case_ref (c); + } + if (md) + { + matrix_mget_commit_var (rows, n_rows, d, rowtype_, + mget->state, si, fi, cs, cn); + n_rows = 0; + case_unref (cc); + cc = case_ref (c); + } + + if (n_rows >= allocated_rows) + rows = x2nrealloc (rows, &allocated_rows, sizeof *rows); + rows[n_rows++] = c; + } + matrix_mget_commit_var (rows, n_rows, d, rowtype_, + mget->state, si, fi, cs, cn); + free (rows); + +exit: + casereader_destroy (r); +} + +static bool +matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp) +{ + if (!lex_force_id (s->lexer)) + return false; + + *varp = matrix_var_lookup (s, lex_tokss (s->lexer)); + if (!*varp) + *varp = matrix_var_create (s, lex_tokss (s->lexer)); + lex_get (s->lexer); + return true; +} + +static struct matrix_cmd * +matrix_parse_eigen (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_EIGEN, + .eigen = { .expr = NULL } + }; + + struct eigen_command *eigen = &cmd->eigen; + if (!lex_force_match (s->lexer, T_LPAREN)) + goto error; + eigen->expr = matrix_parse_expr (s); + if (!eigen->expr + || !lex_force_match (s->lexer, T_COMMA) + || !matrix_parse_dst_var (s, &eigen->evec) + || !lex_force_match (s->lexer, T_COMMA) + || !matrix_parse_dst_var (s, &eigen->eval) + || !lex_force_match (s->lexer, T_RPAREN)) + goto error; + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_eigen (struct eigen_command *eigen) +{ + gsl_matrix *A = matrix_expr_evaluate (eigen->expr); + if (!A) + return; + if (!is_symmetric (A)) + { + msg (SE, _("Argument of EIGEN must be symmetric.")); + gsl_matrix_free (A); + return; + } + + gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1); + gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1); + gsl_vector v_eval = to_vector (eval); + gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2); + gsl_eigen_symmv (A, &v_eval, evec, w); + gsl_eigen_symmv_free (w); + + gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC); + + gsl_matrix_free (A); + + gsl_matrix_free (eigen->eval->value); + eigen->eval->value = eval; + + gsl_matrix_free (eigen->evec->value); + eigen->evec->value = evec; +} + +static struct matrix_cmd * +matrix_parse_setdiag (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_SETDIAG, + .setdiag = { .dst = NULL } + }; + + struct setdiag_command *setdiag = &cmd->setdiag; + if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer)) + goto error; + + setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer)); + if (!setdiag->dst) + { + lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer)); + goto error; + } + lex_get (s->lexer); + + if (!lex_force_match (s->lexer, T_COMMA)) + goto error; + + setdiag->expr = matrix_parse_expr (s); + if (!setdiag->expr) + goto error; + + if (!lex_force_match (s->lexer, T_RPAREN)) + goto error; + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_setdiag (struct setdiag_command *setdiag) +{ + gsl_matrix *dst = setdiag->dst->value; + if (!dst) + { + msg (SE, _("SETDIAG destination matrix %s is uninitialized."), + setdiag->dst->name); + return; + } + + gsl_matrix *src = matrix_expr_evaluate (setdiag->expr); + if (!src) + return; + + size_t n = MIN (dst->size1, dst->size2); + if (is_scalar (src)) + { + double d = to_scalar (src); + for (size_t i = 0; i < n; i++) + gsl_matrix_set (dst, i, i, d); + } + else if (is_vector (src)) + { + gsl_vector v = to_vector (src); + for (size_t i = 0; i < n && i < v.size; i++) + gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i)); + } + else + msg (SE, _("SETDIAG argument 2 must be a scalar or a vector but it " + "has dimensions (%zu,%zu)."), + src->size1, src->size2); + gsl_matrix_free (src); +} + +static struct matrix_cmd * +matrix_parse_svd (struct matrix_state *s) +{ + struct matrix_cmd *cmd = xmalloc (sizeof *cmd); + *cmd = (struct matrix_cmd) { + .type = MCMD_SVD, + .svd = { .expr = NULL } + }; + + struct svd_command *svd = &cmd->svd; + if (!lex_force_match (s->lexer, T_LPAREN)) + goto error; + svd->expr = matrix_parse_expr (s); + if (!svd->expr + || !lex_force_match (s->lexer, T_COMMA) + || !matrix_parse_dst_var (s, &svd->u) + || !lex_force_match (s->lexer, T_COMMA) + || !matrix_parse_dst_var (s, &svd->s) + || !lex_force_match (s->lexer, T_COMMA) + || !matrix_parse_dst_var (s, &svd->v) + || !lex_force_match (s->lexer, T_RPAREN)) + goto error; + + return cmd; + +error: + matrix_cmd_destroy (cmd); + return NULL; +} + +static void +matrix_cmd_execute_svd (struct svd_command *svd) +{ + gsl_matrix *m = matrix_expr_evaluate (svd->expr); + if (!m) + return; + + if (m->size1 >= m->size2) + { + gsl_matrix *A = m; + gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2); + gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2); + gsl_vector Sv = gsl_matrix_diagonal (S).vector; + gsl_vector *work = gsl_vector_alloc (A->size2); + gsl_linalg_SV_decomp (A, V, &Sv, work); + gsl_vector_free (work); + + matrix_var_set (svd->u, A); + matrix_var_set (svd->s, S); + matrix_var_set (svd->v, V); + } + else + { + gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1); + gsl_matrix_transpose_memcpy (At, m); + gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2); + gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2); + gsl_vector Stv = gsl_matrix_diagonal (St).vector; + gsl_vector *work = gsl_vector_alloc (At->size2); + gsl_linalg_SV_decomp (At, Vt, &Stv, work); + gsl_vector_free (work); + + matrix_var_set (svd->v, At); + matrix_var_set (svd->s, St); + matrix_var_set (svd->u, Vt); + } +} + +static bool +matrix_cmd_execute (struct matrix_cmd *cmd) +{ + switch (cmd->type) + { + case MCMD_COMPUTE: + matrix_cmd_execute_compute (&cmd->compute); + break; + + case MCMD_PRINT: + matrix_cmd_execute_print (&cmd->print); + break; + + case MCMD_DO_IF: + return matrix_cmd_execute_do_if (&cmd->do_if); + + case MCMD_LOOP: + matrix_cmd_execute_loop (&cmd->loop); + break; + + case MCMD_BREAK: + return false; + + case MCMD_DISPLAY: + matrix_cmd_execute_display (&cmd->display); + break; + + case MCMD_RELEASE: + matrix_cmd_execute_release (&cmd->release); + break; + + case MCMD_SAVE: + matrix_cmd_execute_save (&cmd->save); + break; + + case MCMD_READ: + matrix_cmd_execute_read (&cmd->read); + break; + + case MCMD_WRITE: + matrix_cmd_execute_write (&cmd->write); + break; + + case MCMD_GET: + matrix_cmd_execute_get (&cmd->get); + break; + + case MCMD_MSAVE: + matrix_cmd_execute_msave (&cmd->msave); + break; + + case MCMD_MGET: + matrix_cmd_execute_mget (&cmd->mget); + break; + + case MCMD_EIGEN: + matrix_cmd_execute_eigen (&cmd->eigen); + break; + + case MCMD_SETDIAG: + matrix_cmd_execute_setdiag (&cmd->setdiag); + break; + + case MCMD_SVD: + matrix_cmd_execute_svd (&cmd->svd); + break; + } + + return true; +} + +struct matrix_command_name + { + const char *name; + struct matrix_cmd *(*parse) (struct matrix_state *); + }; + +static const struct matrix_command_name * +matrix_parse_command_name (struct lexer *lexer) +{ + static const struct matrix_command_name commands[] = { + { "COMPUTE", matrix_parse_compute }, + { "CALL EIGEN", matrix_parse_eigen }, + { "CALL SETDIAG", matrix_parse_setdiag }, + { "CALL SVD", matrix_parse_svd }, + { "PRINT", matrix_parse_print }, + { "DO IF", matrix_parse_do_if }, + { "LOOP", matrix_parse_loop }, + { "BREAK", matrix_parse_break }, + { "READ", matrix_parse_read }, + { "WRITE", matrix_parse_write }, + { "GET", matrix_parse_get }, + { "SAVE", matrix_parse_save }, + { "MGET", matrix_parse_mget }, + { "MSAVE", matrix_parse_msave }, + { "DISPLAY", matrix_parse_display }, + { "RELEASE", matrix_parse_release }, + }; + static size_t n = sizeof commands / sizeof *commands; + + for (const struct matrix_command_name *c = commands; c < &commands[n]; c++) + if (lex_match_phrase (lexer, c->name)) + return c; + return NULL; +} + +static struct matrix_cmd * +matrix_parse_command (struct matrix_state *s) +{ + size_t nesting_level = SIZE_MAX; + + struct matrix_cmd *c = NULL; + const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer); + if (!cmd) + lex_error (s->lexer, _("Unknown matrix command.")); + else if (!cmd->parse) + lex_error (s->lexer, _("Matrix command %s is not yet implemented."), + cmd->name); + else + { + nesting_level = output_open_group ( + group_item_create_nocopy (utf8_to_title (cmd->name), + utf8_to_title (cmd->name))); + c = cmd->parse (s); + } + + if (c) + lex_end_of_command (s->lexer); + lex_discard_rest_of_command (s->lexer); + if (nesting_level != SIZE_MAX) + output_close_groups (nesting_level); + + return c; +} + +int +cmd_matrix (struct lexer *lexer, struct dataset *ds) +{ + if (!lex_force_match (lexer, T_ENDCMD)) + return CMD_FAILURE; + + struct matrix_state state = { + .session = dataset_session (ds), + .lexer = lexer, + .vars = HMAP_INITIALIZER (state.vars), + }; + + for (;;) + { + while (lex_match (lexer, T_ENDCMD)) + continue; + if (lex_token (lexer) == T_STOP) + { + msg (SE, _("Unexpected end of input expecting matrix command.")); + break; + } + + if (lex_match_phrase (lexer, "END MATRIX")) + break; + + struct matrix_cmd *c = matrix_parse_command (&state); + if (c) + { + matrix_cmd_execute (c); + matrix_cmd_destroy (c); + } + } + + if (state.common) + { + dict_unref (state.common->dict); + casewriter_destroy (state.common->writer); + free (state.common); + } + + return CMD_SUCCESS; +} diff --git a/src/language/utilities/host.c b/src/language/utilities/host.c index 477f9796a0..d730e02634 100644 --- a/src/language/utilities/host.c +++ b/src/language/utilities/host.c @@ -261,8 +261,7 @@ run_command (const char *command, struct timespec timeout) if (end > output && end[-1] == '\n') end[-1] = '\0'; - output_log ("%s", output); - free (output); + output_log_nocopy (output); } free (locale_output); diff --git a/src/language/utilities/set.c b/src/language/utilities/set.c index a29e5c2ef0..9b18bfb0df 100644 --- a/src/language/utilities/set.c +++ b/src/language/utilities/set.c @@ -669,6 +669,24 @@ show_LOCALE (const struct dataset *ds UNUSED) return xstrdup (get_default_encoding ()); } +static bool +parse_MDISPLAY (struct lexer *lexer) +{ + int mdisplay = force_parse_enum (lexer, + "TEXT", SETTINGS_MDISPLAY_TEXT, + "TABLES", SETTINGS_MDISPLAY_TABLES); + if (mdisplay >= 0) + settings_set_mdisplay (mdisplay); + return mdisplay >= 0; +} + +static char * +show_MDISPLAY (const struct dataset *ds UNUSED) +{ + return xstrdup (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT + ? "TEXT" : "TABLES"); +} + static bool parse_MESSAGES (struct lexer *lexer) { @@ -1177,6 +1195,7 @@ static const struct setting settings[] = { { "JOURNAL", parse_JOURNAL, show_JOURNAL }, { "LENGTH", parse_LENGTH, show_LENGTH }, { "LOCALE", parse_LOCALE, show_LOCALE }, + { "MDISPLAY", parse_MDISPLAY, show_MDISPLAY }, { "MESSAGES", parse_MESSAGES, show_MESSAGES }, { "MEXPAND", parse_MEXPAND, show_MEXPAND }, { "MITERATE", parse_MITERATE, show_MITERATE }, diff --git a/src/libpspp/i18n.c b/src/libpspp/i18n.c index 69162f14f0..ac195cc517 100644 --- a/src/libpspp/i18n.c +++ b/src/libpspp/i18n.c @@ -829,7 +829,15 @@ utf8_hash_case_bytes (const char *s, size_t n, unsigned int basis) unsigned int utf8_hash_case_string (const char *s, unsigned int basis) { - return utf8_hash_case_bytes (s, strlen (s), basis); + return utf8_hash_case_substring (ss_cstr (s), basis); +} + +/* Returns a hash value for UTF-8 string S, with lowercase and uppercase + letters treated as equal, starting from BASIS. */ +unsigned int +utf8_hash_case_substring (struct substring s, unsigned int basis) +{ + return utf8_hash_case_bytes (s.string, s.length, basis); } /* Compares UTF-8 strings A and B case-insensitively. @@ -837,7 +845,13 @@ utf8_hash_case_string (const char *s, unsigned int basis) int utf8_strcasecmp (const char *a, const char *b) { - return utf8_strncasecmp (a, strlen (a), b, strlen (b)); + return utf8_sscasecmp (ss_cstr (a), ss_cstr (b)); +} + +int +utf8_sscasecmp (struct substring a, struct substring b) +{ + return utf8_strncasecmp (a.string, a.length, b.string, b.length); } /* Compares UTF-8 strings A (with length AN) and B (with length BN) diff --git a/src/libpspp/i18n.h b/src/libpspp/i18n.h index 232c5dc166..ee8be2fd63 100644 --- a/src/libpspp/i18n.h +++ b/src/libpspp/i18n.h @@ -18,6 +18,7 @@ #define I18N_H #include "libpspp/compiler.h" +#include "libpspp/str.h" #include #include @@ -74,7 +75,10 @@ const char *uc_name (ucs4_t uc, char buffer[16]); unsigned int utf8_hash_case_bytes (const char *, size_t n, unsigned int basis) WARN_UNUSED_RESULT; unsigned int utf8_hash_case_string (const char *, unsigned int basis) WARN_UNUSED_RESULT; +unsigned int utf8_hash_case_substring (struct substring, unsigned int basis) + WARN_UNUSED_RESULT; int utf8_strcasecmp (const char *, const char *); +int utf8_sscasecmp (struct substring, struct substring); int utf8_strncasecmp (const char *, size_t, const char *, size_t); int utf8_strverscasecmp (const char *, const char *); char *utf8_to_upper (const char *); diff --git a/src/libpspp/string-array.c b/src/libpspp/string-array.c index d4badf4638..d162593b2c 100644 --- a/src/libpspp/string-array.c +++ b/src/libpspp/string-array.c @@ -23,6 +23,7 @@ #include #include "libpspp/array.h" +#include "libpspp/i18n.h" #include "libpspp/str.h" #include "gl/xalloc.h" @@ -253,6 +254,36 @@ string_array_uniq (struct string_array *sa) sa->n = n; } +/* Returns true if A and B contain the same strings in the same order, + false otherwise. */ +bool +string_array_equal (const struct string_array *a, + const struct string_array *b) +{ + if (a->n != b->n) + return false; + + for (size_t i = 0; i < a->n; i++) + if (strcmp (a->strings[i], b->strings[i])) + return false; + return true; +} + +/* Returns true if A and B contain the same strings in the same order, + false otherwise. */ +bool +string_array_equal_case (const struct string_array *a, + const struct string_array *b) +{ + if (a->n != b->n) + return false; + + for (size_t i = 0; i < a->n; i++) + if (utf8_strcasecmp (a->strings[i], b->strings[i])) + return false; + return true; +} + /* Divides STRING into tokens at DELIMITERS and adds each token to SA. */ void string_array_parse (struct string_array *sa, struct substring string, diff --git a/src/libpspp/string-array.h b/src/libpspp/string-array.h index 353ddf5018..b4a6989d6d 100644 --- a/src/libpspp/string-array.h +++ b/src/libpspp/string-array.h @@ -63,6 +63,11 @@ void string_array_shrink (struct string_array *); void string_array_sort (struct string_array *); void string_array_uniq (struct string_array *); +bool string_array_equal (const struct string_array *, + const struct string_array *); +bool string_array_equal_case (const struct string_array *, + const struct string_array *); + void string_array_parse (struct string_array *, struct substring string, struct substring delimiters); char *string_array_join (const struct string_array *, const char *separator); diff --git a/src/libpspp/u8-line.h b/src/libpspp/u8-line.h index 2674e2e275..1da5f8869a 100644 --- a/src/libpspp/u8-line.h +++ b/src/libpspp/u8-line.h @@ -32,6 +32,8 @@ struct u8_line size_t width; /* Display width, in character positions. */ }; +#define U8_LINE_INITIALIZER { .s = DS_EMPTY_INITIALIZER }; + void u8_line_init (struct u8_line *); void u8_line_destroy (struct u8_line *); void u8_line_clear (struct u8_line *); diff --git a/src/output/driver.c b/src/output/driver.c index f40876449b..42842bfc20 100644 --- a/src/output/driver.c +++ b/src/output/driver.c @@ -366,6 +366,12 @@ output_log (const char *format, ...) char *s = xvasprintf (format, args); va_end (args); + output_log_nocopy (s); +} + +void +output_log_nocopy (char *s) +{ output_submit (text_item_create_nocopy (TEXT_ITEM_LOG, s, NULL)); } diff --git a/src/output/driver.h b/src/output/driver.h index 3e5cf193a4..c16774dd99 100644 --- a/src/output/driver.h +++ b/src/output/driver.h @@ -35,6 +35,7 @@ void output_submit (struct output_item *); void output_flush (void); void output_log (const char *, ...) PRINTF_FORMAT (1, 2); +void output_log_nocopy (char *); const char *output_get_title (void); void output_set_title (const char *); diff --git a/tests/automake.mk b/tests/automake.mk index 2ddea5c8b1..3f76cbc32a 100644 --- a/tests/automake.mk +++ b/tests/automake.mk @@ -399,6 +399,7 @@ TESTSUITE_AT = \ tests/language/stats/frequencies.at \ tests/language/stats/glm.at \ tests/language/stats/logistic.at \ + tests/language/stats/matrix.at \ tests/language/stats/means.at \ tests/language/stats/npar.at \ tests/language/stats/oneway.at \ diff --git a/tests/language/lexer/scan.at b/tests/language/lexer/scan.at index 90dea5d346..c877628fdf 100644 --- a/tests/language/lexer/scan.at +++ b/tests/language/lexer/scan.at @@ -170,13 +170,13 @@ LBRACK RBRACK EXP MACRO_PUNCT "%" -MACRO_PUNCT ":" -MACRO_PUNCT ";" +COLON +SEMICOLON MACRO_PUNCT "?" MACRO_PUNCT "_" MACRO_PUNCT "`" -MACRO_PUNCT "{" -MACRO_PUNCT "}" +LCURLY +RCURLY NOT STOP ]) diff --git a/tests/language/stats/matrix.at b/tests/language/stats/matrix.at new file mode 100644 index 0000000000..b86393bd20 --- /dev/null +++ b/tests/language/stats/matrix.at @@ -0,0 +1,2573 @@ +AT_BANNER([MATRIX]) + +AT_SETUP([MATRIX - empty matrices]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE a={}. +PRINT a. +COMPUTE b={a; 1; 2; 3}. +PRINT b. +COMPUTE c={a, 1, 2, 3}. +PRINT c. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +a + +b + 1 + 2 + 3 + +c + 1 2 3 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - submatrices as rvalues - all columns or all rows]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(1, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({1}, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 2}, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 2, 3}, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({1; 3; 2}, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 3, 3}, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(1:2, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(1:3, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}({}, :). + +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1}). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 2}). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 2, 3}). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1; 3; 2}). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 3, 3}). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1:2). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1:3). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {}). + +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, :). + +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(0, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 0). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(4, :). +PRINT {1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 4). + +PRINT {}(:,{}). +PRINT {}({},:). +PRINT {}({},{}). + +PRINT {1, 2, 3, 4}({1, 2; 3, 4}, :). +PRINT {1, 2, 3, 4}(:, {1, 2; 3, 4}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +{1, 2, 3; 4, 5, 6; 7, 8, 9}(1, :) + 1 2 3 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({1}, :) + 1 2 3 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 2}, :) + 1 2 3 + 4 5 6 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 2, 3}, :) + 1 2 3 + 4 5 6 + 7 8 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({1; 3; 2}, :) + 1 2 3 + 7 8 9 + 4 5 6 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({1, 3, 3}, :) + 1 2 3 + 7 8 9 + 7 8 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(1:2, :) + 1 2 3 + 4 5 6 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(1:3, :) + 1 2 3 + 4 5 6 + 7 8 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}({}, :) + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1) + 1 + 4 + 7 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1}) + 1 + 4 + 7 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 2}) + 1 2 + 4 5 + 7 8 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 2, 3}) + 1 2 3 + 4 5 6 + 7 8 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1; 3; 2}) + 1 3 2 + 4 6 5 + 7 9 8 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {1, 3, 3}) + 1 3 3 + 4 6 6 + 7 9 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1:2) + 1 2 + 4 5 + 7 8 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, 1:3) + 1 2 3 + 4 5 6 + 7 8 9 + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, {}) + + + +{1, 2, 3; 4, 5, 6; 7, 8, 9}(:, :) + 1 2 3 + 4 5 6 + 7 8 9 + +matrix.sps:24: error: MATRIX: 0 is not a valid row index for a matrix with +dimensions (3,3). + +matrix.sps:25: error: MATRIX: 0 is not a valid column index for a matrix with +dimensions (3,3). + +matrix.sps:26: error: MATRIX: 4 is not a valid row index for a matrix with +dimensions (3,3). + +matrix.sps:27: error: MATRIX: 4 is not a valid column index for a matrix with +dimensions (3,3). + +{}(:,{}) + +{}({},:) + +{}({},{}) + +matrix.sps:33: error: MATRIX: Matrix row index must be scalar or vector, not a +matrix with dimensions (2,2). + +matrix.sps:34: error: MATRIX: Matrix column index must be scalar or vector, not +a matrix with dimensions (2,2). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - COMPUTE submatrices as lvalues]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE y={1, 2, 3; 4, 5, 6; 7, 8, 9}. + +COMPUTE x1=y. +COMPUTE x1(1, :) = {11, 12, 13}. +PRINT x1. + +COMPUTE x2=y. +COMPUTE x2(2, :) = {14, 15, 16}. +PRINT x2. + +COMPUTE x3=y. +COMPUTE x3(3, :) = {17, 18, 19}. +PRINT x3. + +COMPUTE x4=y. +COMPUTE x4(:, 1) = {11; 14; 17}. +PRINT x4. + +COMPUTE x5=y. +COMPUTE x5(:, 2) = {12; 15; 18}. +PRINT x5. + +COMPUTE x6=y. +COMPUTE x6(:, 3) = {13; 16; 19}. +PRINT x6. + +COMPUTE x7=y. +COMPUTE x7(1, 1) = 11. +PRINT x7. + +COMPUTE x8=y. +COMPUTE x8(1:2, 2:3) = {12, 13; 15, 16}. +PRINT x8. + +COMPUTE x9=y. +COMPUTE x9({3, 1}, {2; 3}) = {18, 19; 12, 13}. +PRINT x9. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +x1 + 11 12 13 + 4 5 6 + 7 8 9 + +x2 + 1 2 3 + 14 15 16 + 7 8 9 + +x3 + 1 2 3 + 4 5 6 + 17 18 19 + +x4 + 11 2 3 + 14 5 6 + 17 8 9 + +x5 + 1 12 3 + 4 15 6 + 7 18 9 + +x6 + 1 2 13 + 4 5 16 + 7 8 19 + +x7 + 11 2 3 + 4 5 6 + 7 8 9 + +x8 + 1 12 13 + 4 15 16 + 7 8 9 + +x9 + 1 12 13 + 4 5 6 + 7 18 19 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - COMPUTE submatrices as lvalues - negative]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE x={1, 2, 3; 4, 5, 6; 7, 8, 9}. +COMPUTE x(1, :) = {}. +COMPUTE x(1, :) = 15. +COMPUTE x(1, :) = {11, 12}. +COMPUTE x(1, :) = {11, 12, 13, 14}. +COMPUTE x(:, 1) = {}. +COMPUTE x(:, 1) = 15. +COMPUTE x(:, 1) = {11, 12}. +COMPUTE x(:, 1) = {11, 12, 13, 14}. +COMPUTE x(:) = 1. +COMPUTE x(0, 1) = 1. +COMPUTE x(1, 0) = 1. +COMPUTE x({1, 0, 2}, 1) = {1; 2; 3}. +COMPUTE x(4, 3) = 1. +COMPUTE x(3, 4) = 1. +COMPUTE x({1, 2; 3, 4}, 5) = 1. +COMPUTE x(3, {1, 2; 3, 4}) = 1. +PRINT x. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +matrix.sps:3: error: MATRIX: Row index vector for assignment to x has 1 +elements but source matrix has 0 rows. + +matrix.sps:4: error: MATRIX: Column index vector for assignment to x has 3 +elements but source matrix has 1 columns. + +matrix.sps:5: error: MATRIX: Column index vector for assignment to x has 3 +elements but source matrix has 2 columns. + +matrix.sps:6: error: MATRIX: Column index vector for assignment to x has 3 +elements but source matrix has 4 columns. + +matrix.sps:7: error: MATRIX: Row index vector for assignment to x has 3 +elements but source matrix has 0 rows. + +matrix.sps:8: error: MATRIX: Row index vector for assignment to x has 3 +elements but source matrix has 1 rows. + +matrix.sps:9: error: MATRIX: Row index vector for assignment to x has 3 +elements but source matrix has 1 rows. + +matrix.sps:10: error: MATRIX: Row index vector for assignment to x has 3 +elements but source matrix has 1 rows. + +matrix.sps:11: error: MATRIX: Can't use vector indexing on matrix x with +dimensions (3,3). + +matrix.sps:12: error: MATRIX: 0 is not a valid row index for a matrix with +dimensions (3,3). + +matrix.sps:13: error: MATRIX: 0 is not a valid column index for a matrix with +dimensions (3,3). + +matrix.sps:14: error: MATRIX: 0 is not a valid row index for a matrix with +dimensions (3,3). + +matrix.sps:15: error: MATRIX: 4 is not a valid row index for a matrix with +dimensions (3,3). + +matrix.sps:16: error: MATRIX: 4 is not a valid column index for a matrix with +dimensions (3,3). + +matrix.sps:17: error: MATRIX: Matrix row index must be scalar or vector, not a +matrix with dimensions (2,2). + +matrix.sps:18: error: MATRIX: Matrix column index must be scalar or vector, not +a matrix with dimensions (2,2). + +x + 1 2 3 + 4 5 6 + 7 8 9 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - subvectors as rvalues]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT {10, 20, 30}({}). +PRINT {10, 20, 30}(2). +PRINT {10, 20, 30}({2}). +PRINT {10, 20, 30}({1,3}). +PRINT {10, 20, 30}({2,3}). +PRINT {10, 20, 30}({1;3}). +PRINT {10, 20, 30}({2;3}). +PRINT {10, 20, 30}(2:3). +PRINT {10, 20, 30}(:). + +PRINT {10; 20; 30}({}). +PRINT {10; 20; 30}(2). +PRINT {10; 20; 30}({2}). +PRINT {10; 20; 30}({1,3}). +PRINT {10; 20; 30}({2,3}). +PRINT {10; 20; 30}({1;3}). +PRINT {10; 20; 30}({2;3}). +PRINT {10; 20; 30}(2:3). +PRINT {10; 20; 30}(:). + +PRINT {}({}). + +PRINT {1, 2; 3, 4}(:). +PRINT {1, 2, 3, 4}({1, 2; 3, 4}). +PRINT {1, 2, 3, 4}(0). +PRINT {1, 2, 3, 4}(5). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +{10, 20, 30}({}) + +{10, 20, 30}(2) + 20 + +{10, 20, 30}({2}) + 20 + +{10, 20, 30}({1,3}) + 10 30 + +{10, 20, 30}({2,3}) + 20 30 + +{10, 20, 30}({1;3}) + 10 30 + +{10, 20, 30}({2;3}) + 20 30 + +{10, 20, 30}(2:3) + 20 30 + +{10, 20, 30}(:) + 10 20 30 + +{10; 20; 30}({}) + +{10; 20; 30}(2) + 20 + +{10; 20; 30}({2}) + 20 + +{10; 20; 30}({1,3}) + 10 + 30 + +{10; 20; 30}({2,3}) + 20 + 30 + +{10; 20; 30}({1;3}) + 10 + 30 + +{10; 20; 30}({2;3}) + 20 + 30 + +{10; 20; 30}(2:3) + 20 + 30 + +{10; 20; 30}(:) + 10 + 20 + 30 + +{}({}) + +matrix.sps:24: error: MATRIX: Vector index operator must be applied to vector, +not a matrix with dimensions (2,2). + +matrix.sps:25: error: MATRIX: Vector index must be scalar or vector, not a +matrix with dimensions (2,2). + +matrix.sps:26: error: MATRIX: Index 0 is out of range for vector with 4 +elements. + +matrix.sps:27: error: MATRIX: Index 5 is out of range for vector with 4 +elements. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - COMPUTE subvectors as lvalues]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE r={1, 2, 3, 4, 5, 6, 7, 8, 9}. + +COMPUTE r1=r. +COMPUTE r1(:) = {11, 12, 13, 14, 15, 16, 17, 18, 19}. +PRINT r1. + +COMPUTE r2=r. +COMPUTE r2(:) = {11; 12; 13; 14; 15; 16; 17; 18; 19}. +PRINT r2. + +COMPUTE r3=r. +COMPUTE r3(1) = 11. +PRINT r3. + +COMPUTE r4=r. +COMPUTE r4(1:2) = {11:12}. +PRINT r4. + +COMPUTE r5=r. +COMPUTE r5({8;9}) = {18:19}. +PRINT r5. + +COMPUTE c={1, 2, 3, 4, 5, 6, 7, 8, 9}. + +COMPUTE c1=c. +COMPUTE c1(:) = {11, 12, 13, 14, 15, 16, 17, 18, 19}. +PRINT c1. + +COMPUTE c2=c. +COMPUTE c2(:) = {11; 12; 13; 14; 15; 16; 17; 18; 19}. +PRINT c2. + +COMPUTE c3=c. +COMPUTE c3(1) = 11. +PRINT c3. + +COMPUTE c4=c. +COMPUTE c4(1:2) = {11:12}. +PRINT c4. + +COMPUTE c5=c. +COMPUTE c5(8:9) = {18:19}. +PRINT c5. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +r1 + 11 12 13 14 15 16 17 18 19 + +r2 + 11 12 13 14 15 16 17 18 19 + +r3 + 11 2 3 4 5 6 7 8 9 + +r4 + 11 12 3 4 5 6 7 8 9 + +r5 + 1 2 3 4 5 6 7 18 19 + +c1 + 11 12 13 14 15 16 17 18 19 + +c2 + 11 12 13 14 15 16 17 18 19 + +c3 + 11 2 3 4 5 6 7 8 9 + +c4 + 11 12 3 4 5 6 7 8 9 + +c5 + 1 2 3 4 5 6 7 18 19 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - COMPUTE subvectors as lvalues - negative]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE r={1, 2, 3, 4, 5, 6, 7, 8, 9}. +COMPUTE r(1:3) = {1, 2; 3, 4}. +COMPUTE r(1:3) = {}. +COMPUTE r(1:3) = {1}. +COMPUTE r(1:3) = {1, 2}. +COMPUTE r(1:3) = {1, 2, 3, 4}. +COMPUTE r(1:3) = {}. +COMPUTE r(1:3) = {1}. +COMPUTE r(1:3) = {1; 2}. +COMPUTE r(1:3) = {1; 2; 3; 4}. +COMPUTE r(:) = {1; 2; 3; 4}. +COMPUTE r(0) = 5. +COMPUTE r(10) = 5. +COMPUTE r({1, 2; 3, 4}) = 1. + +COMPUTE c={1, 2, 3, 4, 5, 6, 7, 8, 9}. +COMPUTE c(1:3) = {1, 2; 3, 4}. +COMPUTE c(1:3) = {}. +COMPUTE c(1:3) = {1}. +COMPUTE c(1:3) = {1, 2}. +COMPUTE c(1:3) = {1, 2, 3, 4}. +COMPUTE c(1:3) = {}. +COMPUTE c(1:3) = {1}. +COMPUTE c(1:3) = {1; 2}. +COMPUTE c(1:3) = {1; 2; 3; 4}. +COMPUTE c(:) = {1; 2; 3; 4}. +COMPUTE c(0) = 5. +COMPUTE c(10) = 5. +COMPUTE c({1, 2; 3, 4}) = 1. + +COMPUTE m = {1, 2; 3, 4}. +COMPUTE m(5) = 1. +COMPUTE m(:) = 1. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +matrix.sps:3: error: MATRIX: Can't assign matrix with dimensions (2,2) to +subvector. + +matrix.sps:4: error: MATRIX: Can't assign vector with 0 elements to subvector +with 3. + +matrix.sps:5: error: MATRIX: Can't assign vector with 1 elements to subvector +with 3. + +matrix.sps:6: error: MATRIX: Can't assign vector with 2 elements to subvector +with 3. + +matrix.sps:7: error: MATRIX: Can't assign vector with 4 elements to subvector +with 3. + +matrix.sps:8: error: MATRIX: Can't assign vector with 0 elements to subvector +with 3. + +matrix.sps:9: error: MATRIX: Can't assign vector with 1 elements to subvector +with 3. + +matrix.sps:10: error: MATRIX: Can't assign vector with 2 elements to subvector +with 3. + +matrix.sps:11: error: MATRIX: Can't assign vector with 4 elements to subvector +with 3. + +matrix.sps:12: error: MATRIX: Can't assign vector with 4 elements to subvector +with 9. + +matrix.sps:13: error: MATRIX: Index 0 is out of range for vector with 9 +elements. + +matrix.sps:14: error: MATRIX: Index 10 is out of range for vector with 9 +elements. + +matrix.sps:15: error: MATRIX: Vector index must be scalar or vector, not a +matrix with dimensions (2,2). + +matrix.sps:18: error: MATRIX: Can't assign matrix with dimensions (2,2) to +subvector. + +matrix.sps:19: error: MATRIX: Can't assign vector with 0 elements to subvector +with 3. + +matrix.sps:20: error: MATRIX: Can't assign vector with 1 elements to subvector +with 3. + +matrix.sps:21: error: MATRIX: Can't assign vector with 2 elements to subvector +with 3. + +matrix.sps:22: error: MATRIX: Can't assign vector with 4 elements to subvector +with 3. + +matrix.sps:23: error: MATRIX: Can't assign vector with 0 elements to subvector +with 3. + +matrix.sps:24: error: MATRIX: Can't assign vector with 1 elements to subvector +with 3. + +matrix.sps:25: error: MATRIX: Can't assign vector with 2 elements to subvector +with 3. + +matrix.sps:26: error: MATRIX: Can't assign vector with 4 elements to subvector +with 3. + +matrix.sps:27: error: MATRIX: Can't assign vector with 4 elements to subvector +with 9. + +matrix.sps:28: error: MATRIX: Index 0 is out of range for vector with 9 +elements. + +matrix.sps:29: error: MATRIX: Index 10 is out of range for vector with 9 +elements. + +matrix.sps:30: error: MATRIX: Vector index must be scalar or vector, not a +matrix with dimensions (2,2). + +matrix.sps:33: error: MATRIX: Can't use vector indexing on matrix m with +dimensions (2,2). + +matrix.sps:34: error: MATRIX: Can't use vector indexing on matrix m with +dimensions (2,2). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - COMPUTE - negative]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE x. +COMPUTE x=. +COMPUTE x(5)=1. +COMPUTE y(5)=1. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +matrix.sps:2.10: error: COMPUTE: Syntax error at end of command: expecting `='. + +matrix.sps:3.11: error: COMPUTE: Syntax error at end of command. + +matrix.sps:4: error: MATRIX: Undefined variable x. + +matrix.sps:5: error: COMPUTE: Undefined variable y. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - elementwise arithmetic operators]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT (-(5)). +PRINT (-{1,2;3,4}). + +PRINT ({1,2;3,4} + {5,6;7,8}). +PRINT ({1,2;3,4} + 5). +PRINT (5 + {5,6;7,8}). +PRINT ({1,2;3,4} + {5,6}). + +PRINT ({1,2;3,4} - {5,6;7,8}). +PRINT ({1,2;3,4} - 5). +PRINT (5 - {5,6;7,8}). +PRINT ({1,2;3,4} - {5,6}). + +PRINT ({1,2;3,4} * 5). +PRINT (5 * {5,6;7,8}). + +PRINT ({2,4;6,8} / 2). +PRINT (12 / {1,2;3,4}). +PRINT ({2,8;18,32} / {1,2;3,4}). + +PRINT ({1,2;3,4} &* {5,6;7,8}). +PRINT ({1,2;3,4} &* 5). +PRINT (5 &* {5,6;7,8}). +PRINT ({1,2;3,4} &* {5,6}). + +PRINT ({2,4;6,8} &/ 2). +PRINT (12 &/ {1,2;3,4}). +PRINT ({2,8;18,32} &/ {1,2;3,4}). + +PRINT ({1,2;3,4} &** 2). +PRINT (2 &** {1,2;3,4}). +PRINT ({1,2;3,4} &** {2,3;4,5}). +PRINT ({1,2;3,4} &** {5,6}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +(-(5)) + -5 + +(-{1,2;3,4}) + -1 -2 + -3 -4 + +({1,2;3,4} + {5,6;7,8}) + 6 8 + 10 12 + +({1,2;3,4} + 5) + 6 7 + 8 9 + +(5 + {5,6;7,8}) + 10 11 + 12 13 + +matrix.sps:8: error: MATRIX: Operands to + must have the same dimensions or one +must be a scalar, not matrices with dimensions (2,2) and (1,2). + +({1,2;3,4} - {5,6;7,8}) + -4 -4 + -4 -4 + +({1,2;3,4} - 5) + -4 -3 + -2 -1 + +(5 - {5,6;7,8}) + 0 -1 + -2 -3 + +matrix.sps:13: error: MATRIX: Operands to - must have the same dimensions or +one must be a scalar, not matrices with dimensions (2,2) and (1,2). + +({1,2;3,4} * 5) + 5 10 + 15 20 + +(5 * {5,6;7,8}) + 25 30 + 35 40 + +({2,4;6,8} / 2) + 1 2 + 3 4 + +(12 / {1,2;3,4}) + 12 6 + 4 3 + +({2,8;18,32} / {1,2;3,4}) + 2 4 + 6 8 + +({1,2;3,4} &* {5,6;7,8}) + 5 12 + 21 32 + +({1,2;3,4} &* 5) + 5 10 + 15 20 + +(5 &* {5,6;7,8}) + 25 30 + 35 40 + +matrix.sps:25: error: MATRIX: Operands to &* must have the same dimensions or +one must be a scalar, not matrices with dimensions (2,2) and (1,2). + +({2,4;6,8} &/ 2) + 1 2 + 3 4 + +(12 &/ {1,2;3,4}) + 12 6 + 4 3 + +({2,8;18,32} &/ {1,2;3,4}) + 2 4 + 6 8 + +({1,2;3,4} &** 2) + 1 4 + 9 16 + +(2 &** {1,2;3,4}) + 2 4 + 8 16 + +({1,2;3,4} &** {2,3;4,5}) + 1 8 + 81 1024 + +matrix.sps:34: error: MATRIX: Operands to &** must have the same dimensions or +one must be a scalar, not matrices with dimensions (2,2) and (1,2). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - relational operators]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT ({1, 1; 2, 2} > {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} > 1). +PRINT (2 > {1, 2; 1, 2}). +PRINT ({1, 2} > {1; 2}). + +PRINT ({1, 1; 2, 2} < {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} < 2). +PRINT (1 < {1, 2; 1, 2}). +PRINT ({1, 2} < {1; 2}). + +PRINT ({1, 1; 2, 2} <> {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} <> 2). +PRINT (1 <> {1, 2; 1, 2}). +PRINT ({1, 2} <> {1; 2}). + +PRINT ({1, 1; 2, 2} >= {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} >= 2). +PRINT (1 >= {1, 2; 1, 2}). +PRINT ({1, 2} >= {1; 2}). + +PRINT ({1, 1; 2, 2} <= {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} <= 2). +PRINT (1 <= {1, 2; 1, 2}). +PRINT ({1, 2} <= {1; 2}). + +PRINT ({1, 1; 2, 2} = {1, 2; 1, 2}). +PRINT ({1, 1; 2, 2} = 2). +PRINT (1 = {1, 2; 1, 2}). +PRINT ({1, 2} = {1; 2}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +({1, 1; 2, 2} > {1, 2; 1, 2}) + 0 0 + 1 0 + +({1, 1; 2, 2} > 1) + 0 0 + 1 1 + +(2 > {1, 2; 1, 2}) + 1 0 + 1 0 + +matrix.sps:5: error: MATRIX: Operands to > must have the same dimensions or one +must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({1, 1; 2, 2} < {1, 2; 1, 2}) + 0 1 + 0 0 + +({1, 1; 2, 2} < 2) + 1 1 + 0 0 + +(1 < {1, 2; 1, 2}) + 0 1 + 0 1 + +matrix.sps:10: error: MATRIX: Operands to < must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({1, 1; 2, 2} <> {1, 2; 1, 2}) + 0 1 + 1 0 + +({1, 1; 2, 2} <> 2) + 1 1 + 0 0 + +(1 <> {1, 2; 1, 2}) + 0 1 + 0 1 + +matrix.sps:15: error: MATRIX: Operands to <> must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({1, 1; 2, 2} >= {1, 2; 1, 2}) + 1 0 + 1 1 + +({1, 1; 2, 2} >= 2) + 0 0 + 1 1 + +(1 >= {1, 2; 1, 2}) + 1 0 + 1 0 + +matrix.sps:20: error: MATRIX: Operands to >= must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({1, 1; 2, 2} <= {1, 2; 1, 2}) + 1 1 + 0 1 + +({1, 1; 2, 2} <= 2) + 1 1 + 1 1 + +(1 <= {1, 2; 1, 2}) + 1 1 + 1 1 + +matrix.sps:25: error: MATRIX: Operands to <= must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({1, 1; 2, 2} = {1, 2; 1, 2}) + 1 0 + 0 1 + +({1, 1; 2, 2} = 2) + 0 0 + 1 1 + +(1 = {1, 2; 1, 2}) + 1 0 + 1 0 + +matrix.sps:30: error: MATRIX: Operands to = must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - logical operators]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT (NOT {-1, 0, 1}). + +PRINT ({-1, 0, 1; -1, 0, 1; -1, 0, 1} AND {-1, -1, -1; 0, 0, 0; 1, 1, 1}). +PRINT ({-1, 0, 1} AND -1). +PRINT ({-1, 0, 1} AND 0). +PRINT ({-1, 0, 1} AND 1). +PRINT ({-1, 0} AND {2; 3}). + +PRINT ({-1, 0, 1; -1, 0, 1; -1, 0, 1} OR {-1, -1, -1; 0, 0, 0; 1, 1, 1}). +PRINT ({-1, 0, 1} OR -1). +PRINT ({-1, 0, 1} OR 0). +PRINT ({-1, 0, 1} OR 1). +PRINT ({-1, 0} OR {2; 3}). + +PRINT ({-1, 0, 1; -1, 0, 1; -1, 0, 1} XOR {-1, -1, -1; 0, 0, 0; 1, 1, 1}). +PRINT ({-1, 0, 1} XOR -1). +PRINT ({-1, 0, 1} XOR 0). +PRINT ({-1, 0, 1} XOR 1). +PRINT ({-1, 0} XOR {2; 3}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +(NOT {-1, 0, 1}) + 1 1 0 + +({-1, 0, 1; -1, 0, 1; -1, 0, 1} AND {-1, -1, -1; 0, 0, 0; 1, 1, 1}) + 0 0 0 + 0 0 0 + 0 0 1 + +({-1, 0, 1} AND -1) + 0 0 0 + +({-1, 0, 1} AND 0) + 0 0 0 + +({-1, 0, 1} AND 1) + 0 0 1 + +matrix.sps:8: error: MATRIX: Operands to AND must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({-1, 0, 1; -1, 0, 1; -1, 0, 1} OR {-1, -1, -1; 0, 0, 0; 1, 1, 1}) + 0 0 1 + 0 0 1 + 1 1 1 + +({-1, 0, 1} OR -1) + 0 0 1 + +({-1, 0, 1} OR 0) + 0 0 1 + +({-1, 0, 1} OR 1) + 1 1 1 + +matrix.sps:14: error: MATRIX: Operands to OR must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). + +({-1, 0, 1; -1, 0, 1; -1, 0, 1} XOR {-1, -1, -1; 0, 0, 0; 1, 1, 1}) + 0 0 1 + 0 0 1 + 1 1 0 + +({-1, 0, 1} XOR -1) + 0 0 1 + +({-1, 0, 1} XOR 0) + 0 0 1 + +({-1, 0, 1} XOR 1) + 1 1 0 + +matrix.sps:20: error: MATRIX: Operands to XOR must have the same dimensions or +one must be a scalar, not matrices with dimensions (1,2) and (2,1). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - matrix operators]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT ({0, 1; 0, 0} * {0, 0; 1, 0}). +PRINT ({0, 0; 1, 0} * {0, 1; 0, 0}). +PRINT ({1, 2, 3; 4, 5, 6} * {7, 8; 9, 10; 11, 12}). +PRINT ({3, 4, 2} * {13, 9, 7, 15; 8, 7, 4, 6; 6, 4, 0, 3}). +COMPUTE m = {0, 1, 0, 0; 1, 0, 1, 0; 0, 1, 0, 1; 0, 0, 1, 0}. +PRINT m**-2. +PRINT m**-1. +PRINT m**0. +PRINT m**1. +PRINT m**2. +PRINT m**3. +PRINT m**5. +PRINT {3, 3.5; 3.2, 3.6}**-1/FORMAT F6.2. + +PRINT ({1, 2, 3} * {1, 2}). +PRINT {1, 2, 3}**2. +PRINT m**{1, 2}. +PRINT m**1.5. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +({0, 1; 0, 0} * {0, 0; 1, 0}) + 1 0 + 0 0 + +({0, 0; 1, 0} * {0, 1; 0, 0}) + 0 0 + 0 1 + +({1, 2, 3; 4, 5, 6} * {7, 8; 9, 10; 11, 12}) + 58 64 + 139 154 + +({3, 4, 2} * {13, 9, 7, 15; 8, 7, 4, 6; 6, 4, 0, 3}) + 83 63 37 75 + +m**-2 + 2 0 -1 0 + 0 1 0 -1 + -1 0 1 0 + 0 -1 0 2 + +m**-1 + 0 1 0 -1 + 1 0 0 0 + 0 0 0 1 + -1 0 1 0 + +m**0 + 1 0 0 0 + 0 1 0 0 + 0 0 1 0 + 0 0 0 1 + +m**1 + 0 1 0 0 + 1 0 1 0 + 0 1 0 1 + 0 0 1 0 + +m**2 + 1 0 1 0 + 0 2 0 1 + 1 0 2 0 + 0 1 0 1 + +m**3 + 0 2 0 1 + 2 0 3 0 + 0 3 0 2 + 1 0 2 0 + +m**5 + 0 5 0 3 + 5 0 8 0 + 0 8 0 5 + 3 0 5 0 + +{3, 3.5; 3.2, 3.6}**-1 + -9.00 8.75 + 8.00 -7.50 + +matrix.sps:16: error: MATRIX: Matrices with dimensions (1,3) and (1,2) are not +conformable for multiplication. + +matrix.sps:17: error: MATRIX: Matrix exponentation with ** requires a square +matrix on the left-hand size, not one with dimensions (1,3). + +matrix.sps:18: error: MATRIX: Matrix exponentiation with ** requires a scalar +on the right-hand side, not a matrix with dimensions (1,2). + +matrix.sps:19: error: MATRIX: Exponent 1.5 in matrix multiplication is non- +integer or outside the valid range. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - sequences and construction]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT {1:3:-1}. +PRINT {1:3}. +PRINT {1:10:2}. +PRINT {1:11:2}. + +PRINT {-1:-3}. +PRINT {-1:-3:-1}. +PRINT {-1:-10:-2}. +PRINT {-1:-11:-2}. + +PRINT {1:3:0}. +PRINT {-1:-3:0}. + +PRINT {1, 2; 3}. +PRINT {{2; 5}, 3}. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +{1:3:-1} + +{1:3} + 1 2 3 + +{1:10:2} + 1 3 5 7 9 + +{1:11:2} + 1 3 5 7 9 11 + +{-1:-3} + +{-1:-3:-1} + -1 -2 -3 + +{-1:-10:-2} + -1 -3 -5 -7 -9 + +{-1:-11:-2} + -1 -3 -5 -7 -9 -11 + +matrix.sps:12: error: MATRIX: The increment operand to : must be nonzero. + +matrix.sps:13: error: MATRIX: The increment operand to : must be nonzero. + +matrix.sps:15: error: MATRIX: All rows in a matrix must have the same number of +columns, but this tries to stack matrices with 2 and 1 columns. + +matrix.sps:16: error: MATRIX: All columns in a matrix must have the same number +of rows, but this tries to paste matrices with 2 and 1 rows. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - comments]) +AT_DATA([matrix.sps], [dnl +MATRIX. +* Comment one. +PRINT (1+2). +COMMENT Comment two. +PRINT (3+4). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +(1+2) + 3 + +(3+4) + 7 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - string matrices]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE m={'This is', 'a string', 'matrix', 'including', 'some', 'long strings'}. +PRINT m/FORMAT=A8. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +m + This is a string matrix includin some long str +]) +AT_CLEANUP + +AT_SETUP([MATRIX - ABS ALL ANY ARSIN ARTAN]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT ABS({-1, 0, 1}). + +PRINT ALL({0, 0, 0}). +PRINT ALL({-1, 1}). +PRINT ALL({-1, 0, 1}). + +PRINT ANY({0, 0, 0}). +PRINT ANY({-1, 1}). +PRINT ANY({-1, 0, 1}). + +PRINT ARSIN({-1, 0, 1})/FORMAT=F5.2. + +PRINT ARTAN({-5, -1, 0, 1, 5})/FORMAT=F5.2. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +ABS({-1, 0, 1}) + 1 0 1 + +ALL({0, 0, 0}) + 0 + +ALL({-1, 1}) + 1 + +ALL({-1, 0, 1}) + 0 + +ANY({0, 0, 0}) + 0 + +ANY({-1, 1}) + 1 + +ANY({-1, 0, 1}) + 1 + +ARSIN({-1, 0, 1}) + -1.57 .00 1.57 + +ARTAN({-5, -1, 0, 1, 5}) + -1.37 -.79 .00 .79 1.37 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - BLOCK CHOL CMAX CMIN COS]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT BLOCK({1, 2; 3, 4}, 5, {7; 8; 9}, {10, 11}). + +COMPUTE b=CHOL({4, 12, -16; 12, 37, -43; -16, -43, 98}). +PRINT b. +PRINT (T(b)*b). + +PRINT CMAX({9, 3, 4; 5, 8, 6; 7, 4, 11}). + +PRINT CMIN({9, 3, 4; 5, 8, 6; 7, 4, 11}). + +PRINT COS({0.785, 1.57; 3.14, 1.57 + 3.14}) /FORMAT=F5.2. + +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +BLOCK({1, 2; 3, 4}, 5, {7; 8; 9}, {10, 11}) + 1 2 0 0 0 0 + 3 4 0 0 0 0 + 0 0 5 0 0 0 + 0 0 0 7 0 0 + 0 0 0 8 0 0 + 0 0 0 9 0 0 + 0 0 0 0 10 11 + +b + 2 6 -8 + 0 1 5 + 0 0 3 + +(T(b)*b) + 4 12 -16 + 12 37 -43 + -16 -43 98 + +CMAX({9, 3, 4; 5, 8, 6; 7, 4, 11}) + 9 8 11 + +CMIN({9, 3, 4; 5, 8, 6; 7, 4, 11}) + 5 3 4 + +COS({0.785, 1.57; 3.14, 1.57 + 3.14}) + .71 .00 + -1.00 .00 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - CSSQ CSUM DESIGN DET DIAG]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT CSSQ({1, 2, 3; 4, 5, 6; 7, 8, 9}). +PRINT CSUM({1, 2, 3; 4, 5, 6; 7, 8, 9}). +PRINT DESIGN({1, 2, 0; 2, 1, 0; 3, 0, 1}). +PRINT DESIGN({1, 2, 0; 2, 2, 0; 3, 2, 1}). +PRINT DET({1, 2, 3; 4, 5, 6; 7, 8, 9}) /FORMAT F4.1. +PRINT DIAG({1, 2, 3, 4; 4, 5, 6, 7; 7, 8, 9, 10}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +CSSQ({1, 2, 3; 4, 5, 6; 7, 8, 9}) + 66 93 126 + +CSUM({1, 2, 3; 4, 5, 6; 7, 8, 9}) + 12 15 18 + +DESIGN({1, 2, 0; 2, 1, 0; 3, 0, 1}) + 1 0 0 0 0 1 1 0 + 0 1 0 0 1 0 1 0 + 0 0 1 1 0 0 0 1 + +warning: Column 2 in DESIGN argument has constant value. + +DESIGN({1, 2, 0; 2, 2, 0; 3, 2, 1}) + 1 0 0 1 0 + 0 1 0 1 0 + 0 0 1 0 1 + +DET({1, 2, 3; 4, 5, 6; 7, 8, 9}) + .0 + +DIAG({1, 2, 3, 4; 4, 5, 6, 7; 7, 8, 9, 10}) + 1 + 5 + 9 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - EVAL EXP GINV GRADE GSCH]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT EVAL({2, 0, 0; 0, 3, 4; 0, 4, 9}). + +PRINT EXP({2, 3; 4, 5})/FORMAT F5.2. + +PRINT GINV({1, 2})/FORMAT F5.2. +COMPUTE a={1, 2, 3; 4, 5, 6; 7, 8, 9}. +COMPUTE g=GINV(a). +PRINT (a*g*a)/FORMAT F5.2. + +PRINT GRADE({1, 0, 3; 3, 1, 2; 3, 0, 5}). +COMPUTE x={26, 690, 323, 208, 671, 818, 732, 711, 585, 792}. +COMPUTE asort=x. +COMPUTE asort(GRADE(asort))=asort. +PRINT asort. +COMPUTE dsort=x. +COMPUTE dsort(GRADE(-dsort))=dsort. +PRINT dsort. + +PRINT (GSCH({3, 2; 1, 2}) * SQRT(10))/FORMAT F5.2. +PRINT (GSCH({0, 3, 6, 2; 0, 1, 2, 2}) * SQRT(10))/FORMAT F5.2. +PRINT GSCH({0; 0}). +PRINT GSCH({0, 0, 0; 0, 0, 0}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +EVAL({2, 0, 0; 0, 3, 4; 0, 4, 9}) + 11.0000000000 + 2.0000000000 + 1.0000000000 + +EXP({2, 3; 4, 5}) + 7.39 20.09 + 54.60 148.4 + +GINV({1, 2}) + .20 + .40 + +(a*g*a) + 1.00 2.00 3.00 + 4.00 5.00 6.00 + 7.00 8.00 9.00 + +GRADE({1, 0, 3; 3, 1, 2; 3, 0, 5}) + 3 1 6 + 7 4 5 + 8 2 9 + +asort + 26 208 323 585 671 690 711 732 792 818 + +dsort + 818 792 732 711 690 671 585 323 208 26 + +(GSCH({3, 2; 1, 2}) * SQRT(10)) + 3.00 -1.00 + 1.00 3.00 + +(GSCH({0, 3, 6, 2; 0, 1, 2, 2}) * SQRT(10)) + 3.00 -1.00 + 1.00 3.00 + +matrix.sps:22: error: MATRIX: GSCH requires its argument to have at least as +many columns as rows, but it has dimensions (2,1). + +matrix.sps:23: error: MATRIX: Argument to GSCH with dimensions (2,3) contains +only 0 linearly independent columns. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - IDENT INV KRONEKER LG10 LN]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT IDENT(1). +PRINT IDENT(2). +PRINT IDENT(3,5). +PRINT IDENT(5,3). + +PRINT INV({3, 3.5; 3.2, 3.6})/FORMAT F8.2. +PRINT INV({4, 7; 2, 6})/FORMAT F8.2. +PRINT (INV({4, -2, 1; 5, 0, 3; -1, 2, 6})*52)/FORMAT F8.2. + +PRINT KRONEKER({1, 2; 3, 4}, {0, 5; 6, 7}). + +PRINT LG10({1, 10, 100, 1000}). + +PRINT LN({1, 2; 3, 4})/FORMAT F5.2. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +IDENT(1) + 1 + +IDENT(2) + 1 0 + 0 1 + +IDENT(3,5) + 1 0 0 0 0 + 0 1 0 0 0 + 0 0 1 0 0 + +IDENT(5,3) + 1 0 0 + 0 1 0 + 0 0 1 + 0 0 0 + 0 0 0 + +INV({3, 3.5; 3.2, 3.6}) + -9.00 8.75 + 8.00 -7.50 + +INV({4, 7; 2, 6}) + .60 -.70 + -.20 .40 + +(INV({4, -2, 1; 5, 0, 3; -1, 2, 6})*52) + -6.00 14.00 -6.00 + -33.00 25.00 -7.00 + 10.00 -6.00 10.00 + +KRONEKER({1, 2; 3, 4}, {0, 5; 6, 7}) + 0 5 0 10 + 6 7 12 14 + 0 15 0 20 + 18 21 24 28 + +LG10({1, 10, 100, 1000}) + 0 1 2 3 + +LN({1, 2; 3, 4}) + .00 .69 + 1.10 1.39 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - MAGIC]) +AT_DATA([matrix.sps], [dnl +MATRIX. + +LOOP n=3 to 10. + COMPUTE m=MAGIC(n). + COMPUTE total=n*(n**2 + 1) / 2. + COMPUTE tb={MSUM(DIAG(T(m))), CSUM(m), MSUM(DIAG(m))} - total. + COMPUTE lr=RSUM(m) - total. + PRINT {tb; lr, m, lr; tb}/FORMAT F4.0. +END LOOP. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +{tb; lr, m, lr; tb} + 0 0 0 0 0 + 0 8 1 6 0 + 0 3 5 7 0 + 0 4 9 2 0 + 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 + 0 1 5 12 16 0 + 0 15 11 6 2 0 + 0 14 8 9 3 0 + 0 4 10 7 13 0 + 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 + 0 17 24 1 8 15 0 + 0 23 5 7 14 16 0 + 0 4 6 13 20 22 0 + 0 10 12 19 21 3 0 + 0 11 18 25 2 9 0 + 0 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 0 + 0 1 5 9 28 32 36 0 + 0 35 30 27 10 7 2 0 + 0 24 14 22 18 17 16 0 + 0 13 23 15 19 20 21 0 + 0 34 31 26 11 6 3 0 + 0 4 8 12 25 29 33 0 + 0 0 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 0 0 + 0 30 39 48 1 10 19 28 0 + 0 38 47 7 9 18 27 29 0 + 0 46 6 8 17 26 35 37 0 + 0 5 14 16 25 34 36 45 0 + 0 13 15 24 33 42 44 4 0 + 0 21 23 32 41 43 3 12 0 + 0 22 31 40 49 2 11 20 0 + 0 0 0 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 0 0 0 + 0 1 9 17 25 40 48 56 64 0 + 0 63 55 47 39 26 18 10 2 0 + 0 3 11 19 27 38 46 54 62 0 + 0 61 53 45 37 28 20 12 4 0 + 0 60 52 44 32 33 21 13 5 0 + 0 6 14 22 30 35 43 51 59 0 + 0 58 50 42 34 31 23 15 7 0 + 0 8 16 24 36 29 41 49 57 0 + 0 0 0 0 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 0 0 0 0 + 0 47 58 69 80 1 12 23 34 45 0 + 0 57 68 79 9 11 22 33 44 46 0 + 0 67 78 8 10 21 32 43 54 56 0 + 0 77 7 18 20 31 42 53 55 66 0 + 0 6 17 19 30 41 52 63 65 76 0 + 0 16 27 29 40 51 62 64 75 5 0 + 0 26 28 39 50 61 72 74 4 15 0 + 0 36 38 49 60 71 73 3 14 25 0 + 0 37 48 59 70 81 2 13 24 35 0 + 0 0 0 0 0 0 0 0 0 0 0 +{tb; lr, m, lr; tb} + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 1 9 17 25 33 68 76 84 92 100 0 + 0 99 91 83 75 67 34 26 18 10 2 0 + 0 3 11 19 27 35 66 74 82 90 98 0 + 0 97 89 81 72 65 36 29 20 12 4 0 + 0 60 42 58 44 56 50 49 53 47 46 0 + 0 41 59 43 57 45 51 52 48 54 55 0 + 0 96 88 80 73 64 37 28 21 13 5 0 + 0 6 14 22 30 38 63 71 79 87 95 0 + 0 94 86 78 70 62 39 31 23 15 7 0 + 0 8 16 24 32 40 61 69 77 85 93 0 + 0 0 0 0 0 0 0 0 0 0 0 0 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - MAKE MDIAG MMAX MMIN MOD]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT MAKE(1, 2, 3). +PRINT MAKE(2, 1, 4). +PRINT MAKE(2, 3, 5). + +PRINT MDIAG({1, 2, 3, 4}). +PRINT MDIAG({1; 2; 3; 4}). +PRINT MDIAG({1, 2; 3, 4}). + +PRINT MMAX({55, 44; 66, 11}). + +PRINT MMIN({55, 44; 66, 11}). + +PRINT MOD({5, 4, 3, 2, 1, 0}, 3). +PRINT MOD({5, 4, 3, 2, 1, 0}, -3). +PRINT MOD({-5, -4, -3, -2, -1, 0}, 3). +PRINT MOD({-5, -4, -3, -2, -1, 0}, -3). +PRINT MOD({5, 4, 3, 2, 1, 0}, 1.5) /FORMAT F5.1. +PRINT MOD({5, 4, 3, 2, 1, 0}, 0). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +MAKE(1, 2, 3) + 3 3 + +MAKE(2, 1, 4) + 4 + 4 + +MAKE(2, 3, 5) + 5 5 5 + 5 5 5 + +MDIAG({1, 2, 3, 4}) + 1 0 0 0 + 0 2 0 0 + 0 0 3 0 + 0 0 0 4 + +MDIAG({1; 2; 3; 4}) + 1 0 0 0 + 0 2 0 0 + 0 0 3 0 + 0 0 0 4 + +matrix.sps:8: error: MATRIX: Function MDIAG argument 1 must be a vector, but it +has dimensions (2,2). + +MMAX({55, 44; 66, 11}) + 66 + +MMIN({55, 44; 66, 11}) + 11 + +MOD({5, 4, 3, 2, 1, 0}, 3) + 2 1 0 2 1 0 + +MOD({5, 4, 3, 2, 1, 0}, -3) + 2 1 0 2 1 0 + +MOD({-5, -4, -3, -2, -1, 0}, 3) + -2 -1 0 -2 -1 0 + +MOD({-5, -4, -3, -2, -1, 0}, -3) + -2 -1 0 -2 -1 0 + +MOD({5, 4, 3, 2, 1, 0}, 1.5) + .5 1.0 .0 .5 1.0 .0 + +matrix.sps:19: error: MATRIX: Divisor argument to MOD function must be nonzero. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - MSSQ MSUM NCOL NROW RANK]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT MSSQ({1, 0, 1; -2, -3, 1; 3, 3, 0}). + +PRINT MSUM({1, 0, 1; -2, -3, 1; 3, 3, 0}). + +PRINT NCOL({1, 0; -2, -3; 3, 3}). + +PRINT NROW({1, 0; -2, -3; 3, 3}). + +PRINT RANK({1, 0, 1; -2, -3, 1; 3, 3, 0}). +PRINT RANK({1, 1, 0, 2; -1, -1, 0, -2}). +PRINT RANK({1, -1; 1, -1; 0, 0; 2, -2}). +PRINT RANK({1, 2, 1; -2, -3, 1; 3, 5, 0}). +PRINT RANK({1, 0, 2; 2, 1, 0; 3, 2, 1}). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +MSSQ({1, 0, 1; -2, -3, 1; 3, 3, 0}) + 34 + +MSUM({1, 0, 1; -2, -3, 1; 3, 3, 0}) + 4 + +NCOL({1, 0; -2, -3; 3, 3}) + 2 + +NROW({1, 0; -2, -3; 3, 3}) + 3 + +RANK({1, 0, 1; -2, -3, 1; 3, 3, 0}) + 2 + +RANK({1, 1, 0, 2; -1, -1, 0, -2}) + 1 + +RANK({1, -1; 1, -1; 0, 0; 2, -2}) + 1 + +RANK({1, 2, 1; -2, -3, 1; 3, 5, 0}) + 2 + +RANK({1, 0, 2; 2, 1, 0; 3, 2, 1}) + 3 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - RESHAPE RMAX RMIN RND RNKORDER]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 1, 12). +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 2, 6). +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 3, 4). +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 4, 3). +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 6, 2). +PRINT RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 12, 1). + +PRINT RMAX({1, 0, 1; -2, -3, 1; 3, 3, 0}). + +PRINT RMIN({1, 0, 1; -2, -3, 1; 3, 3, 0}). + +PRINT RND({-1.6, -1.5, -1.4; + -.6, -.5, -.4; + .4, .5, .6; + 1.4, 1.5, 1.6})/FORMAT F5.1. + +PRINT RNKORDER({1, 0, 3; 3, 1, 2; 3, 0, 5}) /FORMAT F5.1. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 1, 12) + 1 2 3 4 5 6 7 8 9 10 11 12 + +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 2, 6) + 1 2 3 4 5 6 + 7 8 9 10 11 12 + +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 3, 4) + 1 2 3 4 + 5 6 7 8 + 9 10 11 12 + +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 4, 3) + 1 2 3 + 4 5 6 + 7 8 9 + 10 11 12 + +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 6, 2) + 1 2 + 3 4 + 5 6 + 7 8 + 9 10 + 11 12 + +RESHAPE({1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12}, 12, 1) + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + +RMAX({1, 0, 1; -2, -3, 1; 3, 3, 0}) + 1 + 1 + 3 + +RMIN({1, 0, 1; -2, -3, 1; 3, 3, 0}) + 0 + -3 + 0 + +RND({-1.6, -1.5, -1.4; + -.6, -.5, -.4; + .4, .5, .6; + 1.4, 1.5, 1.6}) + -2.0 -2.0 -1.0 + -1.0 .0 .0 + .0 .0 1.0 + 1.0 2.0 2.0 + +RNKORDER({1, 0, 3; 3, 1, 2; 3, 0, 5}) + 3.5 1.5 7.0 + 7.0 3.5 5.0 + 7.0 1.5 9.0 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - RSSQ RSUM SIN SOLVE SQRT]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT RSSQ({1, 2, 3; 4, 5, 6; 7, 8, 9}). +PRINT RSUM({1, 2, 3; 4, 5, 6; 7, 8, 9}). + +PRINT SIN({0, .78, 1.57, 2.35, 3.14}) /FORMAT F5.2. + +PRINT SOLVE({2, 3; 4, 9}, {6, 2; 15, 5}) /FORMAT=F6.2. +PRINT SOLVE({1, 3, -2; 3, 5, 6; 2, 4, 3}, {5; 7; 8}) /FORMAT=F6.2. +PRINT SOLVE({2, 1, -1; -3, -1, 2; -2, 1, 2}, {8; -11; -3}) /FORMAT=F6.2. +PRINT SOLVE({1, 2; 3, 4}, {1, 2}). + +PRINT SQRT({0, 1, 2, 3, 4, 9, 81}) /FORMAT=F5.2. +PRINT SQRT(-1). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +RSSQ({1, 2, 3; 4, 5, 6; 7, 8, 9}) + 14 + 77 + 194 + +RSUM({1, 2, 3; 4, 5, 6; 7, 8, 9}) + 6 + 15 + 24 + +SIN({0, .78, 1.57, 2.35, 3.14}) + .00 .70 1.00 .71 .00 + +SOLVE({2, 3; 4, 9}, {6, 2; 15, 5}) + 1.50 .50 + 1.00 .33 + +SOLVE({1, 3, -2; 3, 5, 6; 2, 4, 3}, {5; 7; 8}) + -15.00 + 8.00 + 2.00 + +SOLVE({2, 1, -1; -3, -1, 2; -2, 1, 2}, {8; -11; -3}) + 2.00 + 3.00 + -1.00 + +matrix.sps:10: error: MATRIX: SOLVE requires its arguments to have the same +number of rows, but the first argument has dimensions (2,2) and the second +(1,2). + +SQRT({0, 1, 2, 3, 4, 9, 81}) + .00 1.00 1.41 1.73 2.00 3.00 9.00 + +matrix.sps:13: error: MATRIX: Argument to SQRT must be nonnegative. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - SSCP SVAL SWEEP TRACE TRANSPOS TRUNC]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE m={1, 2, 3; 4, 5, 6} +COMPUTE sscp1=SSCP(m). +COMPUTE sscp2=T(m)*m. +PRINT sscp1. +PRINT (sscp1 <> sscp2). + +PRINT SVAL({1, 1; 0, 0})/FORMAT F5.2. +PRINT SVAL({1, 0, 1; 0, 1, 1; 0, 0, 0})/FORMAT F5.2. +PRINT SVAL({1, 0, 0, 0, 2; 0, 0, 3, 0, 0; 0, 0, 0, 0, 0; 0, 2, 0, 0, 0}) + /FORMAT F5.2. +PRINT SVAL({2, 4; 1, 3; 0, 0; 0, 0})/FORMAT F5.2. + +COMPUTE s0 = {6, 12, 0, 12; 12, 28, 0, 25; 0, 0, 6, 2; 12, 25, 2, 28}. +PRINT SWEEP(s0, 1)/FORMAT F5.2. +PRINT SWEEP(SWEEP(s0, 1), 2)/FORMAT F5.2. +PRINT SWEEP(SWEEP(SWEEP(s0, 1), 2), 3)/FORMAT F5.2. + +COMPUTE s1 = {6, 12, 0, 12; 12, 0, 0, 25; 0, 0, 6, 2; 12, 25, 2, 28}. +PRINT SWEEP(s1, 2). + +PRINT TRACE(s0). + +PRINT T(s0). +PRINT TRANSPOS(s0). +PRINT ALL(T(T(s0)) = s0). + +PRINT TRUNC(SVAL({2, 4; 1, 3; 0, 0; 0, 0})). +PRINT TRUNC(-SVAL({2, 4; 1, 3; 0, 0; 0, 0})). +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +sscp1 + 17 22 27 + 22 29 36 + 27 36 45 + +(sscp1 <> sscp2) + 0 0 0 + 0 0 0 + 0 0 0 + +SVAL({1, 1; 0, 0}) + 1.41 + .00 + +SVAL({1, 0, 1; 0, 1, 1; 0, 0, 0}) + 1.73 + 1.00 + .00 + +SVAL({1, 0, 0, 0, 2; 0, 0, 3, 0, 0; 0, 0, 0, 0, 0; 0, 2, 0, 0, 0}) + 3.00 + 2.24 + 2.00 + .00 + +SVAL({2, 4; 1, 3; 0, 0; 0, 0}) + 5.46 + .37 + +SWEEP(s0, 1) + .17 2.00 .00 2.00 + -2.00 4.00 .00 1.00 + .00 .00 6.00 2.00 + -2.00 1.00 2.00 4.00 + +SWEEP(SWEEP(s0, 1), 2) + 1.17 -.50 .00 1.50 + -.50 .25 .00 .25 + .00 .00 6.00 2.00 + -1.50 -.25 2.00 3.75 + +SWEEP(SWEEP(SWEEP(s0, 1), 2), 3) + 1.17 -.50 .00 1.50 + -.50 .25 .00 .25 + .00 .00 .17 .33 + -1.50 -.25 -.33 3.08 + +SWEEP(s1, 2) + 6 0 0 12 + 0 0 0 0 + 0 0 6 2 + 12 0 2 28 + +TRACE(s0) + 68 + +T(s0) + 6 12 0 12 + 12 28 0 25 + 0 0 6 2 + 12 25 2 28 + +TRANSPOS(s0) + 6 12 0 12 + 12 28 0 25 + 0 0 6 2 + 12 25 2 28 + +ALL(T(T(s0)) = s0) + 1 + +TRUNC(SVAL({2, 4; 1, 3; 0, 0; 0, 0})) + 5 + 0 + +TRUNC(-SVAL({2, 4; 1, 3; 0, 0; 0, 0})) + -5 + 0 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - UNIFORM]) +AT_DATA([matrix.sps], [dnl +SET SEED=10. +MATRIX. +PRINT (UNIFORM(4, 5)*10)/FORMAT F5.2. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +(UNIFORM(4, 5)*10) + 7.71 2.99 .21 4.95 6.34 + 4.43 7.49 8.32 4.99 5.83 + 2.25 .25 1.98 7.09 7.61 + 2.66 1.69 2.64 .88 1.50 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - CALL SETDIAG]) +AT_DATA([matrix.sps], [dnl +MATRIX. +COMPUTE x={1, 2, 3; 4, 5, 6; 7, 8, 9}. + +COMPUTE x1=x. +CALL SETDIAG(x1, 10). +PRINT x1. + +COMPUTE x2=x. +CALL SETDIAG(x2, {10, 11}). +PRINT x2. + +COMPUTE x3=x. +CALL SETDIAG(x3, {10, 11, 12}). +PRINT x3. + +COMPUTE x4=x. +CALL SETDIAG(x4, {10, 11, 12, 13}). +PRINT x4. + +COMPUTE x5=x. +CALL SETDIAG(x5, {10, 11; 12, 13}). +PRINT x5. + +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +x1 + 10 2 3 + 4 10 6 + 7 8 10 + +x2 + 10 2 3 + 4 11 6 + 7 8 9 + +x3 + 10 2 3 + 4 11 6 + 7 8 12 + +x4 + 10 2 3 + 4 11 6 + 7 8 12 + +matrix.sps:21: error: MATRIX: SETDIAG argument 2 must be a scalar or a vector +but it has dimensions (2,2). + +x5 + 1 2 3 + 4 5 6 + 7 8 9 +]) +AT_CLEANUP + +dnl I have some doubts about the correctness of the results below. +AT_SETUP([MATRIX - CALL EIGEN]) +AT_DATA([matrix.sps], [dnl +MATRIX. +CALL EIGEN({1, 0; 0, 1}, evec, eval). +PRINT evec. +PRINT eval. + +CALL EIGEN({3, 2, 4; 2, 0, 2; 4, 2, 3}, evec, eval). +PRINT evec. +PRINT eval. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +evec + 1 0 + 0 1 + +eval + 1 + 1 + +evec + -.6666666667 .0000000000 .7453559925 + -.3333333333 -.8944271910 -.2981423970 + -.6666666667 .4472135955 -.5962847940 + +eval + 8.0000000000 + -1.0000000000 + -1.0000000000 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - CALL SVD]) +AT_DATA([matrix.sps], [dnl +MATRIX. +CALL SVD({3, 2, 2; 2, 3, -2}, u, s, v). +PRINT (u * s * T(v))/FORMAT F5.1. + +CALL SVD({2, 4; 1, 3; 0, 0; 0, 0}, u, s, v). +PRINT (u*s*T(v))/FORMAT F5.1. + +CALL SVD({-3, 1; 6, -2; 6, -2}, u, s, v). +PRINT (u*s*T(v))/FORMAT F5.1. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +(u * s * T(v)) + 3.0 2.0 2.0 + 2.0 3.0 -2.0 + +(u*s*T(v)) + 2.0 4.0 + 1.0 3.0 + .0 .0 + .0 .0 + +(u*s*T(v)) + -3.0 1.0 + 6.0 -2.0 + 6.0 -2.0 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - PRINT]) +AT_DATA([matrix.sps], [dnl +MATRIX. +PRINT/TITLE="title 1". +PRINT/SPACE=2/TITLE="title 2". + +COMPUTE m={1, 2, 3; 3, 4, 5; 6, 7, 8}. +PRINT m/RLABELS=123, a b c, long name. +PRINT m/RNAMES={'123', 'a b c', 'long name'}. +PRINT m/CLABELS=col1, col2, long column name. +PRINT m/CNAMES={'col1', 'col2', 'long column name'}. +PRINT m/RLABELS=123, a b c, long name + /CLABELS=col1, col2, long column name. +PRINT m/RNAMES={'123', 'a b c', 'long name'} + /CNAMES={'col1', 'col2', 'long column name'}. +PRINT {123e10, 456e10, 500}. +END MATRIX. +]) + +AT_DATA([matrix-tables.sps], [dnl +SET MDISPLAY=TABLES. +INCLUDE 'matrix.sps'. +]) + +AT_CHECK([pspp matrix.sps], [0], [dnl +title 1 + + + +title 2 + +m +123 1 2 3 +a b c 3 4 5 +long nam 6 7 8 + +m +123 1 2 3 +a b c 3 4 5 +long nam 6 7 8 + +m + col1 col2 long col + 1 2 3 + 3 4 5 + 6 7 8 + +m + col1 col2 long col + 1 2 3 + 3 4 5 + 6 7 8 + +m + col1 col2 long col +123 1 2 3 +a b c 3 4 5 +long nam 6 7 8 + +m + col1 col2 long col +123 1 2 3 +a b c 3 4 5 +long nam 6 7 8 + +{123e10, 456e10, 500} + 10 ** 12 X + 1.2300000000 4.5600000000 .0000000005 +]) + +AT_CHECK([pspp matrix-tables.sps], [0], [dnl +title 1 + + + +title 2 + + m ++---------+-----+ +|123 |1 2 3| +|a b c |3 4 5| +|long name|6 7 8| ++---------+-----+ + + m ++--------+-----+ +|123 |1 2 3| +|a b c |3 4 5| +|long nam|6 7 8| ++--------+-----+ + + m ++----+----+----------------+ +|col1|col2|long column name| ++----+----+----------------+ +| 1| 2| 3| +| 3| 4| 5| +| 6| 7| 8| ++----+----+----------------+ + + m ++----+----+--------+ +|col1|col2|long col| ++----+----+--------+ +| 1| 2| 3| +| 3| 4| 5| +| 6| 7| 8| ++----+----+--------+ + + m ++---------+----+----+----------------+ +| |col1|col2|long column name| ++---------+----+----+----------------+ +|123 | 1| 2| 3| +|a b c | 3| 4| 5| +|long name| 6| 7| 8| ++---------+----+----+----------------+ + + m ++--------+----+----+--------+ +| |col1|col2|long col| ++--------+----+----+--------+ +|123 | 1| 2| 3| +|a b c | 3| 4| 5| +|long nam| 6| 7| 8| ++--------+----+----+--------+ + + {123e10, 456e10, 500} ++----------------------------------------------+ +|1.2300000000[[a]] 4.5600000000[[a]] .0000000005[[a]]| ++----------------------------------------------+ +a. × 10**12 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - DO IF]) +AT_DATA([matrix.sps], [dnl +MATRIX. +DO IF 1. +PRINT/TITLE '1'. +END IF. + +DO IF 0. +PRINT/TITLE '2'. +ELSE IF 1. +PRINT/TITLE '3'. +END IF. + +DO IF -1. +PRINT/TITLE '4'. +ELSE IF 0. +PRINT/TITLE '5'. +ELSE. +PRINT/TITLE '6'. +END IF. + +DO IF {1, 2}. +END IF. + +DO IF 0. +ELSE IF {}. +END IF. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +1 + +3 + +6 + +matrix.sps:21: error: MATRIX: Expression for DO IF must evaluate to scalar, not +a matrix with dimensions (1,2). + +matrix.sps:25: error: MATRIX: Expression for ELSE IF must evaluate to scalar, +not a matrix with dimensions (0,0). +]) +AT_CLEANUP + +AT_SETUP([MATRIX - unbounded LOOP]) +AT_DATA([matrix.sps], [dnl +MATRIX. +* Truly unbounded loop. +COMPUTE x=0. +COMPUTE y={}. +LOOP. +COMPUTE x=x+1. +COMPUTE y={y, x}. +END LOOP. +PRINT x. +PRINT y. + +* Unbounded loop terminates with BREAK. +COMPUTE x=0. +COMPUTE y={}. +LOOP. +COMPUTE x=x+1. +COMPUTE y={y, x}. +DO IF x >= 20. + BREAK. +END IF. +END LOOP. +PRINT x. +PRINT y. + +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +x + 40 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 +40 + +x + 20 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 +]) +AT_CLEANUP + +AT_SETUP([MATRIX - indexed or conditional LOOP]) +AT_DATA([matrix.sps], [dnl +MATRIX. +* Indexed loop terminates based on index. +COMPUTE y={}. +LOOP x=1 TO 20. +COMPUTE y={y, x}. +END LOOP. +PRINT x. +PRINT y. + +* Indexed loop terminates based on MXLOOPS. +COMPUTE y={}. +LOOP x=1 TO 50. +COMPUTE y={y, x}. +END LOOP. +PRINT x. +PRINT y. + +* Indexed loop terminates with BREAK. +COMPUTE y={}. +LOOP x=1 TO 50. +COMPUTE y={y, x}. +DO IF x >= 20. + BREAK. +END IF. +END LOOP. +PRINT x. +PRINT y. + +* Indexed loop terminates with top IF. +COMPUTE y={}. +LOOP x=1 TO 50 IF NCOL(y) < 15. +COMPUTE y={y, x}. +END LOOP. +PRINT x. +PRINT y. + +* Indexed loop terminates with bottom IF. +COMPUTE y={}. +LOOP x=1 TO 50. +COMPUTE y={y, x}. +END LOOP IF NCOL(y) >= 22. +PRINT x. +PRINT y. + +* Index behavior. +COMPUTE indexing={ + 1, 10, 1; + 1, 10, 2; + 1, 10, 3; + 1, 10, -1; + 1, 10, 0; + 10, 1, -1; + 10, 1, -2; + 10, 1, -3; + 10, 1, 1; + 10, 1, 0 +}. +LOOP i=1 TO NROW(indexing). + COMPUTE y={}. + LOOP j=indexing(i, 1) TO indexing(i, 2) BY indexing(i, 3). + COMPUTE y={y, j}. + END LOOP. + PRINT {indexing(i, :), y}. +END LOOP. + +LOOP i={} TO 5. +END LOOP. + +LOOP i=5 TO {}. +END LOOP. + +LOOP i=5 TO 8 BY {}. +END LOOP. + +LOOP IF {}. +END LOOP. + +LOOP. +END LOOP IF {}. + +LOOP i=1e100 to 1e200. +END LOOP. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +x + 20 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 + +x + 40 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 +40 + +x + 20 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 + +x + 16 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 + +x + 22 + +y + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 +20 21 22 + +{indexing(i, :), y} + 1 10 1 1 2 3 4 5 6 7 8 9 10 +{indexing(i, :), y} + 1 10 2 1 3 5 7 9 +{indexing(i, :), y} + 1 10 3 1 4 7 10 +{indexing(i, :), y} + 1 10 -1 +{indexing(i, :), y} + 1 10 0 +{indexing(i, :), y} + 10 1 -1 10 9 8 7 6 5 4 3 2 1 +{indexing(i, :), y} + 10 1 -2 10 8 6 4 2 +{indexing(i, :), y} + 10 1 -3 10 7 4 1 +{indexing(i, :), y} + 10 1 1 +{indexing(i, :), y} + 10 1 0 + +matrix.sps:67: error: MATRIX: Expression for LOOP must evaluate to scalar, not +a matrix with dimensions (0,0). + +matrix.sps:70: error: MATRIX: Expression for TO must evaluate to scalar, not a +matrix with dimensions (0,0). + +matrix.sps:73: error: MATRIX: Expression for BY must evaluate to scalar, not a +matrix with dimensions (0,0). + +matrix.sps:76: error: MATRIX: Expression for LOOP IF must evaluate to scalar, +not a matrix with dimensions (0,0). + +matrix.sps:79: error: MATRIX: Expression for END LOOP IF must evaluate to +scalar, not a matrix with dimensions (0,0). + +matrix.sps:82: error: MATRIX: Expression for LOOP is outside the integer range. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - BREAK outside LOOP]) +AT_DATA([matrix.sps], [dnl +MATRIX. +BREAK. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [1], [dnl +matrix.sps:2: error: BREAK: BREAK not inside LOOP. +]) +AT_CLEANUP + +AT_SETUP([MATRIX - READ]) +AT_DATA([matrix.txt], [dnl +1 2 3 +4 5 6 +7 8 9 +10 11 12,13 +14, 15 ,16 , 17 +18 +19 +20 21 22 23 + 12 34 +5 6 + 78 89 +10 11 +]) +AT_DATA([matrix2.txt], [dnl +2, 3, 5, 7 +11, 13, 17, 19 +23, 29, 31, 37 +41, 43, 47, 53 +]) +AT_DATA([matrix.sps], [dnl +MATRIX. +READ x/FILE='matrix.txt'/SIZE={3,3}/FIELD=1 TO 80. +PRINT x. +READ x/SIZE={2,4}/FIELD=1 TO 80. +PRINT x. +READ x(:,2)/FILE='matrix.txt'/FIELD=1 TO 80. +PRINT x. +READ x(1,:)/SIZE={1,4}/FIELD=1 TO 80. +PRINT x. + +READ x/SIZE={2,6}/FIELD=1 TO 20 BY 5. +PRINT x. + +COMPUTE y={}. +LOOP IF NOT EOF('matrix2.txt'). +READ x/FILE='matrix2.txt'/SIZE={1,4}/FIELD=1 TO 80. +COMPUTE y={y; x}. +END LOOP. +PRINT y. +END MATRIX. +]) +AT_CHECK([pspp matrix.sps], [0], [dnl +x + 1 2 3 + 4 5 6 + 7 8 9 + +x + 10 11 12 13 + 14 15 16 17 + +x + 10 18 12 13 + 14 19 16 17 + +x + 20 21 22 23 + 14 19 16 17 + +x + 1 2 3 4 5 6 + 7 8 8 9 10 11 + +y + 2 3 5 7 + 11 13 17 19 + 23 29 31 37 + 41 43 47 53 +]) +AT_CLEANUP