From 1e194db8ad16820a84dcb0ca5c58fb8907f63ddb Mon Sep 17 00:00:00 2001 From: Ben Pfaff Date: Sun, 24 Oct 2021 21:37:06 -0700 Subject: [PATCH 1/1] More matrix examples. --- doc/matrices.texi | 190 ++++++++++++++++++++++++++++++--- src/language/stats/matrix.c | 9 +- tests/language/stats/matrix.at | 24 ++++- 3 files changed, 198 insertions(+), 25 deletions(-) diff --git a/doc/matrices.texi b/doc/matrices.texi index a95e8e29f0..8eb4842316 100644 --- a/doc/matrices.texi +++ b/doc/matrices.texi @@ -1549,6 +1549,8 @@ output in which @var{v} becomes 1 and other values become 0. @deffn {Matrix Function} DET (@var{M}) Returns the determinant of square matrix @var{M}. + +@t{DET(@{3, 7; 1, -4@}) @result{} -19} @end deffn @deffn {Matrix Function} EVAL (@var{M}) @@ -1558,6 +1560,8 @@ Returns a column vector containing the eigenvalues of symmetric matrix Use @code{CALL EIGEN} (@pxref{CALL EIGEN}) to compute eigenvalues and eigenvectors of a matrix. + +@t{EVAL(@{2, 0, 0; 0, 3, 4; 0, 4, 9@}) @result{} @{11; 2; 1@}} @end deffn @deffn {Matrix Function} GINV (@var{M}) @@ -1566,14 +1570,18 @@ Returns the @math{@var{k}@times{}@var{n}} matrix @var{A} that is the @var{M}, defined such that @math{@var{M}@times{}@var{A}@times{}@var{M}=@var{M}} and @math{@var{A}@times{}@var{M}@times{}@var{A}=@var{A}}. -@end deffn +@t{GINV(@{1, 2@}) @result{} @{.2; .4@}} (approximately) @* +@t{@{1:9@} * GINV(1:9) * @{1:9@} @result{} @{1:9@}} (approximately) +@end deffn @deffn {Matrix Function} GSCH (@var{M}) @var{M} must be a @math{@var{n}@times{}@var{m}} matrix, @math{@var{m} @geq{} @var{n}}, with rank @var{n}. Returns an @math{@var{n}@times{}@var{n}} orthonormal basis for @var{M}, obtained using the Gram-Schmidt process. + +@t{GSCH(@{3, 2; 1, 2@}) * SQRT(10) @result{} @{3, -1; 1, 3@}} (approximately) @end deffn @deffn {Matrix Function} INV (@var{M}) @@ -1582,6 +1590,8 @@ inverse of @math{@var{n}@times{}@var{n}} matrix @var{M}, defined such that @math{@var{M}@times{}@var{A} = @var{A}@times{}@var{M} = I}, where @var{I} is the identity matrix. @var{M} must not be singular, that is, @math{\det(@var{M}) @ne{} 0}. + +@t{INV(@{4, 7; 2, 6@}) @result{} @{.6, -.7; -.2, .4@}} (approximately) @end deffn @deffn {Matrix Function} KRONEKER (@var{Ma}, @var{Mb}) @@ -1593,12 +1603,28 @@ view @var{P} as the concatenation of multiple product of @var{Mb} by a different element of @var{Ma}. For example, when @code{A} is a @math{2@times{}2} matrix, @code{KRONEKER(A, B)} is equivalent to @code{@{A(1,1)*B, A(1,2)*B; A(2,1)*B, A(2,2)*B@}}. + +@format +@t{KRONEKER(@{1, 2; 3, 4@}, @{0, 5; 6, 7@}) @result{} + 0 5 0 10 + 6 7 12 14 + 0 15 0 20 + 18 21 24 28} +@end format @end deffn @deffn {Matrix Function} RANK (@var{M}) Returns the rank of matrix @var{M}, a integer scalar whose value is the dimension of the vector space spanned by its columns or, equivalently, by its rows. + +@format +@t{RANK(@{1, 0, 1; -2, -3, 1; 3, 3, 0@}) @result{} 2 +RANK(@{1, 1, 0, 2; -1, -1, 0, -2@}) @result{} 1 +RANK(@{1, -1; 1, -1; 0, 0; 2, -2@}) @result{} 1 +RANK(@{1, 2, 1; -2, -3, 1; 3, 5, 0@}) @result{} 2 +RANK(@{1, 0, 2; 2, 1, 0; 3, 2, 1@}) @result{} 3} +@end format @end deffn @deffn {Matrix Function} SOLVE (@var{Ma}, @var{Mb}) @@ -1607,6 +1633,22 @@ equivalently, by its rows. @math{@var{n}@times{}@var{k}} matrix. Returns an @math{@var{n}@times{}@var{k}} matrix @var{X} such that @math{@var{Ma} @times{} @var{X} = @var{Mb}}. + +All of the following examples show approximate results: + +@format +@t{SOLVE(@{2, 3; 4, 9@}, @{6, 2; 15, 5@}) @result{} + 1.50 .50 + 1.00 .33 +SOLVE(@{1, 3, -2; 3, 5, 6; 2, 4, 3@}, @{5; 7; 8@}) @result{} + -15.00 + 8.00 + 2.00 +SOLVE(@{2, 1, -1; -3, -1, 2; -2, 1, 2@}, @{8; -11; -3@}) @result{} + 2.00 + 3.00 + -1.00} +@end format @end deffn @deffn {Matrix Function} SVAL (@var{M}) @@ -1618,6 +1660,12 @@ singular values of @var{M} in descending order. Use @code{CALL SVD} (@pxref{CALL SVD}) to compute the full singular value decomposition of a matrix. + +@format +@t{SVAL(@{1, 1; 0, 0@}) @result{} @{1.41; .00@} +SVAL(@{1, 0, 1; 0, 1, 1; 0, 0, 0@}) @result{} @{1.73; 1.00; .00@} +SVAL(@{2, 4; 1, 3; 0, 0; 0, 0@}) @result{} @{5.46; .37@}} +@end format @end deffn @deffn {Matrix Function} SWEEP (@var{M}, @var{nk}) @@ -1632,7 +1680,7 @@ If @math{@var{M}_{kk} @ne{} 0}, then: @math{@var{A}_{kk} = 1/@var{M}_{kk}}, @math{@var{A}_{ik} = -@var{M}_{ik}/@var{M}_{kk} @r{for} i @ne{} k}, @math{@var{A}_{kj} = @var{M}_{kj}/@var{M}_{kk} @r{for} j @ne{} k, @r{and}} -@math{@var{A}_{ij} = @var{M}_{ij} - (@var{M}_{ik} * @var{M}_{kj}) / @var{M}_{kk} @r{for} i @ne{} k @r{and} j @ne{} k}. +@math{@var{A}_{ij} = @var{M}_{ij} - @var{M}_{ik}@var{M}_{kj}/@var{M}_{kk} @r{for} i @ne{} k @r{and} j @ne{} k}. @end display If @math{@var{M}_{kk} = 0}, then: @@ -1641,21 +1689,58 @@ If @math{@var{M}_{kk} = 0}, then: @math{@var{A}_{ik} = @var{A}_{ki} = 0 @r{and}} @math{@var{A}_{ij} = @var{M}_{ij}, @r{for} i @ne{} k @r{and} j @ne{} k}. @end display + +Given @t{M = @{0, 1, 2; 3, 4, 5; 6, 7, 8@}}, then (approximately): + +@format +@t{SWEEP(M, 1) @result{} + .00 .00 .00 + .00 4.00 5.00 + .00 7.00 8.00 +SWEEP(M, 2) @result{} + -.75 -.25 .75 + .75 .25 1.25 + .75 -1.75 -.75 +SWEEP(M, 3) @result{} + -1.50 -.75 -.25 + -.75 -.38 -.63 + .75 .88 .13} +@end format @end deffn @node Matrix IO -@subsubsection I/O +@subsubsection I/O Function + +This function works with files being used on the @code{READ} statement. @deffn {Matrix Function} EOF (@var{file}) @anchor{EOF Matrix Function} Given a file handle or file name @var{file}, returns an integer scalar -that indicates whether the last record in the file has been read. -Determining this requires attempting reading past the current record, +1 if the last line in the file has been read or 0 if more lines are +available. Determining this requires attempting to read another line, which means that @code{REREAD} on the next @code{READ} command -following @code{EOF} on the same file will be ineffective. +following @code{EOF} on the same file will be ineffective. @end deffn +The @code{EOF} function gives a matrix program the flexibility to read +a file with text data without knowing the length of the file in +advance. For example, the following program will read all the lines +of data in @file{data.txt}, each consisting of three numbers, as rows +in matrix @code{data}: + +@verbatim +MATRIX. +COMPUTE data={}. +LOOP IF NOT EOF('data.txt'). + READ row/FILE='data.txt'/FIELD=1 TO 1000/SIZE={1,3}. + COMPUTE data={data; row}. +END LOOP. +PRINT data. +END MATRIX. + +@end verbatim + @node Matrix COMPUTE Command @subsection The @code{COMPUTE} Command @@ -1697,6 +1782,40 @@ descending order to @var{n}-element column vector @var{eval}. Use the @code{EVAL} function (@pxref{EVAL}) to compute just the eigenvalues of a symmetric matrix. +For example, the following matrix language commands: +@example +CALL EIGEN(@{1, 0; 0, 1@}, evec, eval). +PRINT evec. +PRINT eval. + +CALL EIGEN(@{3, 2, 4; 2, 0, 2; 4, 2, 3@}, evec2, eval2). +PRINT evec2. +PRINT eval2. +@end example + +@noindent +yield this output: + +@example +evec + 1 0 + 0 1 + +eval + 1 + 1 + +evec2 + -.6666666667 .0000000000 .7453559925 + -.3333333333 -.8944271910 -.2981423970 + -.6666666667 .4472135955 -.5962847940 + +eval2 + 8.0000000000 + -1.0000000000 + -1.0000000000 +@end example + @item @t{CALL SVD(@var{M}, @var{U}, @var{S}, @var{V})} @anchor{CALL SVD} @@ -1709,6 +1828,22 @@ diagonal of @var{Q} contains the singular values of @var{M}. Use the @code{SVAL} function (@pxref{SVAL}) to compute just the singular values of a matrix. + +For example, the following matrix program: + +@example +CALL SVD(@{3, 2, 2; 2, 3, -2@}, u, s, v). +PRINT (u * s * T(v))/FORMAT F5.1. +@end example + +@noindent +yields this output: + +@example +(u * s * T(v)) + 3.0 2.0 2.0 + 2.0 3.0 -2.0 +@end example @end table The final procedure is implemented via @code{CALL} to allow it to @@ -1729,6 +1864,24 @@ extra elements of @var{V} are ignored. Use the @code{MDIAG} function (@pxref{MDIAG}) to construct a new matrix with a specified main diagonal. + +For example, this matrix program: + +@example +COMPUTE x=@{1, 2, 3; 4, 5, 6; 7, 8, 9@}. +CALL SETDIAG(x, 10). +PRINT x. +@end example + +@noindent +outputs the following: + +@example +x + 10 2 3 + 4 10 6 + 7 8 10 +@end example @end table @node Matrix PRINT Command @@ -1756,21 +1909,24 @@ the largest-magnitude element in the matrix to be displayed: @enumerate @item -If the matrix's elements are all integers, then, if possible, @pspp{} -chooses the narrowest @code{F} format with width 12 or less that fits -@var{m} plus a sign. +If @math{@var{m} < 10^{11}} and the matrix's elements are all +integers, @pspp{} chooses the narrowest @code{F} format that fits +@var{m} plus a sign. For example, if the matrix is @t{@{1:10@}}, then +@math{m = 10}, which fits in 3 columns with room for a sign, the +format is @code{F3.0}. @item Otherwise, if @math{@var{m} @geq{} 10^9} or @math{@var{m} @leq{} 10^{-4}}, @pspp{} scales all of the numbers in the matrix by -@math{10^x}, where @var{x} is the exponent used when @var{m} is -displayed in scientific notation, and displays the scaled value in -format @code{F13.10}. @pspp{} adds a note to the output to indicate -the scale factor. +@math{10^x}, where @var{x} is the exponent that would be used to +display @var{m} in scientific notation. For example, for +@math{@var{m} = 5.123@times{}10^{20}}, the scale factor is +@math{10^{20}}. @pspp{} displays the scaled values in format +@code{F13.10} and notes the scale factor in the output. @item -Otherwise, @pspp{} displays the value, without scaling, in format -@code{F13.10}. +Otherwise, @pspp{} displays the matrix values, without scaling, in +format @code{F13.10}. @end enumerate The optional @code{TITLE} subcommand specifies a title for the output @@ -1797,8 +1953,6 @@ the labels are truncated to 8 bytes. The @code{CLABELS} and @code{CNAMES} subcommands work for labeling columns as @code{RLABELS} and @code{RNAMES} do for labeling rows. -@subsubheading Text Output - When the @var{expression} is omitted, @code{PRINT} does not output a matrix. Instead, it outputs only the text specified on @code{TITLE}, if any, preceded by any space specified on the @code{SPACE} @@ -1806,6 +1960,8 @@ subcommand, if any. Any other subcommands are ignored, and the command acts as if @code{MDISPLAY} is set to @code{TEXT} regardless of its actual setting. + + @node Matrix DO IF Command @subsection The @code{DO IF} Command diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c index 636c242ee1..c59b3fd2db 100644 --- a/src/language/stats/matrix.c +++ b/src/language/stats/matrix.c @@ -2381,7 +2381,7 @@ matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_, return NULL; } - long int n = (end > start && by > 0 ? (end - start + by) / by + 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); @@ -4811,7 +4811,8 @@ matrix_cmd_execute_read (struct read_command *read) if (!is_vector (m)) { - msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector")); + msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, " + "not a %zu×%zu matrix."), m->size1, m->size2); gsl_matrix_free (m); free (iv0.indexes); free (iv1.indexes); @@ -4832,7 +4833,9 @@ matrix_cmd_execute_read (struct read_command *read) } else { - msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector")); + msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, " + "not a %zu×%zu matrix."), + m->size1, m->size2), gsl_matrix_free (m); free (iv0.indexes); free (iv1.indexes); diff --git a/tests/language/stats/matrix.at b/tests/language/stats/matrix.at index cc111c4952..71e8d280dc 100644 --- a/tests/language/stats/matrix.at +++ b/tests/language/stats/matrix.at @@ -1122,6 +1122,9 @@ PRINT {-1:-3:-1}. PRINT {-1:-10:-2}. PRINT {-1:-11:-2}. +PRINT {1:1}. +PRINT {1:1:-1}. + PRINT {1:3:0}. PRINT {-1:-3:0}. @@ -1152,6 +1155,12 @@ AT_CHECK([pspp matrix.sps], [1], [dnl {-1:-11:-2} -1 -3 -5 -7 -9 -11 +{1:1} + 1 + +{1:1:-1} + 1 + matrix.sps:12: error: MATRIX: The increment operand to : must be nonzero. matrix.sps:13: error: MATRIX: The increment operand to : must be nonzero. @@ -1857,6 +1866,11 @@ 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). +COMPUTE s2 = {0, 1, 2; 3, 4, 5; 6, 7, 8}. +PRINT SWEEP(s2, 1). +PRINT SWEEP(s2, 2). +PRINT SWEEP(s2, 3). + PRINT TRACE(s0). PRINT T(s0). @@ -2031,9 +2045,9 @@ 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. +CALL EIGEN({3, 2, 4; 2, 0, 2; 4, 2, 3}, evec2, eval2). +PRINT evec2. +PRINT eval2. END MATRIX. ]) AT_CHECK([pspp matrix.sps], [0], [dnl @@ -2045,12 +2059,12 @@ eval 1 1 -evec +evec2 -.6666666667 .0000000000 .7453559925 -.3333333333 -.8944271910 -.2981423970 -.6666666667 .4472135955 -.5962847940 -eval +eval2 8.0000000000 -1.0000000000 -1.0000000000 -- 2.30.2