diff --git a/CITATION.cff b/CITATION.cff index f52783c..a08d501 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -7,16 +7,16 @@ authors: title: matrico type: software abstract: "A flonum matrix module for CHICKEN Scheme." -version: 0.4 -date-released: 2023-06-01 +version: 0.5 +date-released: 2023-06-06 commit: license: zlib-Acknowledgement repository: https://github.com/gramian/matrico url: https://numerical-schemer.xyz keywords: + - "Lisp" - "Scheme" - "Chicken Scheme" - - "Matrix Calculations" - "Matrix Computations" - "Linear Algebra" - + - "Numerics" diff --git a/Makefile b/Makefile index 57164a2..8e5635d 100644 --- a/Makefile +++ b/Makefile @@ -8,6 +8,11 @@ CHICKEN_INSTALL = chicken-install CHICKEN_PROFILE = chicken-profile TEST_NEW_EGG = test-new-egg # homebrew: /opt/homebrew/Cellar/chicken/5.3.0_1/bin/test-new-egg +# DEV-TEMP +TAR = gtar +SED = gsed +TEST_NEW_EGG = /opt/homebrew/Cellar/chicken/5.3.0_1/bin/test-new-egg + CLARG = LEVEL = -O3 DEBUG = -d0 @@ -47,7 +52,7 @@ test_install: CHICKEN_REPOSITORY_PATH="/tmp/matrico:`$(CHICKEN_INSTALL) -repository`" $(CSI) -e "(import matrico) (matrico 'citation)" iprofile: - make matmul LEVEL='-O3' FLAGS='-profile' DIM=100 + make linpack LEVEL='-O3' FLAGS='-profile' DIM=1000 ls -1 PROFILE.* | sort -n | head -n1 | $(CHICKEN_PROFILE) sprofile: @@ -100,4 +105,4 @@ clean: rm -f test.csv test.mx linpack.txt matmul.txt \ matrico.so matrico.import.scm matrico.import.so matrico.static.o \ matrico.build.sh matrico.install.sh matrico.link matrico.static.so \ - matrico-$(VERSION).tar.gz PROFILE.* + matrico-$(VERSION).tar.gz PROFILE.* heat.csv flame.csv diff --git a/README.md b/README.md index 45e181c..37cd447 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -![matrico 0.4](res/matrico-logo.svg) matrico +![matrico 0.5](res/matrico-logo.svg) matrico ======================================== * **Project**: matrico ([Esperanto for "matrix"](https://translate.google.com/?sl=eo&tl=en&text=matrico&op=translate)) @@ -9,7 +9,7 @@ * **License**: [zlib-acknowledgement](https://spdx.org/licenses/zlib-acknowledgement.html) -* **Version**: 0.4 (2023-06-01) +* **Version**: 0.5 (2023-06-06) * **Depends**: [CHICKEN Scheme](http://call-cc.org) (>= 5.1) @@ -37,6 +37,12 @@ `matrico` is a _Scheme_ module for numerical matrix computations encapsulated in a _CHICKEN Scheme_ egg. +### Clone and Try `matrico` Code + +```shell +./matrico.sh +``` + ### Install and Test `matrico` Egg ```shell @@ -49,26 +55,18 @@ CHICKEN_INSTALL_REPOSITORY="/my/egg/directory/" chicken-install -test matrico CHICKEN_REPOSITORY_PATH="`chicken-install -repository`:/my/egg/directory/" csi ``` -### Clone and Try `matrico` Code - -```scheme -(load "matrico.scm") +### Run Demo Codes -(import matrico) -``` - -### Run Demo Code - -```scheme -(load "RUNME.scm") +```shell +csi RUNME.scm ``` -```scheme -(load "demos/heat.scm") +```shell +csi demos/heat.scm ``` -```scheme -(load "demos/flame.scm") +```shell +csi demos/flame.scm ``` ### Minimal Explanation @@ -250,8 +248,6 @@ CHICKEN_REPOSITORY_PATH="`chicken-install -repository`:/my/egg/directory/" csi * `(mx-where pred x y)` returns **matrix** of entries of **matrix**es `x` or `y` based on predicate **procedure** `pred`. -* `(mx*+ a x y)` returns **matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`, aka axpy. - #### Matrix Mappers ##### Elementary Functions @@ -422,6 +418,8 @@ CHICKEN_REPOSITORY_PATH="`chicken-install -repository`:/my/egg/directory/" csi * `(mx-transpose mat)` returns **matrix** of entries of **matrix** `mat` with swapped row and column indices. +* `(mx-axpy a x y)` returns **matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`, aka "a times x plus y". + * `(mx-sympart mat)` returns **matrix** being symmetric part of square **matrix** `mat`. * `(mx-skewpart mat)` returns **matrix** being skey-symmetric part of square **matrix** `mat`, aka anti-symmetric part. @@ -611,6 +609,8 @@ Defines the matrix type (record) as column-major list-of-columns and provides ge * `(matrix-transpose mat)` returns **matrix** of entries of **matrix** `mat` with swapped row and column indices. +* `(matrix-axpy a x y)` returns **matrix** resulting from scaling **matrix** `x` by **any** `a` and add **matrix** `y`. + * `(matrix-scalar xt y)` returns **any** resulting from the scalar product of column-**matrix**es `xt` and `y`. * `(matrix-dot* xt y)` returns **matrix** resulting from matrix multiplication of transposed of **matrix** `xt` and **matrix** `y`. @@ -655,6 +655,8 @@ Provides homogeneous vector transformations analogous to vectors. * `(f64vector-foreach-index fun . vecs)` returns **void**, applies **procedure** `fun` to index and all corresponding **f64vector**(s) `vecs` elements. +* `(f64vector-axpy a x y)` returns **f64vector** resulting from applying fused-multiply-add to the **flonum** `a` and to all **f64vector**s `x`, `y` elements. + * `(f64vector-fold fun ini . vecs)` returns **any** resulting from applying **procedure** `fun` to `ini` initialized accumulator and sequentially to all **f64vector**(s) `vecs` elements from left to right. * `(f64vector-fold* fun ini . vecs)` returns **any** resulting from applying **procedure** `fun` to `ini` initialized accumulator and sequentially to all **f64vector**(s) `vecs` elements from right to left. @@ -678,7 +680,7 @@ Provides additional mathematical functions, extending CHICKEN's `flonum` module. * `(fprec x)` returns **flonum** reciprocal of **flonum** `x`. -* `(fp*+ z x y)` returns **flonum** sum with product: `x * y + z` of **flonum**s `x`, `y`, `z`. +* `(fp*+ x y z)` returns **flonum** sum with product: `x * y + z` of **flonum**s `x`, `y`, `z`. * `(fptau)` returns **flonum** circle constant Tau via fraction. @@ -844,8 +846,8 @@ make mips * `SYS:` MacOS Monterey (12.6) * `SCM:` CHICKEN Scheme (5.3) -* MATMUL: `1230` Megaflops -* LINPACK: `185` Megaflops +* MATMUL: `267` Megaflops +* LINPACK: `319` Megaflops * BOGOMIPS: `275` Mips ## Development @@ -858,13 +860,21 @@ make mips * [`mx`] Add rank-revealing QR and pseudo-inverse via QR -* [`mx`] Add Eigenvalue decomposition and singular value decomposition +* [`mx`] Add Eigenvalue decomposition and singular value decomposition via QR -* [`mx`] Add [`asciichart`](https://github.com/kroitor/asciichart) interface or functionality +* [`mx`] Add [`UnicodePlot`](https://github.com/JuliaPlots/UnicodePlots.jl)-like lineplot functionality ### Changelog -0.4 (2023-06-01) +0.5 (2023-06-06) + + * **ADDED** `f64vector-axpy` + * **ADDED** `matrix-axpy` + * **CHANGED** `mx*+` to `mx-axpy` + * **IMPROVED** `mx-qr` + * **IMPROVED** `mx-solver` + +
0.4 (2023-06-01) * **ADDED** `mx-angle` * **ADDED** `mx-var` @@ -875,6 +885,8 @@ make mips * **REMOVED** `mx-logb` * ... and many minor updates and fixes. +
+
0.3 (2022-09-16) * **ADDED** `matrico` function @@ -937,12 +949,4 @@ make mips * `matrico` can be build with `-O5` if `matrico.scm` is included into the source. -### Ideas - -* Alias for `.bashrc`: -``` -alias matrico='CHICKEN_REPOSITORY_PATH="`chicken-install -repository`:/my/egg/directory/" csi -R matrico' -``` - ## [`matrico`](https://git.io/matrico) — a (:chicken: λ) :egg: for numerical schemers! - diff --git a/RUNME.scm b/RUNME.scm index ade54eb..49489a4 100644 --- a/RUNME.scm +++ b/RUNME.scm @@ -1,7 +1,7 @@ ;;;; RUNME.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: demo code diff --git a/demos/flame.scm b/demos/flame.scm index 333e365..9c2cdb9 100644 --- a/demos/flame.scm +++ b/demos/flame.scm @@ -1,7 +1,7 @@ ;;;; demo-flame.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: flame demo code @@ -19,6 +19,7 @@ (define 2/h (/ 0.5 h)) (define A^T (mx-tridiag N (* 1.0 2/h) (* -4.0 2/h) (* 3.0 2/h))) +(mx-set! A^T (sub1 N) N (* 4.0 2/h)) (define B (mx-set (mx N 1 0.0) 1 1 (* 1.0 2/h))) (define C^T (mx-set (mx N 1 0.0) N 1 1.0)) diff --git a/demos/heat.scm b/demos/heat.scm index 017544f..6167bfe 100644 --- a/demos/heat.scm +++ b/demos/heat.scm @@ -1,7 +1,7 @@ ;;;; demo-heat.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: cooling demo code diff --git a/doc/matrico.wiki b/doc/matrico.wiki index e2db1c6..01e73eb 100644 --- a/doc/matrico.wiki +++ b/doc/matrico.wiki @@ -165,9 +165,6 @@ returns ''matrix'' of entry-wise exponentiation of ''matrix''es {{x}} to the {{y (mx-where pred x y) returns ''matrix'' of entries of ''matrix''es {{x}} or {{y}} based on predicate ''procedure'' {{pred}}. -(mx*+ a x y) -returns ''matrix'' of entry-wise generalized addition of ''flonum'' {{a}} times ''matrix'' {{x}} plus ''matrix'' {{y}}. - == Matrix Mappers === Entry-Wise Elementary Functions @@ -400,6 +397,9 @@ returns ''matrix'' of vertically concatenated columns of ''matrix'' {{mat}}, aka (mx-transpose mat) returns ''matrix'' of entries of ''matrix'' {{mat}} with swapped row and column indices. +(mx-axpy a x y) +returns ''matrix'' of entry-wise generalized addition of ''flonum'' {{a}} times ''matrix'' {{x}} plus ''matrix'' {{y}}. + (mx-sympart mat) returns ''matrix'' being symmetric part of square ''matrix'' {{mat}}. @@ -536,8 +536,10 @@ Copyright (c) 2022 Christian Himpe. [[https://spdx.org/licenses/zlib-acknowledge ; 0.1 : Initial Release -; 0.2 : Major Update ([[https://github.com/gramian/matrico#changelog|details]]) +; 0.2 : Major Update + +; 0.3 : Major Update -; 0.3 : Major Update ([[https://github.com/gramian/matrico#changelog|details]]) +; 0.4 : Minor Update -; 0.4 : Minor Update ([[https://github.com/gramian/matrico#changelog|details]]) +; 0.5 : Minor Update ([[https://github.com/gramian/matrico#changelog|details]]) \ No newline at end of file diff --git a/matrico.egg b/matrico.egg index 5c1d311..2403a28 100644 --- a/matrico.egg +++ b/matrico.egg @@ -1,6 +1,6 @@ ;;;; matrico.egg -((version "0.4") +((version "0.5") (synopsis "A flonum matrix module for CHICKEN Scheme.") (author "Christian Himpe") (category math) diff --git a/matrico.release-info b/matrico.release-info index a742b78..c313695 100644 --- a/matrico.release-info +++ b/matrico.release-info @@ -14,3 +14,5 @@ (release "0.4rc1") (release "0.4rel") + +(release "0.5rc1") diff --git a/matrico.scm b/matrico.scm index 18cce63..82477f2 100644 --- a/matrico.scm +++ b/matrico.scm @@ -1,7 +1,7 @@ ;;;; matrico.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: A CHICKEN Scheme flonum matrix module. @@ -21,7 +21,7 @@ mx-samecols? mx-samerows? mx-samedims? mx-any? mx-all? mx=? mx-ref11 mx-ref mx-set mx-set! mx-col mx-row mx-diag mx-submatrix - mx+ mx* mx- mx/ mx*2 mx^2 mx^ mx-where mx*+ + mx+ mx* mx- mx/ mx*2 mx^2 mx^ mx-where mx-round mx-floor mx-ceil mx-abs mx-sign mx-delta mx-heaviside mx-sin mx-cos mx-tan @@ -43,6 +43,7 @@ mx-rownorm mx-colnorm mx-norm mx-horcat mx-vercat mx-vec mx-transpose + mx-axpy mx-sympart mx-skewpart mx-diagonal mx-qr mx-solver mx-solve mx-orth mx-absdet mx-logdet @@ -67,7 +68,7 @@ [else]) ;;@assigns: matrico version number as pair. -(define-constant version '(0 . 4)) +(define-constant version '(0 . 5)) ;;@returns (define (matrico . s) diff --git a/src/dense.scm b/src/dense.scm index 3493e54..a1a812c 100644 --- a/src/dense.scm +++ b/src/dense.scm @@ -1,7 +1,7 @@ ;;;; dense.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: dense column function aliases @@ -28,6 +28,7 @@ (f64vector-map-index column-map-index) (f64vector-foreach column-foreach) (f64vector-foreach-index column-foreach-index) + (f64vector-axpy column-axpy) (f64vector-fold column-fold) (f64vector-fold* column-fold*) (f64vector-dot column-dot) diff --git a/src/f64vector.scm b/src/f64vector.scm index ed6c86f..e3d3fab 100644 --- a/src/f64vector.scm +++ b/src/f64vector.scm @@ -1,7 +1,7 @@ ;;;; f64vector.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: homogeneous flonum vector library @@ -15,7 +15,9 @@ f64vector-any? f64vector-all? f64vector-map f64vector-map-index f64vector-foreach f64vector-foreach-index - f64vector-fold f64vector-fold* f64vector-dot) + f64vector-axpy + f64vector-fold f64vector-fold* + f64vector-dot) (import scheme (chicken base) (chicken module) srfi-4 utils fpmath) @@ -122,6 +124,16 @@ (apply fun idx (map (cut f64vector-ref <> idx) vecs)) (rho (fx+1 idx))))))) +;;@returns: **f64vector** resulting from applying fused-multiply-add to the **flonum** `a` and to all **f64vector**s `x`, `y` elements. +(define (f64vector-axpy a x y) ;@marker: Hot + (define dim (f64vector-length x)) + (define ret (make-f64vector dim)) + (let rho [(idx 0)] + (if (fx= idx dim) ret + (begin + (f64vector-set! ret idx (fp*+ a (f64vector-ref x idx) (f64vector-ref y idx))) + (rho (fx+1 idx)))))) + ;;; Vector Reducers ############################################################ ;;@returns: **any** resulting from applying **procedure** `fun` to `ini` initialized accumulator and sequentially to all **f64vector**(s) `vecs` elements from left to right. @@ -141,23 +153,11 @@ ;;@returns: **flonum** resulting from applying fused-multiply-add to zero initialized accumulator and sequentially to all **f64vector**s `x`, `y` elements from left to right. (define (f64vector-dot x y) ;@marker: Hot - (let [(dot (cond-expand - [compiling - (begin - (import (chicken foreign)) - (foreign-lambda* double ((f64vector x) (f64vector y) (size_t dim)) -" double res = 0.0; - for(ptrdiff_t i = 0; i < dim; ++i) - //res = fma(x[i], y[i], res); - res += x[i] * y[i]; - C_return(res);"))] - [else - (lambda (x y dim) - (let rho [(ret 0.0) - (idx 0)] - (if (fx= idx dim) ret - (rho (fp*+ (f64vector-ref x idx) (f64vector-ref y idx) ret) (fx+1 idx)))))]))] - (dot x y (f64vector-length x)))) + (define dim (f64vector-length x)) + (let rho [(idx 0) + (ret 0.0)] + (if (fx= idx dim) ret + (rho (fx+1 idx) (fp*+ (f64vector-ref x idx) (f64vector-ref y idx) ret))))) );end module diff --git a/src/fpmath.scm b/src/fpmath.scm index e77a489..afc819b 100644 --- a/src/fpmath.scm +++ b/src/fpmath.scm @@ -1,7 +1,7 @@ ;;;; fpmath.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: floating-point add-on module @@ -116,30 +116,48 @@ ;;; Hyperbolic Functions ####################################################### ;;@returns: **flonum** hyperbolic sine of **flonum** `x`: `sinh(x) = 0.5 * (exp(x) - exp(-x))`. -(define (fpsinh x) - (fp* 0.5 (fp- (fpexp x) (fpexp (fpneg x))))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fpsinh x) + (fp* 0.5 (fp- (fpexp x) (fpexp (fpneg x)))))] + [else]) ;;@returns: **flonum** hyperbolic cosine of **flonum** `x`: `cosh(x) = 0.5 * (exp(x) + exp(-x))`. -(define (fpcosh x) - (fp* 0.5 (fp+ (fpexp x) (fpexp (fpneg x))))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fpcosh x) + (fp* 0.5 (fp+ (fpexp x) (fpexp (fpneg x)))))] + [else]) ;;@returns: **flonum** hyperbolic tangent of **flonum** `x`: `tanh(x) = 1 - 2 / (exp(2 * x) + 1)`. -(define (fptanh x) - (fp- 1.0 (fp/ 2.0 (fp+ (fpexp (fp*2 x)) 1.0)))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fptanh x) + (fp- 1.0 (fp/ 2.0 (fp+ (fpexp (fp*2 x)) 1.0))))] + [else]) ;;; Inverse Hyperbolic Functions ############################################### ;;@returns: **flonum** area hyperbolic sine of **flonum** `x`: `asinh(x) = log(x + sqrt(x^2 + 1))`. -(define (fpasinh x) - (fp* (fpsign x) (fplog (fp+ (fpabs x) (fpsqrt (fp+ (fp^2 x) 1.0)))))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fpasinh x) + (fp* (fpsign x) (fplog (fp+ (fpabs x) (fpsqrt (fp+ (fp^2 x) 1.0))))))] + [else]) ;;@returns: **flonum** area hyperbolic cosine of **flonum** `x`: `acosh(x) = log(x + sqrt(x^2 - 1))`. -(define (fpacosh x) - (fplog (fp+ x (fpsqrt (fp- (fp^2 x) 1.0))))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fpacosh x) + (fplog (fp+ x (fpsqrt (fp- (fp^2 x) 1.0)))))] + [else]) ;;@returns: **flonum** area hyperbolic tangent of **flonum** `x`: `atanh(x) = 0.5 * log((1 + x) / (1 - x))`. -(define (fpatanh x) - (fp* 0.5 (fplog (fp/ (fp+ 1.0 x) (fp- 1.0 x))))) +(cond-expand + [(or chicken-5.0 chicken-5.1 chicken-5.2 chicken-5.3) + (define (fpatanh x) + (fp* 0.5 (fplog (fp/ (fp+ 1.0 x) (fp- 1.0 x)))))] + [else]) ;;; Haversed Trigonometric Functions ########################################### diff --git a/src/matrix.scm b/src/matrix.scm index 2d65660..9f599bf 100644 --- a/src/matrix.scm +++ b/src/matrix.scm @@ -1,7 +1,7 @@ ;;;; matrix.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: matrix type back-end via list-of-vectors, see @1, @2, @3, @4, @5, @6, @7, @8. @@ -15,7 +15,7 @@ column-ref subcolumn column-set! column-any? column-all? - column-map column-map-index column-foreach column-foreach-index + column-map column-map-index column-foreach column-foreach-index column-axpy column-fold column-fold* column-dot list->column))) @@ -30,7 +30,7 @@ matrix-any? matrix-all? matrix-colfold matrix-rowfold matrix-allfold matrix-map matrix-broadcast - matrix-vec matrix-transpose matrix-scalar matrix-dot* matrix-explode matrix-implode + matrix-vec matrix-transpose matrix-axpy matrix-scalar matrix-dot* matrix-explode matrix-implode matrix-print matrix-export matrix-save matrix-load) (import scheme (chicken base) (chicken module) (only (chicken memory representation) object-copy) utils column) @@ -318,6 +318,12 @@ nil (matrix-data mat)))) +;;@returns: **matrix** resulting from scaling **matrix** `x` by **any** `a` and add **matrix** `y`. +(define (matrix-axpy a x y) + (matrix-map** (lambda (xcol ycol) + (column-axpy a xcol ycol)) + x y)) + ;;@returns: **any** resulting from the scalar product of column-**matrix**es `xt` and `y`. (define (matrix-scalar xt y) (column-dot (head (matrix-data xt)) (head (matrix-data y)))) diff --git a/src/mx.scm b/src/mx.scm index 0ddd451..e2e3bfe 100644 --- a/src/mx.scm +++ b/src/mx.scm @@ -1,7 +1,7 @@ ;;;; mx.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: matrix type front-end @@ -361,12 +361,6 @@ (must-be (matrix-or-flonum? x) (matrix-or-flonum? y)) (matrix-broadcast (lambda (x y) (if (pred x y) x y)) (ensure-mx x) (ensure-mx y))) -;;@returns: **matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`, aka "axpy" (a times x plus y), see @2. -(define* (mx*+ a x y) - (returns "**matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`.") - (must-be (scalar-or-flonum? a) (matrix-or-flonum? x) (matrix-or-flonum? y)) - (matrix-broadcast (cute fp*+ (ensure-fp a) <> <>) (ensure-mx x) (ensure-mx y))) - ;;; Matrix Mappers ############################################################# ;; Rounding Functions @@ -786,6 +780,12 @@ (must-be (matrix? mat)) (matrix-transpose mat)) +;;@returns: **matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`, aka "a times x plus y", see @2. +(define* (mx-axpy a x y) + (returns "**matrix** of entry-wise generalized addition of **flonum** `a` times **matrix** `x` plus **matrix** `y`.") + (must-be (scalar-or-flonum? a) (matrix-samedims? x y)) + (matrix-axpy (ensure-fp a) (ensure-mx x) (ensure-mx y))) + ;;@returns: **matrix** being symmetric part of square **matrix** `mat`. (define* (mx-sympart mat) (returns "**matrix** being symmetric part of square **matrix** `mat`.") @@ -811,47 +811,41 @@ ; TODO: sort by Rkk value and keep largest ;;@returns: **pair** of an orthogonal **matrix** (Q) and upper-right triangular **matrix** (R) factoring full column rank **matrix** `mat`, via column-wise modified Gram-Schmidt, see @5, @6. -(define* (mx-qr mat) +(define* (mx-qr mat . rev) (returns "**pair** of an orthogonal **matrix** (Q) and upper-right triangular **matrix** (R) factoring full column rank **matrix** `mat`.") (must-be (matrix? mat)) - (define cols (matrix-cols mat)) - (define rows (matrix-rows mat)) + (define qrows (matrix-rows mat)) + (define rrows (fxmin (matrix-cols mat) qrows)) (let rho [(Ak (matrix-explode mat)) (Q nil) (R nil)] - (if (empty? Ak) (cons (matrix-implode Q) (matrix-implode (reverse R))) - (let [(Rk (mx (fxmin cols rows) 1 0.0))] + (if (empty? Ak) (cons (matrix-implode Q) (matrix-implode (if (optional rev #t) (reverse R) R))) + (let [(Rk (make-matrix* rrows 1 0.0))] (rho (tail Ak) (let cmgs [(i 0) (Qk (head Ak)) (Qi Q)] - (cond [(fx>= i rows) Q] - [(empty? Qi) (let [(Rkk (mx-norm Qk 'fro))] - (matrix-set! Rk i 0 Rkk) - (append* Q (mx/ Qk Rkk)))] - [else (let [(Rik (matrix-scalar (head Qi) Qk))] - (matrix-set! Rk i 0 Rik) - (cmgs (fx+1 i) (mx*+ (fpneg Rik) (head Qi) Qk) (tail Qi)))])) + (cond [(fx>= i qrows) Q] + [(empty? Qi) (append* Q (mx/ Qk (matrix-set! Rk i 0 (mx-norm Qk 'fro))))] + [else (cmgs (fx+1 i) (mx-axpy (fpneg (matrix-set! Rk i 0 (matrix-scalar (head Qi) Qk))) (head Qi) Qk) (tail Qi))])) (cons Rk R)))))) ;;@returns: **function** returning column-**matrix** solving the linear (least-squares) problem of **matrix** `mat`, given a column-**matrix** `vec` via QR decomposition. (define* (mx-solver mat) (returns "**function** returning column-**matrix** solving the linear (least-squares) problem of **matrix** `mat`, given a column-**matrix** `vec` via QR decomposition.") (must-be (matrix? mat)) - (let [(qr (mx-qr mat))] + (let* [(qr (mx-qr mat #f))] (lambda (vec) (let [(x (mx-dot* (head qr) vec))] - (let rho [(rcol (matrix-explode (tail qr)))] - (if (empty? rcol) (fx-1 (matrix-cols (tail qr))) - (let* [(cur (rho (tail rcol))) - (val (fpneg (matrix-set! x cur 0 (fp/ (matrix-ref*0 x cur) (matrix-ref*0 (head rcol) cur))))) - (nxt (fx-1 cur))] - (let sho [(idx nxt)] - (if (fx<0? idx) nxt + (let rho [(rcol (matrix-explode (tail qr))) + (xrow (fx-1 (matrix-cols (tail qr))))] + (if (empty? rcol) x + (let [(val (fpneg (matrix-set! x xrow 0 (fp/ (matrix-ref*0 x xrow) (matrix-ref*0 (head rcol) xrow)))))] + (let sho [(idx (fx-1 xrow))] + (if (fx<0? idx) (rho (tail rcol) (fx-1 xrow)) (begin (matrix-set! x idx 0 (fp*+ val (matrix-ref*0 (head rcol) idx) (matrix-ref*0 x idx))) - (sho (fx-1 idx)))))))) - x)))) + (sho (fx-1 idx)))))))))))) ;;@returns: column-**matrix** solving the linear (least-squares) problem of **matrix** `mat` and column-**matrix** `vec` via QR decomposition. (define* (mx-solve mat vec) @@ -1010,7 +1004,7 @@ (define* (mx-trapz mat) (returns "column-**matrix** trapezoid approximate integral of **matrix** `mat` being columns data-points of rows-dimensional function.") (must-be (matrix? mat)) - (mx*+ -0.5 (mx-col mat -1) (mx*+ -0.5 (mx-col mat 1) (mx-rowsum mat)))) + (mx-axpy -0.5 (mx-col mat -1) (mx-axpy -0.5 (mx-col mat 1) (mx-rowsum mat)))) ; TODO: type checking ;;@returns: states-times-steps **matrix** trajectory solving an ordinary differential equation, by a 2nd order hyperbolic Runge-Kutta method of **fixnum** `num` stages, with vector field **procedure** `fun`, initial state column-**matrix** `x0`, time step **flonum** `dt`, and time horizon **flonum** `tf`, see @7. @@ -1033,7 +1027,7 @@ (let rho [(cur coeff) (ret xk)] (if (empty? cur) ret - (rho (tail cur) (mx*+ (fp* dt (head cur)) (vf (fp*+ dt (head cur) tk) ret) xk))))))] + (rho (tail cur) (mx-axpy (fp* dt (head cur)) (vf (fp*+ dt (head cur) tk) ret) xk))))))] (time-stepper hyp sys tim x0))) ; TODO: type checking @@ -1047,8 +1041,8 @@ (let rho [(cur (fx-1 num)) (c tk) (ret xk)] - (if (fx=0? cur) (mx/ (mx*+ dt (vf c ret) (mx*+ s-1 ret xk)) (fp num)) - (rho (fx-1 cur) (fp+ tk dt/s-1) (mx*+ dt/s-1 (vf c ret) ret)))))))] + (if (fx=0? cur) (mx/ (mx-axpy dt (vf c ret) (mx-axpy s-1 ret xk)) (fp num)) + (rho (fx-1 cur) (fp+ tk dt/s-1) (mx-axpy dt/s-1 (vf c ret) ret)))))))] (time-stepper ssp sys tim x0))) ;;; Matrix Utilities ########################################################### diff --git a/src/utils.scm b/src/utils.scm index 6142049..904e369 100644 --- a/src/utils.scm +++ b/src/utils.scm @@ -1,7 +1,7 @@ ;;;; utils.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: helper utilities module diff --git a/tests/check.scm b/tests/check.scm index 45018c2..b4b1bad 100644 --- a/tests/check.scm +++ b/tests/check.scm @@ -1,7 +1,7 @@ ;;;; check.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: testing helpers @@ -21,4 +21,3 @@ (let [(pass? (equal=? (handle-exceptions exn 'xfail (eval `,(cons exe (caar tst)))) (cdar tst)))] ; `eval`, not `apply` because macros (print* " " (if pass? "\x1b[32mOK" "\x1b[31mFAIL") "\x1b[0m") (rho (cdr tst) (ok? (and (ok?) pass?))))))) - diff --git a/tests/run.scm b/tests/run.scm index a4ec2c5..ee5f6d5 100644 --- a/tests/run.scm +++ b/tests/run.scm @@ -1,7 +1,7 @@ ;;;; run.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: test runner diff --git a/tests/test-f64vector.scm b/tests/test-f64vector.scm index 7bdba11..7a70ef2 100644 --- a/tests/test-f64vector.scm +++ b/tests/test-f64vector.scm @@ -1,7 +1,7 @@ ;;;; test-f64vector.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: f64vector module unit tests @@ -24,7 +24,7 @@ ( (#f64(1.0 2.0)) . #f64(1.0 2.0)) ( (#f64()) . #f64()))) -;; f64vector-or? +;; f64vector-any? (check 'f64vector-any? '(( (zero? #f64(0.0 1.0 2.0)) . #t) ( (zero? #f64(1.0 2.0 0.0)) . #t) ( (zero? #f64(0.0 0.0 0.0)) . #t) @@ -32,7 +32,7 @@ ( (zero? #f64(1.0 2.0 3.0)) . #f) ( (zero? #f64()) . #f))) -;; f64vector-and? +;; f64vector-all? (check 'f64vector-all? '(( (zero? #f64(0.0 1.0 2.0)) . #f) ( (zero? #f64(1.0 2.0 0.0)) . #f) ( (zero? #f64(0.0 0.0 0.0)) . #t) @@ -55,6 +55,12 @@ ; TODO f64vector-foreach-index +;; f64vector-axpy +(check 'f64vector-axpy '(((2.0 #f64(1.0 2.0 3.0) #f64(1.0 2.0 3.0)) . #f64(3.0 6.0 9.0)) + ((0.0 #f64(1.0 2.0 3.0) #f64(1.0 2.0 3.0)) . #f64(1.0 2.0 3.0)) + ((2.0 #f64(0.0 0.0 0.0) #f64(1.0 2.0 3.0)) . #f64(1.0 2.0 3.0)) + ((2.0 #f64(1.0 2.0 3.0) #f64(0.0 0.0 0.0)) . #f64(2.0 4.0 6.0)))) + ;; f64vector-fold (check 'f64vector-fold '(((+ 0.0 #f64(1.0 2.0 3.0)) . 6.0) ( (+ 0.0 #f64(1.0)) . 1.0) @@ -65,7 +71,7 @@ ( (+ 0.0 #f64(1.0)) . 1.0) ( (+ 0.0 #f64()) . 0.0))) +;; f64vector-dot (check 'f64vector-dot '(((#f64(2.0) #f64(3.0)) . 6.0) ((#f64(1.0 2.0 3.0) #f64(2.0 3.0 4.0)) . 20.0) ((#f64() #f64()) . 0.0))) - diff --git a/tests/test-fpmath.scm b/tests/test-fpmath.scm index 0c1f545..4e26675 100644 --- a/tests/test-fpmath.scm +++ b/tests/test-fpmath.scm @@ -1,7 +1,7 @@ ;;;; test-fpmath.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: fpmath module unit tests @@ -246,4 +246,3 @@ ((+nan.0) . " NaN ") ((+inf.0) . " +∞ ") ((-inf.0) . " -∞ "))) - diff --git a/tests/test-matrico.scm b/tests/test-matrico.scm index d866ee6..a5a9df8 100644 --- a/tests/test-matrico.scm +++ b/tests/test-matrico.scm @@ -1,7 +1,7 @@ ;;;; test-matrico.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: matrico (compounding matrix/dense/mx) module unit tests @@ -547,32 +547,6 @@ (3.0 1.0 3.0) (4.0 1.0 4.0)))))) -;; mx+* -(check 'mx*+ `(( (0.0 m33 m33) . ,m33) - ((-1.0 m33 m33) . ,(mx 3 3 0.0)) - ( (0.0 m31 m31) . ,m31) - ((-1.0 m31 m31) . ,(mx 3 1 0.0)) - ( (0.0 m13 m13) . ,m13) - ((-1.0 m13 m13) . ,(mx 1 3 0.0)) - ( (0.0 m11 m11) . ,m11) - ((-1.0 m11 m11) . ,(mx 1 1 0.0)) - ( (1.0 1.0 m33) . ,(mx% '((9.0 2.0 7.0) - (4.0 6.0 8.0) - (5.0 10.0 3.0)))) - ( (1.0 1.0 m31) . ,(mx% '((9.0) - (4.0) - (5.0)))) - ( (1.0 1.0 m13) . ,(mx% '((9.0 2.0 7.0)))) - ( (1.0 1.0 m11) . ,(mx% '((9.0)))) - ( (1.0 m33 1.0) . ,(mx% '((9.0 2.0 7.0) - (4.0 6.0 8.0) - (5.0 10.0 3.0)))) - ( (1.0 m31 1.0) . ,(mx% '((9.0) - (4.0) - (5.0)))) - ( (1.0 m13 1.0) . ,(mx% '((9.0 2.0 7.0)))) - ( (1.0 m11 1.0) . ,(mx% '((9.0)))))) - ;; mx-rowsum (check 'mx-rowsum `(((m33) . ,(mx% '((15.0) (15.0) @@ -841,6 +815,16 @@ ((m11) . ,m11) ((0.0) . ,'xfail))) +;; mx-axpy +(check 'mx-axpy `(( (0.0 m33 m33) . ,m33) + ((-1.0 m33 m33) . ,(mx 3 3 0.0)) + ( (0.0 m31 m31) . ,m31) + ((-1.0 m31 m31) . ,(mx 3 1 0.0)) + ( (0.0 m13 m13) . ,m13) + ((-1.0 m13 m13) . ,(mx 1 3 0.0)) + ( (0.0 m11 m11) . ,m11) + ((-1.0 m11 m11) . ,(mx 1 1 0.0)))) + ;; mx-sympart (check 'mx-sympart `(((m33) . ,(mx% '((8.0 2.0 5.0) (2.0 5.0 8.0) @@ -1032,4 +1016,3 @@ (mx-export "test.csv" m33) (check 'mx-load `((("test.mx") . ,m33))) - diff --git a/tests/test-utils.scm b/tests/test-utils.scm index 9c871e3..404c3d2 100644 --- a/tests/test-utils.scm +++ b/tests/test-utils.scm @@ -1,7 +1,7 @@ ;;;; test-utils.scm ;;@project: matrico (numerical-schemer.xyz) -;;@version: 0.4 (2023-06-01) +;;@version: 0.5 (2023-06-06) ;;@authors: Christian Himpe (0000-0003-2194-6754) ;;@license: zlib-acknowledgement (spdx.org/licenses/zlib-acknowledgement.html) ;;@summary: utils module unit tests @@ -44,7 +44,6 @@ ( (zero? '(1)) . #f) ( (zero? '()) . #t))) - ;; factorial (check 'factorial '(( (0) . 1) ( (1) . 1) @@ -72,4 +71,3 @@ (( 0 -1) . 0) ((-1 0) . 0) ((-1 -1) . 0))) - diff --git a/version.scm b/version.scm index 9bacaa5..64a0a57 100644 --- a/version.scm +++ b/version.scm @@ -1,2 +1 @@ -(define (matrico-version) "0.4rel") - +(define (matrico-version) "0.5rc1")