diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ce6d0e4..3fbcf92 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -50,6 +50,7 @@ jobs: echo "$filters" >> $GITHUB_OUTPUT echo "EOF" >> $GITHUB_OUTPUT + changes: name: Detect Changes runs-on: ubuntu-latest @@ -70,6 +71,7 @@ jobs: select(. != "integration")]') echo "crates=$crates" >> $GITHUB_OUTPUT + required: name: Required Checks runs-on: ubuntu-latest @@ -106,14 +108,21 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: Build + env: + LIBRARY_PATH: /usr/lib/x86_64-linux-gnu/openblas-pthread:/usr/lib/x86_64-linux-gnu/lapack run: cargo build --workspace --all-targets + lint: name: 'Lint, Format and Docs - Workspace' runs-on: ubuntu-latest steps: - name: Checkout uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: Run Clippy run: cargo clippy --workspace --all-targets -- -D warnings - name: Check Formatting @@ -122,6 +131,7 @@ jobs: run: >- RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps --document-private-items + test: name: Test - Crate runs-on: ubuntu-latest @@ -134,8 +144,16 @@ jobs: crate: '${{ fromJSON(needs.changes.outputs.crates) }}' steps: - uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: 'Run ${{ matrix.crate }} tests' + env: + LIBRARY_PATH: /usr/lib/x86_64-linux-gnu/openblas-pthread:/usr/lib/x86_64-linux-gnu/lapack run: 'cargo test -p ${{ matrix.crate }}' + + # TODO: Not needed for now, as the integration tests are minor. + # This should be added back once there are more substantial integration tests, and the project + # is more mature. integration-test: name: Integration Test - Crate runs-on: ubuntu-latest @@ -149,5 +167,6 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 - - name: Run integration tests - run: timeout 240s cargo netlib -- --nocapture + # TODO: No-op + # - name: Run integration tests + # run: timeout 240s cargo netlib -- --nocapture diff --git a/Cargo.lock b/Cargo.lock index 4e6301e..d0a712b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,12 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + [[package]] name = "anstream" version = "0.6.21" @@ -53,13 +59,10 @@ dependencies = [ ] [[package]] -name = "approx" -version = "0.5.1" +name = "anyhow" +version = "1.0.102" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" -dependencies = [ - "num-traits", -] +checksum = "7f202df86484c868dbad7eaa557ef785d5c66295e41b460ef922eca0723b842c" [[package]] name = "atomic-waker" @@ -69,9 +72,9 @@ checksum = "1505bd5d3d116872e7271a6d4e16d81d0c8570876c8de68093a09ac269d8aac0" [[package]] name = "autocfg" -version = "1.5.0" +version = "1.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" [[package]] name = "aws-lc-rs" @@ -113,18 +116,32 @@ version = "3.19.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510" -[[package]] -name = "bytemuck" -version = "1.25.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c8efb64bd706a16a1bdde310ae86b351e4d21550d98d056f22f8a7f7a2183fec" - [[package]] name = "bytes" version = "1.11.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1e748733b7cbc798e1434b6ac524f0c1ff2ab456fe201501e6497c8417a4fc33" +[[package]] +name = "cblas" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7c834e5b5c32a6dbcdb917788d0c38473e1ae7076ff2f4011eb96c5fec20eef" +dependencies = [ + "cblas-sys", + "libc", + "num-complex", +] + +[[package]] +name = "cblas-sys" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8354a31ad0d37f735e74c1973343c99dde4461ea8626a2b8d4fea9e06b987255" +dependencies = [ + "libc", +] + [[package]] name = "cc" version = "1.2.55" @@ -249,7 +266,10 @@ dependencies = [ name = "cnvx-math" version = "0.0.1" dependencies = [ - "nalgebra", + "cblas", + "lapacke", + "openblas-src", + "pkg-config", ] [[package]] @@ -339,6 +359,36 @@ version = "0.8.7" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "dirs" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3e8aa94d75141228480295a7d0e7feb620b1a5ad9f12bc40be62411e38cce4e" +dependencies = [ + "dirs-sys", +] + +[[package]] +name = "dirs-sys" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e01a3366d27ee9890022452ee61b2b63a67e6f13f58900b651ff5665f0bb1fab" +dependencies = [ + "libc", + "option-ext", + "redox_users", + "windows-sys 0.61.2", +] + [[package]] name = "displaydoc" version = "0.2.5" @@ -387,12 +437,32 @@ version = "2.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37909eebbb50d72f9059c3b6d82c0463f2ff062c9e95845c43a6c9c0355411be" +[[package]] +name = "filetime" +version = "0.2.29" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c287a33c7f0a620c38e641e7f60827713987b3c0f26e8ddc9462cc69cf75759" +dependencies = [ + "cfg-if", + "libc", +] + [[package]] name = "find-msvc-tools" version = "0.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582" +[[package]] +name = "flate2" +version = "1.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + [[package]] name = "fnv" version = "1.0.7" @@ -491,114 +561,6 @@ dependencies = [ "wasm-bindgen", ] -[[package]] -name = "glam" -version = "0.14.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "333928d5eb103c5d4050533cec0384302db6be8ef7d3cebd30ec6a35350353da" - -[[package]] -name = "glam" -version = "0.15.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3abb554f8ee44336b72d522e0a7fe86a29e09f839a36022fa869a7dfe941a54b" - -[[package]] -name = "glam" -version = "0.16.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4126c0479ccf7e8664c36a2d719f5f2c140fbb4f9090008098d2c291fa5b3f16" - -[[package]] -name = "glam" -version = "0.17.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e01732b97afd8508eee3333a541b9f7610f454bb818669e66e90f5f57c93a776" - -[[package]] -name = "glam" -version = "0.18.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "525a3e490ba77b8e326fb67d4b44b4bd2f920f44d4cc73ccec50adc68e3bee34" - -[[package]] -name = "glam" -version = "0.19.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b8509e6791516e81c1a630d0bd7fbac36d2fa8712a9da8662e716b52d5051ca" - -[[package]] -name = "glam" -version = "0.20.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f43e957e744be03f5801a55472f593d43fabdebf25a4585db250f04d86b1675f" - -[[package]] -name = "glam" -version = "0.21.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "518faa5064866338b013ff9b2350dc318e14cc4fcd6cb8206d7e7c9886c98815" - -[[package]] -name = "glam" -version = "0.22.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "12f597d56c1bd55a811a1be189459e8fad2bbc272616375602443bdfb37fa774" - -[[package]] -name = "glam" -version = "0.23.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e4afd9ad95555081e109fe1d21f2a30c691b5f0919c67dfa690a2e1eb6bd51c" - -[[package]] -name = "glam" -version = "0.24.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" - -[[package]] -name = "glam" -version = "0.25.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" - -[[package]] -name = "glam" -version = "0.27.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9e05e7e6723e3455f4818c7b26e855439f7546cf617ef669d1adedb8669e5cb9" - -[[package]] -name = "glam" -version = "0.28.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "779ae4bf7e8421cf91c0b3b64e7e8b40b862fba4d393f59150042de7c4965a94" - -[[package]] -name = "glam" -version = "0.29.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee" - -[[package]] -name = "glam" -version = "0.30.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9" - -[[package]] -name = "glam" -version = "0.31.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "556f6b2ea90b8d15a74e0e7bb41671c9bdf38cd9f78c284d750b9ce58a2b5be7" - -[[package]] -name = "glam" -version = "0.32.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f70749695b063ecbf6b62949ccccde2e733ec3ecbbd71d467dca4e5c6c97cca0" - [[package]] name = "h2" version = "0.4.13" @@ -914,12 +876,41 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "lapacke" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8425aee3cfc69f1e907f8487a291d1ddebfd0db2d74493b4e443be0618648744" +dependencies = [ + "lapacke-sys", + "libc", + "num-complex", +] + +[[package]] +name = "lapacke-sys" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2f7d0817c6f4a6029f3b153de01d6498dcf9df659a7536c58bd8df5cd3ccaa6e" +dependencies = [ + "libc", +] + [[package]] name = "libc" version = "0.2.181" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "459427e2af2b9c839b132acb702a1c654d95e10f8c326bfc2ad11310e458b1c5" +[[package]] +name = "libredox" +version = "0.1.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f02ab6bace2054fb888a3c16f990117b579d14a3088e472d63c6011fa185c9d3" +dependencies = [ + "libc", +] + [[package]] name = "linux-raw-sys" version = "0.11.0" @@ -944,16 +935,6 @@ version = "0.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "112b39cec0b298b6c1999fee3e31427f74f676e4cb9879ed1a121b43661a4154" -[[package]] -name = "matrixmultiply" -version = "0.3.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08" -dependencies = [ - "autocfg", - "rawpointer", -] - [[package]] name = "memchr" version = "2.8.0" @@ -966,6 +947,16 @@ version = "0.3.17" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6877bb514081ee2a7ff5ef9de3281f14a4dd4bceac4c09388074a6b5df8a139a" +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + [[package]] name = "mio" version = "1.1.1" @@ -977,61 +968,6 @@ dependencies = [ "windows-sys 0.61.2", ] -[[package]] -name = "nalgebra" -version = "0.34.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "df76ea0ff5c7e6b88689085804d6132ded0ddb9de5ca5b8aeb9eeadc0508a70a" -dependencies = [ - "approx", - "glam 0.14.0", - "glam 0.15.2", - "glam 0.16.0", - "glam 0.17.3", - "glam 0.18.0", - "glam 0.19.0", - "glam 0.20.5", - "glam 0.21.3", - "glam 0.22.0", - "glam 0.23.0", - "glam 0.24.2", - "glam 0.25.0", - "glam 0.27.0", - "glam 0.28.0", - "glam 0.29.3", - "glam 0.30.10", - "glam 0.31.1", - "glam 0.32.1", - "matrixmultiply", - "nalgebra-macros", - "num-complex", - "num-rational", - "num-traits", - "simba", - "typenum", -] - -[[package]] -name = "nalgebra-macros" -version = "0.3.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "num-bigint" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" -dependencies = [ - "num-integer", - "num-traits", -] - [[package]] name = "num-complex" version = "0.4.6" @@ -1041,26 +977,6 @@ dependencies = [ "num-traits", ] -[[package]] -name = "num-integer" -version = "0.1.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-rational" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" -dependencies = [ - "num-bigint", - "num-integer", - "num-traits", -] - [[package]] name = "num-traits" version = "0.2.19" @@ -1082,6 +998,32 @@ version = "1.70.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "384b8ab6d37215f3c5301a95a4accb5d64aa607f1fcb26a11b5303878451b4fe" +[[package]] +name = "openblas-build" +version = "0.10.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb9c85e9e7dd5acdc67b9f3f0c99656b550df716bc63540c6a224a920754a5c2" +dependencies = [ + "anyhow", + "cc", + "flate2", + "tar", + "thiserror 2.0.18", + "ureq", +] + +[[package]] +name = "openblas-src" +version = "0.10.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1a81a5e467f1861ad6ac32d5ec1690ad097d19854753b7424250fe27da46b98" +dependencies = [ + "dirs", + "openblas-build", + "pkg-config", + "vcpkg", +] + [[package]] name = "openssl-probe" version = "0.2.1" @@ -1089,10 +1031,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7c87def4c32ab89d880effc9e097653c8da5d6ef28e6b539d313baaacfbafcbe" [[package]] -name = "paste" -version = "1.0.15" +name = "option-ext" +version = "0.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" +checksum = "04744f49eae99ab78e0d5c0b603ab218f515ea8cfe5a456d7629ad883a3b6e7d" [[package]] name = "percent-encoding" @@ -1112,6 +1054,12 @@ version = "0.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8b870d8c151b6f2fb93e84a13146138f05d02ed11c7e7c54f8826aaaf7c9f184" +[[package]] +name = "pkg-config" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19f132c84eca552bf34cab8ec81f1c1dcc229b811638f9d283dceabe58c5569e" + [[package]] name = "portable-atomic" version = "1.13.1" @@ -1304,10 +1252,15 @@ dependencies = [ ] [[package]] -name = "rawpointer" -version = "0.2.1" +name = "redox_users" +version = "0.5.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" +checksum = "a4e608c6638b9c18977b00b475ac1f28d14e84b27d8d42f70e0bf1e3dec127ac" +dependencies = [ + "getrandom 0.2.17", + "libredox", + "thiserror 2.0.18", +] [[package]] name = "reqwest" @@ -1389,7 +1342,9 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c665f33d38cea657d9614f766881e4d510e0eda4239891eea56b4cadcf01801b" dependencies = [ "aws-lc-rs", + "log", "once_cell", + "ring", "rustls-pki-types", "rustls-webpki", "subtle", @@ -1463,15 +1418,6 @@ version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d" -[[package]] -name = "safe_arch" -version = "0.7.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" -dependencies = [ - "bytemuck", -] - [[package]] name = "same-file" version = "1.0.6" @@ -1572,17 +1518,10 @@ dependencies = [ ] [[package]] -name = "simba" -version = "0.9.1" +name = "simd-adler32" +version = "0.3.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95" -dependencies = [ - "approx", - "num-complex", - "num-traits", - "paste", - "wide", -] +checksum = "703d5c7ef118737c72f1af64ad2f6f8c5e1921f818cdcb97b8fe6fc69bf66214" [[package]] name = "slab" @@ -1676,6 +1615,17 @@ dependencies = [ "libc", ] +[[package]] +name = "tar" +version = "0.4.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f6221d9a6003c78398e3b239969f352578258df48c8eb051caadae0015bc840" +dependencies = [ + "filetime", + "libc", + "xattr", +] + [[package]] name = "target-lexicon" version = "0.13.5" @@ -1900,12 +1850,6 @@ version = "0.2.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e421abadd41a4225275504ea4d6566923418b7f05506fbc9c0fe86ba7396114b" -[[package]] -name = "typenum" -version = "1.19.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "562d481066bde0658276a35467c4af00bdc6ee726305698a55b86e61d7ad82bb" - [[package]] name = "unicode-ident" version = "1.0.22" @@ -1918,6 +1862,34 @@ version = "0.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8ecb6da28b8a351d773b68d5825ac39017e680750f980f3a1a85cd8dd28a47c1" +[[package]] +name = "ureq" +version = "3.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dea7109cdcd5864d4eeb1b58a1648dc9bf520360d7af16ec26d0a9354bafcfc0" +dependencies = [ + "base64", + "log", + "percent-encoding", + "rustls", + "rustls-pki-types", + "ureq-proto", + "utf8-zero", + "webpki-roots", +] + +[[package]] +name = "ureq-proto" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e994ba84b0bd1b1b0cf92878b7ef898a5c1760108fe7b6010327e274917a808c" +dependencies = [ + "base64", + "http", + "httparse", + "log", +] + [[package]] name = "url" version = "2.5.8" @@ -1930,6 +1902,12 @@ dependencies = [ "serde", ] +[[package]] +name = "utf8-zero" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b8c0a043c9540bae7c578c88f91dda8bd82e59ae27c21baca69c8b191aaf5a6e" + [[package]] name = "utf8_iter" version = "1.0.4" @@ -1942,6 +1920,12 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" +[[package]] +name = "vcpkg" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "accd4ea62f7bb7a82fe23066fb0957d48ef677f6eeb8215f372f52e48bb32426" + [[package]] name = "venial" version = "0.5.0" @@ -2086,13 +2070,12 @@ dependencies = [ ] [[package]] -name = "wide" -version = "0.7.33" +name = "webpki-roots" +version = "1.0.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +checksum = "bf85cb06032201fa7c6f829d7db5a7e5aa45bcc0655327713065f6f0576731bf" dependencies = [ - "bytemuck", - "safe_arch", + "rustls-pki-types", ] [[package]] @@ -2373,6 +2356,16 @@ version = "0.6.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9edde0db4769d2dc68579893f2306b26c6ecfbe0ef499b013d731b7b9247e0b9" +[[package]] +name = "xattr" +version = "1.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32e45ad4206f6d2479085147f02bc2ef834ac85886624a23575ae137c8aa8156" +dependencies = [ + "libc", + "rustix", +] + [[package]] name = "yoke" version = "0.8.1" diff --git a/Cargo.toml b/Cargo.toml index de768ee..83ab4d0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,7 +47,13 @@ cnvx-math = { path = "crates/cnvx-math", version = "0.0.1" } cnvx-graph = { path = "crates/cnvx-graph", version = "0.0.1" } cnvx-parse = { path = "crates/cnvx-parse", version = "0.0.1" } clap = { version = "4.5.57", features = ["derive"] } -nalgebra = "0.34.2" +cblas = { version = "0.5.0" } +lapacke = { version = "0.5.0" } +openblas-src = { version = "0.10.16", features = [ + "cblas", + "lapacke", + "system", +] } [workspace.lints] diff --git a/bindings/typst/.gitignore b/bindings/typst/.gitignore new file mode 100644 index 0000000..91d782c --- /dev/null +++ b/bindings/typst/.gitignore @@ -0,0 +1 @@ +examples/*.wasm diff --git a/bindings/typst/examples/cnvx.wasm b/bindings/typst/examples/cnvx.wasm index 6ada7b6..14d7d37 100755 Binary files a/bindings/typst/examples/cnvx.wasm and b/bindings/typst/examples/cnvx.wasm differ diff --git a/bindings/typst/examples/power.pdf b/bindings/typst/examples/power.pdf index 24a37ca..390febd 100644 Binary files a/bindings/typst/examples/power.pdf and b/bindings/typst/examples/power.pdf differ diff --git a/bindings/typst/justfile b/bindings/typst/justfile index ef8c2e5..50995b5 100644 --- a/bindings/typst/justfile +++ b/bindings/typst/justfile @@ -1,3 +1,17 @@ +# List tasks available +default: + @just --list --list-prefix " - " + +# Build the project using cargo to create a WebAssembly module build: cargo build --release --target wasm32-unknown-unknown cp ../../target/wasm32-unknown-unknown/release/cnvx.wasm examples/cnvx.wasm + +# Compile an example file +compile file="examples/power.typ": build + typst compile {{ file }} + +# Clean the build artifacts +clean: + cargo clean + rm -rf examples/cnvx.wasm diff --git a/crates/cnvx-lp/src/primal_simplex.rs b/crates/cnvx-lp/src/primal_simplex.rs index 8901962..28f9a37 100644 --- a/crates/cnvx-lp/src/primal_simplex.rs +++ b/crates/cnvx-lp/src/primal_simplex.rs @@ -264,7 +264,7 @@ impl PrimalSimplexState { /// Attempt to directly run phase 2 if the initial basis is feasible. fn try_phase2(&mut self, max_iter: usize, tol: f64) -> Result { let mut bmat = self.build_bmat(); - match self.compute_basic_solution(&mut bmat) { + match self.compute_basic_solution(&bmat) { Ok(xb) if xb.iter().all(|&v| v >= -tol) => { self.x_b = xb; self.remove_artificial_from_basis(&mut bmat, self.a.cols()) @@ -366,19 +366,16 @@ impl PrimalSimplexState { pub fn build_bmat(&self) -> A { let m = self.a.rows(); let mut bmat = A::new(m, m); - for i in 0..m { - for j in 0..m { - bmat.set(i, j, self.a.get(i, self.basis[j])); - } + for j in 0..m { + let col = self.a.get_col(self.basis[j]); + bmat.set_col(j, &col).unwrap(); } bmat } /// Compute the values of the basic variables by solving `B x_B = b`. - pub fn compute_basic_solution(&self, bmat: &mut A) -> Result, String> { - let mut xb = self.b.clone(); - bmat.mldivide(&mut xb).map_err(|e| format!("gauss failed: {e}"))?; - Ok(xb) + pub fn compute_basic_solution(&self, bmat: &A) -> Result, String> { + bmat.mldivide(&self.b).map_err(|e| format!("gauss failed: {e}")) } /// Run the main simplex iteration loop. @@ -423,19 +420,11 @@ impl PrimalSimplexState { /// Compute dual variables for the current basis. fn compute_duals(&self, bmat: &A) -> Result, SolveError> { let m = bmat.rows(); - let mut pi = (0..m).map(|i| self.c[self.basis[i]]).collect::>(); - - let mut bt = A::new(m, m); - for i in 0..m { - for j in 0..m { - bt.set(i, j, bmat.get(j, i)); - } - } - - bt.mldivide(&mut pi) - .map_err(|e| SolveError::Other(format!("dual solve failed: {e}")))?; + let pi_input: Vec = (0..m).map(|i| self.c[self.basis[i]]).collect(); - Ok(pi) + let bt = bmat.transpose(); + bt.mldivide(&pi_input) + .map_err(|e| SolveError::Other(format!("dual solve failed: {e}"))) } /// Choose entering variable using reduced costs. @@ -455,15 +444,12 @@ impl PrimalSimplexState { /// Compute the simplex direction `d = B^{-1} A_j`. fn compute_direction( &self, - bmat: &mut A, + bmat: &A, entering: usize, ) -> Result, SolveError> { - let mut d = (0..bmat.rows()).map(|i| self.a.get(i, entering)).collect::>(); - - bmat.mldivide(&mut d) - .map_err(|e| SolveError::Other(format!("direction solve failed: {e}")))?; - - Ok(d) + let d_input = self.a.get_col(entering); + bmat.mldivide(&d_input) + .map_err(|e| SolveError::Other(format!("direction solve failed: {e}"))) } /// Choose leaving variable using minimum ratio test. @@ -497,9 +483,8 @@ impl PrimalSimplexState { self.basis[leave_row] = entering; self.non_basis[enter_pos] = leaving; - for i in 0..bmat.rows() { - bmat.set(i, leave_row, self.a.get(i, entering)); - } + let col = self.a.get_col(entering); + bmat.set_col(leave_row, &col).unwrap(); } /// Update the current objective value. @@ -552,10 +537,9 @@ impl PrimalSimplexState { self.x_b = b_aug; let mut bmat = A::new(m, m); - for i in 0..m { - for j in 0..m { - bmat.set(i, j, self.a.get(i, self.basis[j])); - } + for j in 0..m { + let col = self.a.get_col(self.basis[j]); + bmat.set_col(j, &col).unwrap(); } (orig_a, orig_c, bmat) @@ -582,9 +566,8 @@ impl PrimalSimplexState { let leaving = self.basis[row]; self.basis[row] = j; self.non_basis[nb_pos] = leaving; - for i in 0..m { - bmat.set(i, row, self.a.get(i, j)); - } + let col = self.a.get_col(j); + bmat.set_col(row, &col).unwrap(); } else if self.x_b[row].abs() > 1e-12 { return Err( "artificial variable left in basis with non-zero value".into() @@ -595,9 +578,8 @@ impl PrimalSimplexState { let leaving = self.basis[row]; self.basis[row] = j; self.non_basis[nb_pos] = leaving; - for i in 0..m { - bmat.set(i, row, self.a.get(i, j)); - } + let col = self.a.get_col(j); + bmat.set_col(row, &col).unwrap(); break; } } diff --git a/crates/cnvx-math/Cargo.toml b/crates/cnvx-math/Cargo.toml index f2918e2..7c67db5 100644 --- a/crates/cnvx-math/Cargo.toml +++ b/crates/cnvx-math/Cargo.toml @@ -10,7 +10,14 @@ keywords = { workspace = true } readme = { workspace = true } [dependencies] -nalgebra = { workspace = true } + +[target.'cfg(not(target_arch = "wasm32"))'.dependencies] +cblas = { workspace = true } +lapacke = { workspace = true } +openblas-src = { workspace = true } # Links with openblas, which provides BLAS and LAPACK implementations + +[build-dependencies] +pkg-config = { version = "0.3.33" } [lints] workspace = true diff --git a/crates/cnvx-math/build.rs b/crates/cnvx-math/build.rs new file mode 100644 index 0000000..4427856 --- /dev/null +++ b/crates/cnvx-math/build.rs @@ -0,0 +1,27 @@ +fn main() { + let target_arch = std::env::var("CARGO_CFG_TARGET_ARCH").unwrap(); + + if cfg!(target_os = "windows") { + println!("cargo:error=Windows is not currently supported."); + std::process::exit(1); + } + + if target_arch == "wasm32" { + println!( + "cargo:warning=WebAssembly target detected, skipping lapack probe, using custom implementations." + ); + } else { + probe_lapack(); + } +} + +fn probe_lapack() { + // Try lapacke first, fall back to lapack. + if pkg_config::probe_library("lapacke").is_err() { + pkg_config::probe_library("lapack") + .expect("could not find lapack or lapacke via pkg-config"); + } + + // Add blas here as a "just in case" + pkg_config::probe_library("blas").expect("could not find blas via pkg-config"); +} diff --git a/crates/cnvx-math/src/lib.rs b/crates/cnvx-math/src/lib.rs index 5769665..91cbeac 100644 --- a/crates/cnvx-math/src/lib.rs +++ b/crates/cnvx-math/src/lib.rs @@ -1,13 +1,16 @@ //! # CNVX Math //! //! Linear algebra utilities for LP solvers and numerical algorithms. -//! Provides matrix types and traits used in simplex computations and -//! other numerical routines. +//! Provides native matrix types and traits used in simplex computations and +//! other numerical routines without relying on external heavy dependencies. //! //! # Modules //! //! - [`matrix`]: Defines [`DenseMatrix`] and the [`Matrix`] trait for linear algebra operations. +#[cfg(not(target_arch = "wasm32"))] +extern crate openblas_src; + pub mod matrix; -pub use matrix::{DenseMatrix, MatrixWrapper as Matrix, SparseMatrix}; +pub use matrix::{DenseMatrix, Matrix, SparseMatrix}; diff --git a/crates/cnvx-math/src/matrix/dense.rs b/crates/cnvx-math/src/matrix/dense.rs index e30e2dd..3180e82 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -1,45 +1,372 @@ -use nalgebra::DMatrix; +use crate::matrix::Matrix; -use crate::matrix::MatrixWrapper; +#[cfg(not(target_arch = "wasm32"))] +use cblas::{Layout, Transpose, dgemm, dgemv}; -#[derive(Debug, Clone)] -pub struct ExposedDenseMatrix { - pub inner: DMatrix, +/// A dense matrix stored in row-major order using a flat 1D vector. +#[derive(Debug, Clone, PartialEq)] +pub struct DenseMatrix { + rows: usize, + cols: usize, + data: Vec, } -impl MatrixWrapper for ExposedDenseMatrix { +impl DenseMatrix { + /// Helper to calculate the 1D index from 2D coordinates. + #[inline(always)] + fn index(&self, row: usize, col: usize) -> usize { + row * self.cols + col + } + + /// Expose the underlying data slice to internal LAPACK solvers + #[cfg(not(target_arch = "wasm32"))] + pub(crate) fn data(&self) -> &[f64] { + &self.data + } +} + +impl Matrix for DenseMatrix { fn new(rows: usize, cols: usize) -> Self { - Self { inner: DMatrix::zeros(rows, cols) } + Self { rows, cols, data: vec![0.0; rows * cols] } } + #[inline(always)] fn rows(&self) -> usize { - self.inner.nrows() + self.rows } + #[inline(always)] fn cols(&self) -> usize { - self.inner.ncols() + self.cols } + #[inline(always)] fn get(&self, row: usize, col: usize) -> f64 { - self.inner[(row, col)] + assert!(row < self.rows && col < self.cols, "Index out of bounds"); + self.data[self.index(row, col)] } + #[inline(always)] fn set(&mut self, row: usize, col: usize, value: f64) { - self.inner[(row, col)] = value; + assert!(row < self.rows && col < self.cols, "Index out of bounds"); + let idx = self.index(row, col); + self.data[idx] = value; } - // TODO: This is a very naive implementation, and will need to be improved for performance depending on the matrix shape and sparsity. - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { - let a = &self.inner; - let b = DMatrix::from_column_slice(rhs.len(), 1, rhs); - match a.clone().lu().solve(&b) { - Some(solution) => { - for i in 0..rhs.len() { - rhs[i] = solution[(i, 0)]; + fn add(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for addition".to_string()); + } + let data = self.data.iter().zip(other.data.iter()).map(|(a, b)| a + b).collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn sub(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for subtraction".to_string()); + } + let data = self.data.iter().zip(other.data.iter()).map(|(a, b)| a - b).collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn mul(&self, other: &Self) -> Result { + if self.cols != other.rows { + return Err("Incompatible dimensions for matrix multiplication".to_string()); + } + + let mut result = Self::new(self.rows, other.cols); + + #[cfg(not(target_arch = "wasm32"))] + unsafe { + let m = self.rows as i32; + let k = self.cols as i32; + let n = other.cols as i32; + + dgemm( + Layout::RowMajor, + Transpose::None, + Transpose::None, + m, + n, + k, + 1.0, + &self.data, + k, // lda + &other.data, + n, // ldb + 0.0, + &mut result.data, + n, // ldc + ); + } + + #[cfg(target_arch = "wasm32")] + { + for i in 0..self.rows { + for j in 0..other.cols { + let mut sum = 0.0; + for k in 0..self.cols { + sum += self.get(i, k) * other.get(k, j); + } + result.set(i, j, sum); } - Ok(()) } - None => Err("Matrix is singular or not square".to_string()), + } + + Ok(result) + } + + #[cfg(not(target_arch = "wasm32"))] + fn mul_vec(&self, rhs: &[f64]) -> Result, String> { + if self.cols != rhs.len() { + return Err("Matrix columns must match vector length".to_string()); + } + + let m = self.rows as i32; + let n = self.cols as i32; + let mut result = vec![0.0; self.rows]; + + unsafe { + dgemv( + Layout::RowMajor, + Transpose::None, + m, + n, + 1.0, + &self.data, + n, // lda + rhs, + 1, // incx + 0.0, + &mut result, + 1, // incy + ); + } + Ok(result) + } + + #[cfg(target_arch = "wasm32")] + fn mul_vec(&self, rhs: &[f64]) -> Result, String> { + if self.cols != rhs.len() { + return Err("Matrix columns must match vector length".to_string()); + } + + let mut result = vec![0.0; self.rows]; + for i in 0..self.rows { + let mut sum = 0.0; + for j in 0..self.cols { + sum += self.get(i, j) * rhs[j]; + } + result[i] = sum; + } + Ok(result) + } + + fn add_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x + scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn sub_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x - scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn mul_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x * scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn div_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x / scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn mul_elementwise(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for element-wise multiplication" + .to_string()); + } + let data = self + .data + .iter() + .zip(other.data.iter()) + .map(|(&a, &b)| a * b) + .collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn get_row(&self, row: usize) -> Vec { + assert!(row < self.rows, "Row index out of bounds"); + let start = row * self.cols; + self.data[start..start + self.cols].to_vec() + } + + fn get_col(&self, col: usize) -> Vec { + assert!(col < self.cols, "Column index out of bounds"); + (0..self.rows).map(|i| self.data[i * self.cols + col]).collect() + } + + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String> { + if row >= self.rows { + return Err("Row index out of bounds".to_string()); + } + if vec.len() != self.cols { + return Err("Vector length must match matrix columns".to_string()); + } + let start = row * self.cols; + self.data[start..start + self.cols].copy_from_slice(vec); + Ok(()) + } + + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String> { + if col >= self.cols { + return Err("Column index out of bounds".to_string()); + } + if vec.len() != self.rows { + return Err("Vector length must match matrix rows".to_string()); + } + for (i, &val) in vec.iter().enumerate() { + self.data[i * self.cols + col] = val; + } + Ok(()) + } + + fn swap_rows(&mut self, row1: usize, row2: usize) { + assert!(row1 < self.rows && row2 < self.rows, "Row index out of bounds"); + if row1 == row2 { + return; + } + let (r1, r2) = if row1 < row2 { (row1, row2) } else { (row2, row1) }; + let (first, second) = self.data.split_at_mut(r2 * self.cols); + + let row1_slice = &mut first[r1 * self.cols..(r1 + 1) * self.cols]; + let row2_slice = &mut second[0..self.cols]; + row1_slice.swap_with_slice(row2_slice); + } + + fn norm_inf(&self) -> f64 { + let mut max_sum = 0.0_f64; + for i in 0..self.rows { + let start = i * self.cols; + let row_sum: f64 = + self.data[start..start + self.cols].iter().map(|&x| x.abs()).sum(); + if row_sum > max_sum { + max_sum = row_sum; + } + } + max_sum + } + + fn transpose(&self) -> Self { + let mut result = Self::new(self.cols, self.rows); + for i in 0..self.rows { + for j in 0..self.cols { + result.set(j, i, self.get(i, j)); + } + } + result + } + + fn identity(size: usize) -> Self { + let mut m = Self::new(size, size); + for i in 0..size { + m.set(i, i, 1.0); + } + m + } + + fn from_diagonal(vec: &[f64]) -> Self { + let size = vec.len(); + let mut m = Self::new(size, size); + for (i, &val) in vec.iter().enumerate() { + m.data[i * size + i] = val; // More efficient than set() for diagonal + } + m + } + + fn diagonal(&self) -> Vec { + let size = self.rows.min(self.cols); + (0..size).map(|i| self.get(i, i)).collect() + } + + fn mldivide(&self, rhs: &[f64]) -> Result, String> { + // solvers::mldivide_dense(self, rhs) + #[cfg(not(target_arch = "wasm32"))] + { + use crate::matrix::solvers; + #[allow(clippy::needless_return)] + return solvers::mldivide_dense(self, rhs); + } + + #[cfg(target_arch = "wasm32")] + { + if self.rows != self.cols { + return Err("Matrix must be square to solve linear systems".to_string()); + } + if self.rows != rhs.len() { + return Err("RHS vector length must match matrix dimensions".to_string()); + } + + let n = self.rows; + let mut lu = self.data.clone(); + let mut p: Vec = (0..n).collect(); + let b = rhs.to_vec(); + + // LU Decomposition with Partial Pivoting + for i in 0..n { + let mut max_a = 0.0; + let mut imax = i; + + for k in i..n { + let abs_a = lu[k * n + i].abs(); + if abs_a > max_a { + max_a = abs_a; + imax = k; + } + } + + if max_a < 1e-12 { + return Err("Matrix is singular or nearly singular".to_string()); + } + + if imax != i { + for k in 0..n { + lu.swap(i * n + k, imax * n + k); + } + p.swap(i, imax); + } + + for j in (i + 1)..n { + lu[j * n + i] /= lu[i * n + i]; + for k in (i + 1)..n { + lu[j * n + k] -= lu[j * n + i] * lu[i * n + k]; + } + } + } + + // Apply permutation to RHS + let mut x = vec![0.0; n]; + for i in 0..n { + x[i] = b[p[i]]; + } + + // Forward substitution (Ly = Pb) + for i in 0..n { + for k in 0..i { + x[i] -= lu[i * n + k] * x[k]; + } + } + + // Backward substitution (Ux = y) + for i in (0..n).rev() { + for k in (i + 1)..n { + x[i] -= lu[i * n + k] * x[k]; + } + x[i] /= lu[i * n + i]; + } + + Ok(x) } } } diff --git a/crates/cnvx-math/src/matrix/mod.rs b/crates/cnvx-math/src/matrix/mod.rs index d878895..c46597d 100644 --- a/crates/cnvx-math/src/matrix/mod.rs +++ b/crates/cnvx-math/src/matrix/mod.rs @@ -1,27 +1,19 @@ mod dense; mod sparse; -pub use dense::ExposedDenseMatrix as DenseMatrix; -pub use sparse::ExposedSparseMatrix as SparseMatrix; +#[cfg(not(target_arch = "wasm32"))] +mod solvers; // Only include LAPACK-based solvers on non-WASM targets + +pub use dense::DenseMatrix; +pub use sparse::SparseMatrix; /// A generic matrix trait for linear algebra operations. /// /// This trait defines the interface for matrices used in the LP solver, -/// including creation, element access, and solving linear systems. -/// -/// Implementors must provide row/column indexing, element access, and -/// a method for solving linear systems. -pub trait MatrixWrapper: Clone { +/// including creation, element access, arithmetic, and solving linear systems. +pub trait Matrix: Clone { /// Create a new matrix with the given number of rows and columns, /// initialized with zeros. - /// - /// # Example - /// ``` - /// # use cnvx_math::{DenseMatrix, Matrix}; - /// let m = DenseMatrix::new(2, 3); - /// assert_eq!(m.rows(), 2); - /// assert_eq!(m.cols(), 3); - /// ``` fn new(rows: usize, cols: usize) -> Self where Self: Sized; @@ -44,29 +36,92 @@ pub trait MatrixWrapper: Clone { /// Panics if `row` or `col` are out of bounds. fn set(&mut self, row: usize, col: usize, value: f64); - /// Solve a square linear system `Ax = rhs` - /// - /// On success, `rhs` is overwritten with the solution vector `x`. - /// - /// # Errors - /// Returns an `Err(String)` if the system cannot be solved, e.g., if the matrix - /// is singular. - /// - /// # Example - /// ``` - /// # use cnvx_math::{DenseMatrix, Matrix}; - /// let mut a = DenseMatrix::new(2, 2); - /// a.set(0, 0, 2.0); - /// a.set(0, 1, 1.0); - /// a.set(1, 0, 1.0); - /// a.set(1, 1, 3.0); - /// - /// let mut rhs = vec![3.0, 7.0]; - /// a.mldivide(&mut rhs).unwrap(); - /// assert!((rhs[0] - 0.4).abs() < 1e-6); - /// assert!((rhs[1] - 2.2).abs() < 1e-6); - /// ``` - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> + /// Adds another matrix to this matrix. + fn add(&self, other: &Self) -> Result + where + Self: Sized; + + /// Subtracts another matrix from this matrix. + fn sub(&self, other: &Self) -> Result + where + Self: Sized; + + /// Multiplies this matrix by another matrix. + fn mul(&self, other: &Self) -> Result + where + Self: Sized; + + /// Multiplies this matrix by a column vector. + fn mul_vec(&self, rhs: &[f64]) -> Result, String>; + + /// Adds a scalar to every element in the matrix. + fn add_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Subtracts a scalar from every element in the matrix. + fn sub_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Multiplies every element in the matrix by a scalar. + fn mul_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Divides every element in the matrix by a scalar. + fn div_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Performs element-wise (Hadamard) multiplication with another matrix. + fn mul_elementwise(&self, other: &Self) -> Result + where + Self: Sized; + + /// Extracts a specific row as a standard vector. + fn get_row(&self, row: usize) -> Vec; + + /// Extracts a specific column as a standard vector. + fn get_col(&self, col: usize) -> Vec; + + /// Overwrites a specific row with a new vector. + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String>; + + /// Overwrites a specific column with a new vector. + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String>; + + /// Swaps two rows in place. + fn swap_rows(&mut self, row1: usize, row2: usize); + + /// Calculates the L_infinity norm (maximum absolute row sum). + fn norm_inf(&self) -> f64; + + /// Returns the transpose of the matrix. + fn transpose(&self) -> Self + where + Self: Sized; + + /// Creates a square identity matrix of the given size. + fn identity(size: usize) -> Self + where + Self: Sized; + + /// Creates a square diagonal matrix from a vector of diagonal elements. + fn from_diagonal(vec: &[f64]) -> Self where Self: Sized; + + /// Extracts the diagonal elements of the matrix as a vector. + fn diagonal(&self) -> Vec; + + /// Solves Ax = b. + /// Returns a dynamically allocated `Vec` to support both square + /// and rectangular (least-squares) solutions. + fn mldivide(&self, rhs: &[f64]) -> Result, String>; +} + +/// Helper function to calculate the L2 (Euclidean) norm of a standard vector. +pub fn vector_norm_l2(vec: &[f64]) -> f64 { + vec.iter().map(|&x| x * x).sum::().sqrt() } diff --git a/crates/cnvx-math/src/matrix/solvers/cholesky.rs b/crates/cnvx-math/src/matrix/solvers/cholesky.rs new file mode 100644 index 0000000..33989fa --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/cholesky.rs @@ -0,0 +1,30 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dposv}; + +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let n = a.rows() as i32; + + let mut a_copy = a.data().to_vec(); + let mut x = b.to_vec(); + + let info = unsafe { + dposv( + Layout::RowMajor, + b'U', // Check Upper triangle + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is not strictly positive definite".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/lu.rs b/crates/cnvx-math/src/matrix/solvers/lu.rs new file mode 100644 index 0000000..40fe348 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/lu.rs @@ -0,0 +1,36 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dgesv}; + +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let n = a.rows() as i32; + + // LAPACK mutates A to store the LU factors, so we must clone it. + let mut a_copy = a.data().to_vec(); + + // LAPACK mutates b to store the solution x. + let mut x = b.to_vec(); + + // Pivot indices array + let mut ipiv = vec![0; n as usize]; + + let info = unsafe { + dgesv( + Layout::RowMajor, + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut ipiv, + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is singular (U(i,i) is exactly zero)".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/mod.rs b/crates/cnvx-math/src/matrix/solvers/mod.rs new file mode 100644 index 0000000..93e92c4 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/mod.rs @@ -0,0 +1,52 @@ +pub mod cholesky; +pub mod lu; +pub mod qr; +pub mod triangular; + +use crate::matrix::{DenseMatrix, Matrix}; + +/// MATLAB-style mldivide dispatcher for full matrices +/// See +pub fn mldivide_dense(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let m = a.rows(); + let n = a.cols(); + + // 1. Is matrix square? If not, use QR for least squares. + if m != n { + return qr::solve(a, b); + } + + if b.len() != n { + return Err("RHS vector length must match matrix rows".to_string()); + } + + // 2. Is it perfectly triangular? + if let Some(tri_type) = triangular::check_structure(a) { + return triangular::solve(a, b, tri_type); + } + + // 3. Is it symmetric? + if is_symmetric(a) { + // Try Cholesky decomposition (Symmetric Positive Definite) + // If it succeeds, we are done. If it fails, fall through to LU. + if let Ok(x) = cholesky::solve(a, b) { + return Ok(x); + } + } + + // 4. Fallback: LU Decomposition with partial pivoting + lu::solve(a, b) +} + +/// Checks if A is symmetric: A_ij == A_ji +fn is_symmetric(a: &DenseMatrix) -> bool { + let n = a.rows(); + for i in 0..n { + for j in (i + 1)..n { + if (a.get(i, j) - a.get(j, i)).abs() > 1e-12 { + return false; + } + } + } + true +} diff --git a/crates/cnvx-math/src/matrix/solvers/qr.rs b/crates/cnvx-math/src/matrix/solvers/qr.rs new file mode 100644 index 0000000..59c30e1 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/qr.rs @@ -0,0 +1,47 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dgels}; + +/// Solves overdetermined rectangular systems (Least Squares) using QR factorization. +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let m = a.rows() as i32; + let n = a.cols() as i32; + + if m < n { + return Err("Underdetermined systems (m < n) are not yet supported.".into()); + } + if b.len() != m as usize { + return Err("RHS vector length must match matrix rows".into()); + } + + let mut a_copy = a.data().to_vec(); + + // LAPACK's dgels requires the RHS array to have a size of max(m, n) + // because it overwrites the RHS with the solution. + let max_mn = std::cmp::max(m, n) as usize; + let mut x = vec![0.0; max_mn]; + x[..b.len()].copy_from_slice(b); + + let info = unsafe { + dgels( + Layout::RowMajor, + b'N', // No transpose + m, + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is rank deficient".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + // The first `n` elements contain the least squares solution + x.truncate(n as usize); + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/triangular.rs b/crates/cnvx-math/src/matrix/solvers/triangular.rs new file mode 100644 index 0000000..9e58dc7 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/triangular.rs @@ -0,0 +1,75 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use cblas::{Diagonal, Layout, Part, Transpose, dtrsv}; + +pub enum TriType { + Upper, + Lower, + Diagonal, +} + +pub fn check_structure(a: &DenseMatrix) -> Option { + let n = a.rows(); + let mut is_upper = true; + let mut is_lower = true; + + for i in 0..n { + for j in 0..n { + let val = a.get(i, j).abs(); + if i > j && val > 1e-12 { + is_upper = false; + } + if i < j && val > 1e-12 { + is_lower = false; + } + } + } + + if is_upper && is_lower { + Some(TriType::Diagonal) + } else if is_upper { + Some(TriType::Upper) + } else if is_lower { + Some(TriType::Lower) + } else { + None + } +} + +pub fn solve(a: &DenseMatrix, b: &[f64], tri_type: TriType) -> Result, String> { + let n = a.rows() as i32; + let mut x = b.to_vec(); + + // Check for exact zeros on the diagonal to avoid NaN pollution before CBLAS + for i in 0..a.rows() { + if a.get(i, i).abs() < 1e-12 { + return Err("Singular matrix".into()); + } + } + + let uplo = match tri_type { + TriType::Upper => Part::Upper, + TriType::Lower => Part::Lower, + TriType::Diagonal => { + // Pure diagonal division is faster than invoking BLAS + for (i, xi) in x.iter_mut().enumerate() { + *xi /= a.get(i, i); + } + return Ok(x); + } + }; + + unsafe { + dtrsv( + Layout::RowMajor, + uplo, + Transpose::None, + Diagonal::Generic, + n, + a.data(), + n, // lda + &mut x, + 1, // incx + ); + } + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/sparse.rs b/crates/cnvx-math/src/matrix/sparse.rs index efe6207..04764f1 100644 --- a/crates/cnvx-math/src/matrix/sparse.rs +++ b/crates/cnvx-math/src/matrix/sparse.rs @@ -1,32 +1,89 @@ -use crate::matrix::MatrixWrapper; +use crate::matrix::Matrix; -// TODO: This will either use sprs or nalgebra_sparse csr format #[derive(Debug, Clone)] -pub struct ExposedSparseMatrix {} +pub struct SparseMatrix {} #[allow(unused)] -impl MatrixWrapper for ExposedSparseMatrix { +impl Matrix for SparseMatrix { fn new(rows: usize, cols: usize) -> Self { unimplemented!() } - fn rows(&self) -> usize { unimplemented!() } - fn cols(&self) -> usize { unimplemented!() } - fn get(&self, row: usize, col: usize) -> f64 { unimplemented!() } - fn set(&mut self, row: usize, col: usize, value: f64) { unimplemented!() } - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { + fn add(&self, other: &Self) -> Result { + unimplemented!() + } + fn sub(&self, other: &Self) -> Result { + unimplemented!() + } + fn mul(&self, other: &Self) -> Result { + unimplemented!() + } + fn mul_vec(&self, rhs: &[f64]) -> Result, String> { + unimplemented!() + } + + fn add_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn sub_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn mul_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn div_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn mul_elementwise(&self, other: &Self) -> Result { + unimplemented!() + } + + fn get_row(&self, row: usize) -> Vec { + unimplemented!() + } + fn get_col(&self, col: usize) -> Vec { + unimplemented!() + } + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String> { + unimplemented!() + } + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String> { + unimplemented!() + } + fn swap_rows(&mut self, row1: usize, row2: usize) { + unimplemented!() + } + + fn norm_inf(&self) -> f64 { + unimplemented!() + } + + fn transpose(&self) -> Self { + unimplemented!() + } + fn identity(size: usize) -> Self { + unimplemented!() + } + fn from_diagonal(vec: &[f64]) -> Self { + unimplemented!() + } + fn diagonal(&self) -> Vec { + unimplemented!() + } + + fn mldivide(&self, rhs: &[f64]) -> Result, String> { unimplemented!() } } diff --git a/crates/cnvx-math/tests/dense_matrix_tests.rs b/crates/cnvx-math/tests/dense_matrix_tests.rs new file mode 100644 index 0000000..3b3bd4d --- /dev/null +++ b/crates/cnvx-math/tests/dense_matrix_tests.rs @@ -0,0 +1,130 @@ +use cnvx_math::{DenseMatrix, Matrix, matrix::vector_norm_l2}; + +// const EPSILON: f64 = 1e-9; + +#[test] +fn test_scalar_operations() { + let mut m = DenseMatrix::new(2, 2); + m.set(0, 0, 1.0); + m.set(0, 1, 2.0); + m.set(1, 0, 3.0); + m.set(1, 1, 4.0); + + let add_m = m.add_scalar(2.0); + assert_eq!(add_m.get(0, 0), 3.0); + assert_eq!(add_m.get(1, 1), 6.0); + + let mul_m = m.mul_scalar(3.0); + assert_eq!(mul_m.get(0, 1), 6.0); + assert_eq!(mul_m.get(1, 0), 9.0); +} + +#[test] +fn test_elementwise_multiplication() { + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 2.0); + a.set(1, 0, 3.0); + a.set(1, 1, 4.0); + + let mut b = DenseMatrix::new(2, 2); + b.set(0, 0, 2.0); + b.set(0, 1, 3.0); + b.set(1, 0, 4.0); + b.set(1, 1, 5.0); + + let c = a.mul_elementwise(&b).unwrap(); + assert_eq!(c.get(0, 0), 2.0); + assert_eq!(c.get(0, 1), 6.0); + assert_eq!(c.get(1, 0), 12.0); + assert_eq!(c.get(1, 1), 20.0); +} + +#[test] +fn test_row_col_manipulations() { + let mut m = DenseMatrix::new(3, 3); + for i in 0..3 { + for j in 0..3 { + m.set(i, j, (i * 3 + j) as f64); + } + } + // Matrix is: + // [0, 1, 2] + // [3, 4, 5] + // [6, 7, 8] + + assert_eq!(m.get_row(1), vec![3.0, 4.0, 5.0]); + assert_eq!(m.get_col(2), vec![2.0, 5.0, 8.0]); + + m.set_row(0, &[9.0, 9.0, 9.0]).unwrap(); + assert_eq!(m.get(0, 1), 9.0); + + m.set_col(1, &[0.0, 0.0, 0.0]).unwrap(); + assert_eq!(m.get(2, 1), 0.0); + + m.swap_rows(1, 2); + assert_eq!(m.get_row(1), vec![6.0, 0.0, 8.0]); + assert_eq!(m.get_row(2), vec![3.0, 0.0, 5.0]); +} + +#[test] +fn test_norms() { + let mut m = DenseMatrix::new(2, 2); + m.set(0, 0, -3.0); + m.set(0, 1, 4.0); // abs sum = 7 + m.set(1, 0, 1.0); + m.set(1, 1, -8.0); // abs sum = 9 + + assert_eq!(m.norm_inf(), 9.0); + + let v = vec![3.0, -4.0, 0.0]; + assert_eq!(vector_norm_l2(&v), 5.0); +} + +#[test] +fn test_constructors() { + let id = DenseMatrix::identity(3); + assert_eq!(id.get(0, 0), 1.0); + assert_eq!(id.get(1, 1), 1.0); + assert_eq!(id.get(0, 1), 0.0); + + let diag = DenseMatrix::from_diagonal(&[1.5, 2.5]); + assert_eq!(diag.rows(), 2); + assert_eq!(diag.cols(), 2); + assert_eq!(diag.get(0, 0), 1.5); + assert_eq!(diag.get(1, 1), 2.5); + + let mut m = DenseMatrix::new(2, 3); + m.set(0, 0, 7.0); + m.set(1, 1, 8.0); + assert_eq!(m.diagonal(), vec![7.0, 8.0]); +} + +// #[test] +// fn test_mldivide_success() { +// let mut a = DenseMatrix::new(2, 2); +// a.set(0, 0, 2.0); +// a.set(0, 1, 1.0); +// a.set(1, 0, 1.0); +// a.set(1, 1, 3.0); +// +// let mut rhs = vec![3.0, 7.0]; +// a.mldivide(&mut rhs).unwrap(); +// +// assert!((rhs[0] - 0.4).abs() < EPSILON); +// assert!((rhs[1] - 2.2).abs() < EPSILON); +// } + +#[test] +fn test_mldivide_singular_matrix() { + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 1.0); // Singular + + let rhs = vec![2.0, 2.0]; + let result = a.mldivide(&rhs); + + assert!(result.is_err()); +} diff --git a/crates/cnvx-math/tests/mldivide_tests.rs b/crates/cnvx-math/tests/mldivide_tests.rs new file mode 100644 index 0000000..ee5032a --- /dev/null +++ b/crates/cnvx-math/tests/mldivide_tests.rs @@ -0,0 +1,60 @@ +use cnvx_math::{DenseMatrix, Matrix}; + +#[test] +fn test_mldivide_triangular() { + // Upper triangular + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 2.0); + a.set(0, 1, 1.0); + a.set(1, 0, 0.0); + a.set(1, 1, 3.0); + let b = vec![5.0, 6.0]; + let x = a.mldivide(&b).unwrap(); + assert_eq!(x, vec![1.5, 2.0]); +} + +#[test] +fn test_mldivide_cholesky() { + // Symmetric Positive Definite + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 4.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 3.0); + let b = vec![1.0, 2.0]; + let x = a.mldivide(&b).unwrap(); + // Validates via Cholesky decomposition automatically + assert!((x[0] - 0.09090909).abs() < 1e-6); + assert!((x[1] - 0.63636363).abs() < 1e-6); +} + +#[test] +fn test_mldivide_lu() { + // General Square (Requires LU Pivot) + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 0.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 0.0); + let b = vec![2.0, 3.0]; + let x = a.mldivide(&b).unwrap(); + assert_eq!(x, vec![3.0, 2.0]); +} + +#[test] +fn test_mldivide_qr_least_squares() { + // 3x2 Overdetermined system (Line fitting y = mx + c) + let mut a = DenseMatrix::new(3, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 1.0); + a.set(1, 0, 2.0); + a.set(1, 1, 1.0); + a.set(2, 0, 3.0); + a.set(2, 1, 1.0); + let b = vec![1.2, 1.9, 3.2]; + let x = a.mldivide(&b).unwrap(); + + // Result x[0] == slope, x[1] == y-intercept + assert!((x[0] - 1.0).abs() < 1e-5); + assert!((x[1] - 0.1).abs() < 1e-5); +} diff --git a/tests/src/netlib.rs b/tests/src/netlib.rs index feb108e..a489f0d 100644 --- a/tests/src/netlib.rs +++ b/tests/src/netlib.rs @@ -144,7 +144,7 @@ fn run_cnvx(mps: &Path) -> Result { } #[test_case("afiro", Some(-4.6475314286E+02))] -// #[test_case("adlittle", Some(2.2549496316E+05))] +#[test_case("adlittle", Some(2.2549496316E+05))] #[test_case("sc50a", Some(-6.4575077059E+01))] #[test_case("sc50b", Some(-7.0000000000E+01))] // #[test_case("sc105", Some(-5.2202061212E+01))]