diff --git a/.github/release.yml b/.github/release.yml
new file mode 100644
index 00000000..570c4fc0
--- /dev/null
+++ b/.github/release.yml
@@ -0,0 +1,26 @@
+changelog:
+ categories:
+ - title: "Breaking changes"
+ labels:
+ - "breaking change"
+ - title: "New features"
+ labels:
+ - "new feature"
+ - title: "Bug fixes"
+ labels:
+ - bug
+ - title: Changes
+ labels:
+ - change
+ - title: Deprecated
+ labels:
+ - deprecate
+ - title: Removed
+ labels:
+ - remove
+ - title: "Update Dependencies"
+ labels:
+ - dependencies
+ - title: Others
+ labels:
+ - "*"
diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml
index e63fc81f..4123d4a6 100644
--- a/.github/workflows/gh-pages.yml
+++ b/.github/workflows/gh-pages.yml
@@ -1,21 +1,53 @@
+# Based on starter workflow
+# https://github.com/actions/starter-workflows/blob/8217436fdee2338da2d6fd02b7c9fcff634c40e7/pages/static.yml
+#
+# Simple workflow for deploying static content to GitHub Pages
name: "GitHub Pages"
on:
+ # Runs on pushes targeting the default branch
push:
branches:
- master
+ # Allows you to run this workflow manually from the Actions tab
+ workflow_dispatch:
+
+# Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages
+permissions:
+ contents: read
+ pages: write
+ id-token: write
+
+# Allow one concurrent deployment
+concurrency:
+ group: "pages"
+ cancel-in-progress: true
+
jobs:
- pages:
- runs-on: ubuntu-18.04
+ # Single deploy job since we're just deploying
+ deploy:
+ environment:
+ name: github-pages
+ url: ${{ steps.deployment.outputs.page_url }}
+ runs-on: ubuntu-latest
steps:
- - uses: actions/checkout@v2
- - name: Generate code coverage
- run: |
- RUSTDOCFLAGS="--html-in-header katex-header.html" cargo doc --no-deps
- mv target/doc public
- - name: Deploy GitHub Pages
- uses: peaceiris/actions-gh-pages@v3
+ - name: Checkout
+ uses: actions/checkout@v4
+
+ # Generate cargo-doc
+ - name: Generate documentation
+ run: cargo doc --no-deps
+
+ - name: Setup Pages
+ uses: actions/configure-pages@v5
+
+ # Upload target/doc directory
+ - name: Upload artifact
+ uses: actions/upload-pages-artifact@v3
with:
- github_token: ${{ secrets.GITHUB_TOKEN }}
- publish_dir: ./public
+ path: 'target/doc'
+
+ - name: Deploy to GitHub Pages
+ id: deployment
+ uses: actions/deploy-pages@v4
diff --git a/.github/workflows/intel-mkl.yml b/.github/workflows/intel-mkl.yml
index e7e06f2f..14a0b443 100644
--- a/.github/workflows/intel-mkl.yml
+++ b/.github/workflows/intel-mkl.yml
@@ -5,43 +5,16 @@ on:
branches:
- master
pull_request: {}
+ workflow_dispatch:
jobs:
- windows:
- runs-on: windows-2019
+ intel-mkl:
+ strategy:
+ fail-fast: false
+ matrix:
+ system: [ubuntu-22.04, windows-latest]
+ runs-on: ${{ matrix.system }}
steps:
- - uses: actions/checkout@v1
- - uses: actions-rs/cargo@v1
- with:
- command: test
- args: >
- --manifest-path=ndarray-linalg/Cargo.toml
- --no-default-features
- --features=intel-mkl-static
-
- linux:
- runs-on: ubuntu-18.04
- steps:
- - uses: actions/checkout@v1
- - uses: actions-rs/cargo@v1
- name: cargo test
- with:
- command: test
- args: >
- --manifest-path=ndarray-linalg/Cargo.toml
- --no-default-features
- --features=intel-mkl-static
-
- linux-container:
- runs-on: ubuntu-18.04
- container: rustmath/mkl-rust:1.43.0
- steps:
- - uses: actions/checkout@v1
- - uses: actions-rs/cargo@v1
- name: cargo test
- with:
- command: test
- args: >
- --manifest-path=ndarray-linalg/Cargo.toml
- --no-default-features
- --features=intel-mkl-system
+ - uses: actions/checkout@v4
+ - name: cargo test
+ run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=intel-mkl-static --verbose
diff --git a/.github/workflows/netlib.yml b/.github/workflows/netlib.yml
index f278f0ca..70048503 100644
--- a/.github/workflows/netlib.yml
+++ b/.github/workflows/netlib.yml
@@ -5,6 +5,7 @@ on:
branches:
- master
pull_request: {}
+ workflow_dispatch:
jobs:
linux:
@@ -13,17 +14,12 @@ jobs:
matrix:
feature:
- static
- runs-on: ubuntu-18.04
+ runs-on: ubuntu-22.04
steps:
- - uses: actions/checkout@v1
+ - uses: actions/checkout@v4
- name: apt install gfortran
run: |
sudo apt update
sudo apt install -y gfortran
- - uses: actions-rs/cargo@v1
- with:
- command: test
- args: >
- --manifest-path=ndarray-linalg/Cargo.toml
- --no-default-features
- --features=netlib-${{ matrix.feature }}
+ - name: cargo test
+ run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=netlib-${{ matrix.feature }}
diff --git a/.github/workflows/openblas.yml b/.github/workflows/openblas.yml
index 644bf313..5375be70 100644
--- a/.github/workflows/openblas.yml
+++ b/.github/workflows/openblas.yml
@@ -5,12 +5,11 @@ on:
branches:
- master
pull_request: {}
+ workflow_dispatch:
jobs:
linux:
- runs-on: ubuntu-18.04
- container:
- image: rust
+ runs-on: ubuntu-22.04
strategy:
fail-fast: false
matrix:
@@ -18,20 +17,15 @@ jobs:
- static
- system
steps:
- - uses: actions/checkout@v1
+ - uses: actions/checkout@v4
- name: apt install gfortran
run: |
- apt update
- apt install -y gfortran
+ sudo apt update
+ sudo apt install -y gfortran
- name: Install OpenBLAS by apt
run: |
- apt update
- apt install -y libopenblas-dev
+ sudo apt update
+ sudo apt install -y libopenblas-dev
if: ${{ contains(matrix.feature, 'system') }}
- - uses: actions-rs/cargo@v1
- with:
- command: test
- args: >
- --manifest-path=ndarray-linalg/Cargo.toml
- --no-default-features
- --features=openblas-${{ matrix.feature }}
+ - name: cargo test
+ run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=openblas-${{ matrix.feature }}
diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml
index b1c2bc1f..66ca1553 100644
--- a/.github/workflows/rust.yml
+++ b/.github/workflows/rust.yml
@@ -5,34 +5,50 @@ on:
branches:
- master
pull_request: {}
+ workflow_dispatch:
jobs:
check-format:
- runs-on: ubuntu-18.04
+ runs-on: ubuntu-22.04
steps:
- - uses: actions/checkout@v1
- - uses: actions-rs/cargo@v1
- with:
- command: fmt
- args: -- --check
+ - uses: actions/checkout@v4
+ - name: fmt
+ run: cargo fmt -- --check
+
+ check:
+ runs-on: ubuntu-22.04
+ steps:
+ - uses: actions/checkout@v4
+ - name: cargo check
+ run: cargo check --all-targets
+
+ check-doc:
+ runs-on: ubuntu-22.04
+ steps:
+ - uses: actions/checkout@v4
+ - name: cargo doc
+ run: cargo doc --no-deps
clippy:
- runs-on: ubuntu-18.04
+ runs-on: ubuntu-22.04
steps:
- - uses: actions/checkout@v1
- - uses: actions-rs/cargo@v1
- with:
- command: clippy
+ - uses: actions/checkout@v4
+ - name: cargo clippy
+ run: cargo clippy
coverage:
- runs-on: ubuntu-18.04
+ runs-on: ubuntu-22.04
container:
- image: rustmath/mkl-rust:1.43.0
+ image: xd009642/tarpaulin:develop-nightly
options: --security-opt seccomp=unconfined
steps:
- - uses: actions/checkout@v2
+ - uses: actions/checkout@v4
+ - uses: dtolnay/rust-toolchain@nightly
+ - name: Install Cross
+ uses: taiki-e/install-action@v2
+ with:
+ tool: cargo-tarpaulin
- name: Generate code coverage
- run: |
- cargo tarpaulin --verbose --features=intel-mkl --out Xml --manifest-path=ndarray-linalg/Cargo.toml
+ run: cargo +nightly tarpaulin --features=intel-mkl-static --out xml
- name: Upload to codecov.io
- uses: codecov/codecov-action@v1
+ uses: codecov/codecov-action@v5
diff --git a/CHANGELOG.md b/CHANGELOG.md
deleted file mode 100644
index d15ccf81..00000000
--- a/CHANGELOG.md
+++ /dev/null
@@ -1,140 +0,0 @@
-Unreleased
------------
-
-0.13.0 - 20 Feb 2021
-=====================
-
-https://github.com/rust-ndarray/ndarray-linalg/milestone/5
-
-Updated dependencies
----------------------
-- ndarray 0.14 https://github.com/rust-ndarray/ndarray-linalg/pull/258
-- cauchy 0.3.0 (num-complex 0.3.1, rand 0.7.3), lapack 0.17.0 https://github.com/rust-ndarray/ndarray-linalg/pull/260
-
-### optional dependencies
-
-- openblas-src 0.10.2 https://github.com/rust-ndarray/ndarray-linalg/pull/253
-- intel-mkl-src 0.6.0 https://github.com/rust-ndarray/ndarray-linalg/pull/204
-
-Added
-------
-- Split out `ndarray_linalg::lapack` as "lax" crate https://github.com/rust-ndarray/ndarray-linalg/pull/207
- - cargo-workspace https://github.com/rust-ndarray/ndarray-linalg/pull/209
-
-Changed
---------
-- Dual license, MIT or Apache-2.0 License https://github.com/rust-ndarray/ndarray-linalg/pull/262
-- Revise tests for least-square problem https://github.com/rust-ndarray/ndarray-linalg/pull/227
-- Support static link to LAPACK backend https://github.com/rust-ndarray/ndarray-linalg/pull/204
-- Drop LAPACKE dependence, and rewrite them in Rust (see below) https://github.com/rust-ndarray/ndarray-linalg/pull/206
-- Named record like `C { row: i32, lda: i32 }` instead of enum for `MatrixLayout` https://github.com/rust-ndarray/ndarray-linalg/pull/211
-- Split LAPACK error into computational failure and invalid values https://github.com/rust-ndarray/ndarray-linalg/pull/210
-- Use thiserror crate https://github.com/rust-ndarray/ndarray-linalg/pull/208
-
-### LAPACKE rewrite
-
-- Cholesky https://github.com/rust-ndarray/ndarray-linalg/pull/225
-- Eigenvalue for general matrix https://github.com/rust-ndarray/ndarray-linalg/pull/212
-- Eigenvalue for symmetric/Hermitian matrix https://github.com/rust-ndarray/ndarray-linalg/pull/217
-- least squares problem https://github.com/rust-ndarray/ndarray-linalg/pull/220
-- QR decomposition https://github.com/rust-ndarray/ndarray-linalg/pull/224
-- LU decomposition https://github.com/rust-ndarray/ndarray-linalg/pull/213
-- LDL decomposition https://github.com/rust-ndarray/ndarray-linalg/pull/216
-- SVD https://github.com/rust-ndarray/ndarray-linalg/pull/218
-- SVD divid-and-conquer https://github.com/rust-ndarray/ndarray-linalg/pull/219
-- Tridiagonal https://github.com/rust-ndarray/ndarray-linalg/pull/235
-
-Maintenance
------------
-- Coverage report using codecov https://github.com/rust-ndarray/ndarray-linalg/pull/215
-- Fix for clippy, and add CI check https://github.com/rust-ndarray/ndarray-linalg/pull/205
-
-0.12.1 - 28 June 2020
-======================
-
-Added
-------
-- Tridiagonal matrix support https://github.com/rust-ndarray/ndarray-linalg/pull/196
-- KaTeX support in rustdoc https://github.com/rust-ndarray/ndarray-linalg/pull/202
-- Least square problems https://github.com/rust-ndarray/ndarray-linalg/pull/197
-- LOBPCG solver https://github.com/rust-ndarray/ndarray-linalg/pull/184
-
-Changed
--------
-- Grouping and Plot in benchmark https://github.com/rust-ndarray/ndarray-linalg/pull/200
-- `Clone` trait for `LUFactorized` https://github.com/rust-ndarray/ndarray-linalg/pull/192
-
-Maintenance
------------
-- Fix repository URL https://github.com/rust-ndarray/ndarray-linalg/pull/198
-- Use GitHub Actions instead of Azure Pipeline https://github.com/rust-ndarray/ndarray-linalg/pull/193
-- Test cargo-fmt on CI https://github.com/rust-ndarray/ndarray-linalg/pull/194
-
-0.12.0 - 14 Oct 2019
-====================
-
-Added
------
-- SVD by divide-and-conquer https://github.com/rust-ndarray/ndarray-linalg/pull/164
-- Householder reflection https://github.com/rust-ndarray/ndarray-linalg/pull/154
-- Arnoldi iteration https://github.com/rust-ndarray/ndarray-linalg/pull/155
-
-Changed
-----------
-- Replace `operator::Operator*` traits by new `LinearOperator trait` https://github.com/rust-ndarray/ndarray-linalg/pull/159
-- ndarray 0.13.0 https://github.com/rust-ndarray/ndarray-linalg/pull/172
-- blas-src 0.4.0, lapack-src 0.4.0, openblas-src 0.7.0 https://github.com/rust-ndarray/ndarray-linalg/pull/174
-- restore `static` feature flag
-
-0.11.1 - 12 June 2019
-======================
-
-- Hotfix for document generation https://github.com/rust-ndarray/ndarray-linalg/pull/153
-
-0.11.0 - 12 June 2019
-====================
-
-Added
---------
-- Dependency to cauchy 0.2 https://github.com/rust-ndarray/ndarray-linalg/pull/139
-- `generate::random_{unitary,regular}` for debug use https://github.com/rust-ndarray/ndarray-linalg/pull/140
-- `krylov` submodule
- - modified Gram-Schmit https://github.com/rust-ndarray/ndarray-linalg/pull/149 https://github.com/rust-ndarray/ndarray-linalg/pull/150
- - Krylov subspace methods are not implemented yet.
-
-Removed
-----------
-- `static` feature https://github.com/rust-ndarray/ndarray-linalg/pull/136
- - See README for detail
-- `accelerate` feature https://github.com/rust-ndarray/ndarray-linalg/pull/141
-- Dependencies to derive-new, procedurals
-
-Changed
----------
-- Switch CI service: Circle CI -> Azure Pipeline https://github.com/rust-ndarray/ndarray-linalg/pull/141
-- submodule `lapack_traits` is renamed to https://github.com/rust-ndarray/ndarray-linalg/pull/139
-- `ndarray_linalg::Scalar` trait is split into two parts https://github.com/rust-ndarray/ndarray-linalg/pull/139
- - [cauchy::Scalar](https://docs.rs/cauchy/0.2.0/cauchy/trait.Scalar.html) is a refined real/complex common trait
- - `lapack::Lapack` is a trait for primitive types which LAPACK supports
-- Error type becomes simple https://github.com/rust-ndarray/ndarray-linalg/pull/118 https://github.com/rust-ndarray/ndarray-linalg/pull/127
-- Assertions becomes more verbose https://github.com/rust-ndarray/ndarray-linalg/pull/147
-- blas-src 0.3, lapack-src 0.3
- - intel-mkl-src becomes 0.4, which supports Windows! https://github.com/rust-ndarray/ndarray-linalg/pull/146
-
-0.10.0 - 2 Sep 2018
-===================
-
-Update Dependencies
---------------------
-
-- ndarray 0.12
-- rand 0.5
-- num-complex 0.2
-- openblas-src 0.6
-- lapacke 0.2
-
-See also https://github.com/rust-ndarray/ndarray-linalg/pull/110
-
-Added
-------
-- serde-1 feature gate https://github.com/rust-ndarray/ndarray-linalg/pull/99, https://github.com/rust-ndarray/ndarray-linalg/pull/116
diff --git a/README.md b/README.md
index 4478dca3..85b1428c 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,8 @@
ndarray-linalg
===============
-[](https://crates.io/crates/ndarray-linalg)
+[](https://crates.io/crates/ndarray-linalg)
[](https://docs.rs/ndarray-linalg)
+[](https://rust-ndarray.github.io/ndarray-linalg/ndarray_linalg/index.html)
Linear algebra package for Rust with [ndarray](https://github.com/rust-ndarray/ndarray) based on external LAPACK implementations.
@@ -61,7 +62,7 @@ Supported features are following:
### For library developer
-If you creating a library depending on this crate, we encourage you not to link any backend:
+If you are creating a library depending on this crate, we encourage you not to link any backend:
```toml
[dependencies]
@@ -85,25 +86,6 @@ Only x86_64 system is supported currently.
|Netlib |✔️ |- |- |
|Intel MKL|✔️ |✔️ |✔️ |
-Generate document with KaTeX
-------------------------------
-
-You need to set `RUSTDOCFLAGS` explicitly:
-
-```shell
-RUSTDOCFLAGS="--html-in-header katex-header.html" cargo doc --no-deps
-```
-
-This **only** works for `--no-deps` build because `katex-header.html` does not exists for dependent crates.
-If you wish to set `RUSTDOCFLAGS` automatically in this crate, you can put [.cargo/config](https://doc.rust-lang.org/cargo/reference/config.html):
-
-```toml
-[build]
-rustdocflags = ["--html-in-header", "katex-header.html"]
-```
-
-But, be sure that this works only for `--no-deps`. `cargo doc` will fail with this `.cargo/config`.
-
License
--------
diff --git a/katex-header.html b/katex-header.html
deleted file mode 100644
index 6e10c052..00000000
--- a/katex-header.html
+++ /dev/null
@@ -1,16 +0,0 @@
-
-
-
-
-
diff --git a/lax/CHANGELOG.md b/lax/CHANGELOG.md
new file mode 100644
index 00000000..68e828b4
--- /dev/null
+++ b/lax/CHANGELOG.md
@@ -0,0 +1,15 @@
+Unreleased
+-----------
+
+0.2.0 - 17 July 2021
+=====================
+
+Updated dependencies
+---------------------
+- cauchy 0.4 (num-complex 0.4, rand 0.8), lapack 0.18 https://github.com/rust-ndarray/ndarray-linalg/pull/276
+
+Fixed
+-----
+- Fix memory layout of the output of inverse of LUFactorized https://github.com/rust-ndarray/ndarray-linalg/pull/297
+- Fix Eig for column-major arrays with real elements https://github.com/rust-ndarray/ndarray-linalg/pull/298
+- Fix Solve::solve_h_* for complex inputs with standard layout https://github.com/rust-ndarray/ndarray-linalg/pull/296
diff --git a/lax/Cargo.toml b/lax/Cargo.toml
index 812ad126..076d5cf2 100644
--- a/lax/Cargo.toml
+++ b/lax/Cargo.toml
@@ -1,6 +1,6 @@
[package]
name = "lax"
-version = "0.1.0"
+version = "0.17.0"
authors = ["Toshiki Teramura "]
edition = "2018"
@@ -25,31 +25,29 @@ netlib-system = ["netlib-src/system"]
openblas-static = ["openblas-src/static"]
openblas-system = ["openblas-src/system"]
-intel-mkl-static = ["intel-mkl-src/mkl-static-lp64-seq", "intel-mkl-src/download"]
+intel-mkl-static = ["intel-mkl-src/mkl-static-lp64-seq"]
intel-mkl-system = ["intel-mkl-src/mkl-dynamic-lp64-seq"]
[dependencies]
-thiserror = "1.0.23"
-cauchy = "0.3.0"
+thiserror = "2.0.0"
+cauchy = "0.4.0"
num-traits = "0.2.14"
-lapack = "0.17.0"
+lapack-sys = "0.15.0"
+katexit = "0.1.2"
[dependencies.intel-mkl-src]
-version = "0.6.0"
+version = "0.8.1"
default-features = false
optional = true
[dependencies.netlib-src]
-version = "0.8.0"
+version = "0.9.0"
optional = true
features = ["cblas"]
default-features = false
[dependencies.openblas-src]
-version = "0.10.2"
+version = "0.10.4"
optional = true
default-features = false
features = ["cblas"]
-
-[package.metadata.release]
-no-dev-version = true
diff --git a/lax/README.md b/lax/README.md
index ed563735..9dea6b32 100644
--- a/lax/README.md
+++ b/lax/README.md
@@ -1,6 +1,9 @@
Linear Algebra eXtension (LAX)
===============================
+[](https://crates.io/crates/lax)
+[](https://docs.rs/lax)
+
ndarray-free safe Rust wrapper for LAPACK FFI for implementing ndarray-linalg crate.
This crate responsibles for
diff --git a/lax/src/alloc.rs b/lax/src/alloc.rs
new file mode 100644
index 00000000..63458818
--- /dev/null
+++ b/lax/src/alloc.rs
@@ -0,0 +1,78 @@
+use cauchy::*;
+use std::mem::MaybeUninit;
+
+/// Helper for getting pointer of slice
+pub(crate) trait AsPtr: Sized {
+ type Elem;
+ fn as_ptr(vec: &[Self]) -> *const Self::Elem;
+ fn as_mut_ptr(vec: &mut [Self]) -> *mut Self::Elem;
+}
+
+macro_rules! impl_as_ptr {
+ ($target:ty, $elem:ty) => {
+ impl AsPtr for $target {
+ type Elem = $elem;
+ fn as_ptr(vec: &[Self]) -> *const Self::Elem {
+ vec.as_ptr() as *const _
+ }
+ fn as_mut_ptr(vec: &mut [Self]) -> *mut Self::Elem {
+ vec.as_mut_ptr() as *mut _
+ }
+ }
+ };
+}
+impl_as_ptr!(i32, i32);
+impl_as_ptr!(f32, f32);
+impl_as_ptr!(f64, f64);
+impl_as_ptr!(c32, lapack_sys::__BindgenComplex);
+impl_as_ptr!(c64, lapack_sys::__BindgenComplex);
+impl_as_ptr!(MaybeUninit, i32);
+impl_as_ptr!(MaybeUninit, f32);
+impl_as_ptr!(MaybeUninit, f64);
+impl_as_ptr!(MaybeUninit, lapack_sys::__BindgenComplex);
+impl_as_ptr!(MaybeUninit, lapack_sys::__BindgenComplex);
+
+pub(crate) trait VecAssumeInit {
+ type Elem;
+ unsafe fn assume_init(self) -> Vec;
+
+ /// An replacement of unstable API
+ /// https://doc.rust-lang.org/std/mem/union.MaybeUninit.html#method.slice_assume_init_ref
+ unsafe fn slice_assume_init_ref(&self) -> &[Self::Elem];
+
+ /// An replacement of unstable API
+ /// https://doc.rust-lang.org/std/mem/union.MaybeUninit.html#method.slice_assume_init_mut
+ unsafe fn slice_assume_init_mut(&mut self) -> &mut [Self::Elem];
+}
+
+impl VecAssumeInit for Vec> {
+ type Elem = T;
+ unsafe fn assume_init(self) -> Vec {
+ // FIXME use Vec::into_raw_parts instead after stablized
+ // https://doc.rust-lang.org/std/vec/struct.Vec.html#method.into_raw_parts
+ let mut me = std::mem::ManuallyDrop::new(self);
+ Vec::from_raw_parts(me.as_mut_ptr() as *mut T, me.len(), me.capacity())
+ }
+
+ unsafe fn slice_assume_init_ref(&self) -> &[T] {
+ std::slice::from_raw_parts(self.as_ptr() as *const T, self.len())
+ }
+
+ unsafe fn slice_assume_init_mut(&mut self) -> &mut [T] {
+ std::slice::from_raw_parts_mut(self.as_mut_ptr() as *mut T, self.len())
+ }
+}
+
+/// Create a vector without initialization
+///
+/// Safety
+/// ------
+/// - Memory is not initialized. Do not read the memory before write.
+///
+pub(crate) fn vec_uninit(n: usize) -> Vec> {
+ let mut v = Vec::with_capacity(n);
+ unsafe {
+ v.set_len(n);
+ }
+ v
+}
diff --git a/lax/src/cholesky.rs b/lax/src/cholesky.rs
index 8305efe5..785f6e5e 100644
--- a/lax/src/cholesky.rs
+++ b/lax/src/cholesky.rs
@@ -1,27 +1,25 @@
-//! Cholesky decomposition
+//! Factorize positive-definite symmetric/Hermitian matrices using Cholesky algorithm
use super::*;
use crate::{error::*, layout::*};
use cauchy::*;
-pub trait Cholesky_: Sized {
- /// Cholesky: wrapper of `*potrf`
- ///
- /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.**
+/// Compute Cholesky decomposition according to [UPLO]
+///
+/// LAPACK correspondance
+/// ----------------------
+///
+/// | f32 | f64 | c32 | c64 |
+/// |:-------|:-------|:-------|:-------|
+/// | spotrf | dpotrf | cpotrf | zpotrf |
+///
+pub trait CholeskyImpl: Scalar {
fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
-
- /// Wrapper of `*potri`
- ///
- /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.**
- fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
-
- /// Wrapper of `*potrs`
- fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
}
-macro_rules! impl_cholesky {
- ($scalar:ty, $trf:path, $tri:path, $trs:path) => {
- impl Cholesky_ for $scalar {
+macro_rules! impl_cholesky_ {
+ ($s:ty, $trf:path) => {
+ impl CholeskyImpl for $s {
fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
let (n, _) = l.size();
if matches!(l, MatrixLayout::C { .. }) {
@@ -29,7 +27,7 @@ macro_rules! impl_cholesky {
}
let mut info = 0;
unsafe {
- $trf(uplo as u8, n, a, n, &mut info);
+ $trf(uplo.as_ptr(), &n, AsPtr::as_mut_ptr(a), &n, &mut info);
}
info.as_lapack_result()?;
if matches!(l, MatrixLayout::C { .. }) {
@@ -37,7 +35,30 @@ macro_rules! impl_cholesky {
}
Ok(())
}
+ }
+ };
+}
+impl_cholesky_!(c64, lapack_sys::zpotrf_);
+impl_cholesky_!(c32, lapack_sys::cpotrf_);
+impl_cholesky_!(f64, lapack_sys::dpotrf_);
+impl_cholesky_!(f32, lapack_sys::spotrf_);
+
+/// Compute inverse matrix using Cholesky factroization result
+///
+/// LAPACK correspondance
+/// ----------------------
+///
+/// | f32 | f64 | c32 | c64 |
+/// |:-------|:-------|:-------|:-------|
+/// | spotri | dpotri | cpotri | zpotri |
+///
+pub trait InvCholeskyImpl: Scalar {
+ fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+}
+macro_rules! impl_inv_cholesky {
+ ($s:ty, $tri:path) => {
+ impl InvCholeskyImpl for $s {
fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
let (n, _) = l.size();
if matches!(l, MatrixLayout::C { .. }) {
@@ -45,7 +66,7 @@ macro_rules! impl_cholesky {
}
let mut info = 0;
unsafe {
- $tri(uplo as u8, n, a, l.lda(), &mut info);
+ $tri(uplo.as_ptr(), &n, AsPtr::as_mut_ptr(a), &l.lda(), &mut info);
}
info.as_lapack_result()?;
if matches!(l, MatrixLayout::C { .. }) {
@@ -53,7 +74,30 @@ macro_rules! impl_cholesky {
}
Ok(())
}
+ }
+ };
+}
+impl_inv_cholesky!(c64, lapack_sys::zpotri_);
+impl_inv_cholesky!(c32, lapack_sys::cpotri_);
+impl_inv_cholesky!(f64, lapack_sys::dpotri_);
+impl_inv_cholesky!(f32, lapack_sys::spotri_);
+/// Solve linear equation using Cholesky factroization result
+///
+/// LAPACK correspondance
+/// ----------------------
+///
+/// | f32 | f64 | c32 | c64 |
+/// |:-------|:-------|:-------|:-------|
+/// | spotrs | dpotrs | cpotrs | zpotrs |
+///
+pub trait SolveCholeskyImpl: Scalar {
+ fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
+}
+
+macro_rules! impl_solve_cholesky {
+ ($s:ty, $trs:path) => {
+ impl SolveCholeskyImpl for $s {
fn solve_cholesky(
l: MatrixLayout,
mut uplo: UPLO,
@@ -70,7 +114,16 @@ macro_rules! impl_cholesky {
}
}
unsafe {
- $trs(uplo as u8, n, nrhs, a, l.lda(), b, n, &mut info);
+ $trs(
+ uplo.as_ptr(),
+ &n,
+ &nrhs,
+ AsPtr::as_ptr(a),
+ &l.lda(),
+ AsPtr::as_mut_ptr(b),
+ &n,
+ &mut info,
+ );
}
info.as_lapack_result()?;
if matches!(l, MatrixLayout::C { .. }) {
@@ -82,9 +135,8 @@ macro_rules! impl_cholesky {
}
}
};
-} // end macro_rules
-
-impl_cholesky!(f64, lapack::dpotrf, lapack::dpotri, lapack::dpotrs);
-impl_cholesky!(f32, lapack::spotrf, lapack::spotri, lapack::spotrs);
-impl_cholesky!(c64, lapack::zpotrf, lapack::zpotri, lapack::zpotrs);
-impl_cholesky!(c32, lapack::cpotrf, lapack::cpotri, lapack::cpotrs);
+}
+impl_solve_cholesky!(c64, lapack_sys::zpotrs_);
+impl_solve_cholesky!(c32, lapack_sys::cpotrs_);
+impl_solve_cholesky!(f64, lapack_sys::dpotrs_);
+impl_solve_cholesky!(f32, lapack_sys::spotrs_);
diff --git a/lax/src/eig.rs b/lax/src/eig.rs
index 53245de7..f02035bb 100644
--- a/lax/src/eig.rs
+++ b/lax/src/eig.rs
@@ -1,165 +1,297 @@
-//! Eigenvalue decomposition for general matrices
+//! Eigenvalue problem for general matricies
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32 | f64 | c32 | c64 |
+//! |:------|:------|:------|:------|
+//! | sgeev | dgeev | cgeev | zgeev |
+//!
use crate::{error::*, layout::MatrixLayout, *};
use cauchy::*;
use num_traits::{ToPrimitive, Zero};
-/// Wraps `*geev` for general matrices
-pub trait Eig_: Scalar {
- /// Calculate Right eigenvalue
- fn eig(
- calc_v: bool,
- l: MatrixLayout,
- a: &mut [Self],
- ) -> Result<(Vec, Vec)>;
+#[cfg_attr(doc, katexit::katexit)]
+/// Eigenvalue problem for general matrix
+///
+/// To manage memory more strictly, use [EigWork].
+///
+/// Right and Left eigenvalue problem
+/// ----------------------------------
+/// LAPACK can solve both right eigenvalue problem
+/// $$
+/// AV_R = V_R \Lambda
+/// $$
+/// where $V_R = \left( v_R^1, \cdots, v_R^n \right)$ are right eigenvectors
+/// and left eigenvalue problem
+/// $$
+/// V_L^\dagger A = V_L^\dagger \Lambda
+/// $$
+/// where $V_L = \left( v_L^1, \cdots, v_L^n \right)$ are left eigenvectors
+/// and eigenvalues
+/// $$
+/// \Lambda = \begin{pmatrix}
+/// \lambda_1 & & 0 \\\\
+/// & \ddots & \\\\
+/// 0 & & \lambda_n
+/// \end{pmatrix}
+/// $$
+/// which satisfies $A v_R^i = \lambda_i v_R^i$ and
+/// $\left(v_L^i\right)^\dagger A = \lambda_i \left(v_L^i\right)^\dagger$
+/// for column-major matrices, although row-major matrices are not supported.
+/// Since a row-major matrix can be interpreted
+/// as a transpose of a column-major matrix,
+/// this transforms right eigenvalue problem to left one:
+///
+/// $$
+/// A^\dagger V = V Λ ⟺ V^\dagger A = Λ V^\dagger
+/// $$
+///
+#[non_exhaustive]
+pub struct EigWork {
+ /// Problem size
+ pub n: i32,
+ /// Compute right eigenvectors or not
+ pub jobvr: JobEv,
+ /// Compute left eigenvectors or not
+ pub jobvl: JobEv,
+
+ /// Eigenvalues
+ pub eigs: Vec>,
+ /// Real part of eigenvalues used in real routines
+ pub eigs_re: Option>>,
+ /// Imaginary part of eigenvalues used in real routines
+ pub eigs_im: Option>>,
+
+ /// Left eigenvectors
+ pub vc_l: Option>>,
+ /// Left eigenvectors used in real routines
+ pub vr_l: Option>>,
+ /// Right eigenvectors
+ pub vc_r: Option>>,
+ /// Right eigenvectors used in real routines
+ pub vr_r: Option>>,
+
+ /// Working memory
+ pub work: Vec>,
+ /// Working memory with `T::Real`
+ pub rwork: Option>>,
+}
+
+impl EigWork
+where
+ T: Scalar,
+ EigWork: EigWorkImpl,
+{
+ /// Create new working memory for eigenvalues compution.
+ pub fn new(calc_v: bool, l: MatrixLayout) -> Result {
+ EigWorkImpl::new(calc_v, l)
+ }
+
+ /// Compute eigenvalues and vectors on this working memory.
+ pub fn calc(&mut self, a: &mut [T]) -> Result> {
+ EigWorkImpl::calc(self, a)
+ }
+
+ /// Compute eigenvalues and vectors by consuming this working memory.
+ pub fn eval(self, a: &mut [T]) -> Result> {
+ EigWorkImpl::eval(self, a)
+ }
+}
+
+/// Owned result of eigenvalue problem by [EigWork::eval]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigOwned {
+ /// Eigenvalues
+ pub eigs: Vec,
+ /// Right eigenvectors
+ pub vr: Option>,
+ /// Left eigenvectors
+ pub vl: Option>,
}
-macro_rules! impl_eig_complex {
- ($scalar:ty, $ev:path) => {
- impl Eig_ for $scalar {
- fn eig(
- calc_v: bool,
- l: MatrixLayout,
- mut a: &mut [Self],
- ) -> Result<(Vec, Vec)> {
+/// Reference result of eigenvalue problem by [EigWork::calc]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigRef<'work, T: Scalar> {
+ /// Eigenvalues
+ pub eigs: &'work [T::Complex],
+ /// Right eigenvectors
+ pub vr: Option<&'work [T::Complex]>,
+ /// Left eigenvectors
+ pub vl: Option<&'work [T::Complex]>,
+}
+
+/// Helper trait for implementing [EigWork] methods
+pub trait EigWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(calc_v: bool, l: MatrixLayout) -> Result;
+ fn calc<'work>(&'work mut self, a: &mut [Self::Elem]) -> Result>;
+ fn eval(self, a: &mut [Self::Elem]) -> Result>;
+}
+
+macro_rules! impl_eig_work_c {
+ ($c:ty, $ev:path) => {
+ impl EigWorkImpl for EigWork<$c> {
+ type Elem = $c;
+
+ fn new(calc_v: bool, l: MatrixLayout) -> Result {
let (n, _) = l.size();
- // Because LAPACK assumes F-continious array, C-continious array should be taken Hermitian conjugate.
- // However, we utilize a fact that left eigenvector of A^H corresponds to the right eigenvector of A
let (jobvl, jobvr) = if calc_v {
match l {
- MatrixLayout::C { .. } => (b'V', b'N'),
- MatrixLayout::F { .. } => (b'N', b'V'),
+ MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+ MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
}
} else {
- (b'N', b'N')
+ (JobEv::None, JobEv::None)
};
- let mut eigs = unsafe { vec_uninit(n as usize) };
- let mut rwork = unsafe { vec_uninit(2 * n as usize) };
+ let mut eigs = vec_uninit(n as usize);
+ let mut rwork = vec_uninit(2 * n as usize);
- let mut vl = if jobvl == b'V' {
- Some(unsafe { vec_uninit((n * n) as usize) })
- } else {
- None
- };
- let mut vr = if jobvr == b'V' {
- Some(unsafe { vec_uninit((n * n) as usize) })
- } else {
- None
- };
+ let mut vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let mut vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
// calc work size
let mut info = 0;
- let mut work_size = [Self::zero()];
+ let mut work_size = [<$c>::zero()];
unsafe {
$ev(
- jobvl,
- jobvr,
- n,
- &mut a,
- n,
- &mut eigs,
- &mut vl.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut vr.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut work_size,
- -1,
- &mut rwork,
+ jobvl.as_ptr(),
+ jobvr.as_ptr(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs),
+ AsPtr::as_mut_ptr(vc_l.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(vc_r.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ AsPtr::as_mut_ptr(&mut rwork),
&mut info,
)
};
info.as_lapack_result()?;
- // actal ev
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
+ let work: Vec> = vec_uninit(lwork);
+ Ok(Self {
+ n,
+ jobvl,
+ jobvr,
+ eigs,
+ eigs_re: None,
+ eigs_im: None,
+ rwork: Some(rwork),
+ vc_l,
+ vc_r,
+ vr_l: None,
+ vr_r: None,
+ work,
+ })
+ }
+
+ fn calc<'work>(
+ &'work mut self,
+ a: &mut [Self::Elem],
+ ) -> Result> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
unsafe {
$ev(
- jobvl,
- jobvr,
- n,
- &mut a,
- n,
- &mut eigs,
- &mut vl.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut vr.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut work,
- lwork as i32,
- &mut rwork,
+ self.jobvl.as_ptr(),
+ self.jobvr.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.eigs),
+ AsPtr::as_mut_ptr(self.vc_l.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(self.vc_r.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
&mut info,
)
};
info.as_lapack_result()?;
-
// Hermite conjugate
- if jobvl == b'V' {
- for c in vl.as_mut().unwrap().iter_mut() {
- c.im = -c.im
+ if let Some(vl) = self.vc_l.as_mut() {
+ for value in vl {
+ let value = unsafe { value.assume_init_mut() };
+ value.im = -value.im;
}
}
+ Ok(EigRef {
+ eigs: unsafe { self.eigs.slice_assume_init_ref() },
+ vl: self
+ .vc_l
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ vr: self
+ .vc_r
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ })
+ }
- Ok((eigs, vr.or(vl).unwrap_or(Vec::new())))
+ fn eval(mut self, a: &mut [Self::Elem]) -> Result> {
+ let _eig_ref = self.calc(a)?;
+ Ok(EigOwned {
+ eigs: unsafe { self.eigs.assume_init() },
+ vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+ vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+ })
}
}
};
}
-impl_eig_complex!(c64, lapack::zgeev);
-impl_eig_complex!(c32, lapack::cgeev);
-
-macro_rules! impl_eig_real {
- ($scalar:ty, $ev:path) => {
- impl Eig_ for $scalar {
- fn eig(
- calc_v: bool,
- l: MatrixLayout,
- mut a: &mut [Self],
- ) -> Result<(Vec, Vec)> {
+impl_eig_work_c!(c32, lapack_sys::cgeev_);
+impl_eig_work_c!(c64, lapack_sys::zgeev_);
+
+macro_rules! impl_eig_work_r {
+ ($f:ty, $ev:path) => {
+ impl EigWorkImpl for EigWork<$f> {
+ type Elem = $f;
+
+ fn new(calc_v: bool, l: MatrixLayout) -> Result {
let (n, _) = l.size();
- // Because LAPACK assumes F-continious array, C-continious array should be taken Hermitian conjugate.
- // However, we utilize a fact that left eigenvector of A^H corresponds to the right eigenvector of A
let (jobvl, jobvr) = if calc_v {
match l {
- MatrixLayout::C { .. } => (b'V', b'N'),
- MatrixLayout::F { .. } => (b'N', b'V'),
+ MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+ MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
}
} else {
- (b'N', b'N')
- };
- let mut eig_re = unsafe { vec_uninit(n as usize) };
- let mut eig_im = unsafe { vec_uninit(n as usize) };
-
- let mut vl = if jobvl == b'V' {
- Some(unsafe { vec_uninit((n * n) as usize) })
- } else {
- None
- };
- let mut vr = if jobvr == b'V' {
- Some(unsafe { vec_uninit((n * n) as usize) })
- } else {
- None
+ (JobEv::None, JobEv::None)
};
+ let mut eigs_re = vec_uninit(n as usize);
+ let mut eigs_im = vec_uninit(n as usize);
+ let mut vr_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let mut vr_r = jobvr.then(|| vec_uninit((n * n) as usize));
+ let vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
// calc work size
let mut info = 0;
- let mut work_size = [0.0];
+ let mut work_size: [$f; 1] = [0.0];
unsafe {
$ev(
- jobvl,
- jobvr,
- n,
- &mut a,
- n,
- &mut eig_re,
- &mut eig_im,
- vl.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- vr.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut work_size,
- -1,
+ jobvl.as_ptr(),
+ jobvr.as_ptr(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs_re),
+ AsPtr::as_mut_ptr(&mut eigs_im),
+ AsPtr::as_mut_ptr(vr_l.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(vr_r.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
&mut info,
)
};
@@ -167,92 +299,153 @@ macro_rules! impl_eig_real {
// actual ev
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
+ let work = vec_uninit(lwork);
+
+ Ok(Self {
+ n,
+ jobvr,
+ jobvl,
+ eigs: vec_uninit(n as usize),
+ eigs_re: Some(eigs_re),
+ eigs_im: Some(eigs_im),
+ rwork: None,
+ vr_l,
+ vr_r,
+ vc_l,
+ vc_r,
+ work,
+ })
+ }
+
+ fn calc<'work>(
+ &'work mut self,
+ a: &mut [Self::Elem],
+ ) -> Result> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
unsafe {
$ev(
- jobvl,
- jobvr,
- n,
- &mut a,
- n,
- &mut eig_re,
- &mut eig_im,
- vl.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- vr.as_mut().map(|v| v.as_mut_slice()).unwrap_or(&mut []),
- n,
- &mut work,
- lwork as i32,
+ self.jobvl.as_ptr(),
+ self.jobvr.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(self.eigs_re.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(self.eigs_im.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(self.vr_l.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(self.vr_r.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
&mut info,
)
};
info.as_lapack_result()?;
- // reconstruct eigenvalues
- let eigs: Vec = eig_re
- .iter()
- .zip(eig_im.iter())
- .map(|(&re, &im)| Self::complex(re, im))
- .collect();
+ let eigs_re = self
+ .eigs_re
+ .as_ref()
+ .map(|e| unsafe { e.slice_assume_init_ref() })
+ .unwrap();
+ let eigs_im = self
+ .eigs_im
+ .as_ref()
+ .map(|e| unsafe { e.slice_assume_init_ref() })
+ .unwrap();
+ reconstruct_eigs(eigs_re, eigs_im, &mut self.eigs);
- if !calc_v {
- return Ok((eigs, Vec::new()));
+ if let Some(v) = self.vr_l.as_ref() {
+ let v = unsafe { v.slice_assume_init_ref() };
+ reconstruct_eigenvectors(true, eigs_im, v, self.vc_l.as_mut().unwrap());
}
-
- // Reconstruct eigenvectors into complex-array
- // --------------------------------------------
- //
- // From LAPACK API https://software.intel.com/en-us/node/469230
- //
- // - If the j-th eigenvalue is real,
- // - v(j) = VR(:,j), the j-th column of VR.
- //
- // - If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
- // - v(j) = VR(:,j) + i*VR(:,j+1)
- // - v(j+1) = VR(:,j) - i*VR(:,j+1).
- //
- // ```
- // j -> <----pair----> <----pair---->
- // [ ... (real), (imag), (imag), (imag), (imag), ... ] : eigs
- // ^ ^ ^ ^ ^
- // false false true false true : is_conjugate_pair
- // ```
- let n = n as usize;
- let v = vr.or(vl).unwrap();
- let mut eigvecs = unsafe { vec_uninit(n * n) };
- let mut is_conjugate_pair = false; // flag for check `j` is complex conjugate
- for j in 0..n {
- if eig_im[j] == 0.0 {
- // j-th eigenvalue is real
- for i in 0..n {
- eigvecs[i + j * n] = Self::complex(v[i + j * n], 0.0);
- }
- } else {
- // j-th eigenvalue is complex
- // complex conjugated pair can be `j-1` or `j+1`
- if is_conjugate_pair {
- let j_pair = j - 1;
- assert!(j_pair < n);
- for i in 0..n {
- eigvecs[i + j * n] = Self::complex(v[i + j_pair * n], v[i + j * n]);
- }
- } else {
- let j_pair = j + 1;
- assert!(j_pair < n);
- for i in 0..n {
- eigvecs[i + j * n] =
- Self::complex(v[i + j * n], -v[i + j_pair * n]);
- }
- }
- is_conjugate_pair = !is_conjugate_pair;
- }
+ if let Some(v) = self.vr_r.as_ref() {
+ let v = unsafe { v.slice_assume_init_ref() };
+ reconstruct_eigenvectors(false, eigs_im, v, self.vc_r.as_mut().unwrap());
}
- Ok((eigs, eigvecs))
+ Ok(EigRef {
+ eigs: unsafe { self.eigs.slice_assume_init_ref() },
+ vl: self
+ .vc_l
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ vr: self
+ .vc_r
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ })
+ }
+
+ fn eval(mut self, a: &mut [Self::Elem]) -> Result> {
+ let _eig_ref = self.calc(a)?;
+ Ok(EigOwned {
+ eigs: unsafe { self.eigs.assume_init() },
+ vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+ vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+ })
}
}
};
}
+impl_eig_work_r!(f32, lapack_sys::sgeev_);
+impl_eig_work_r!(f64, lapack_sys::dgeev_);
-impl_eig_real!(f64, lapack::dgeev);
-impl_eig_real!(f32, lapack::sgeev);
+/// Reconstruct eigenvectors into complex-array
+///
+/// From LAPACK API https://software.intel.com/en-us/node/469230
+///
+/// - If the j-th eigenvalue is real,
+/// - v(j) = VR(:,j), the j-th column of VR.
+///
+/// - If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
+/// - v(j) = VR(:,j) + i*VR(:,j+1)
+/// - v(j+1) = VR(:,j) - i*VR(:,j+1).
+///
+/// In the C-layout case, we need the conjugates of the left
+/// eigenvectors, so the signs should be reversed.
+pub(crate) fn reconstruct_eigenvectors(
+ take_hermite_conjugate: bool,
+ eig_im: &[T],
+ vr: &[T],
+ vc: &mut [MaybeUninit],
+) {
+ let n = eig_im.len();
+ assert_eq!(vr.len(), n * n);
+ assert_eq!(vc.len(), n * n);
+
+ let mut col = 0;
+ while col < n {
+ if eig_im[col].is_zero() {
+ // The corresponding eigenvalue is real.
+ for row in 0..n {
+ let re = vr[row + col * n];
+ vc[row + col * n].write(T::complex(re, T::zero()));
+ }
+ col += 1;
+ } else {
+ // This is a complex conjugate pair.
+ assert!(col + 1 < n);
+ for row in 0..n {
+ let re = vr[row + col * n];
+ let mut im = vr[row + (col + 1) * n];
+ if take_hermite_conjugate {
+ im = -im;
+ }
+ vc[row + col * n].write(T::complex(re, im));
+ vc[row + (col + 1) * n].write(T::complex(re, -im));
+ }
+ col += 2;
+ }
+ }
+}
+
+/// Create complex eigenvalues from real and imaginary parts.
+fn reconstruct_eigs(re: &[T], im: &[T], eigs: &mut [MaybeUninit]) {
+ let n = eigs.len();
+ assert_eq!(re.len(), n);
+ assert_eq!(im.len(), n);
+ for i in 0..n {
+ eigs[i].write(T::complex(re[i], im[i]));
+ }
+}
diff --git a/lax/src/eig_generalized.rs b/lax/src/eig_generalized.rs
new file mode 100644
index 00000000..ea99dbdb
--- /dev/null
+++ b/lax/src/eig_generalized.rs
@@ -0,0 +1,520 @@
+//! Generalized eigenvalue problem for general matrices
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32 | f64 | c32 | c64 |
+//! |:------|:------|:------|:------|
+//! | sggev | dggev | cggev | zggev |
+//!
+use std::mem::MaybeUninit;
+
+use crate::eig::reconstruct_eigenvectors;
+use crate::{error::*, layout::MatrixLayout, *};
+use cauchy::*;
+use num_traits::{ToPrimitive, Zero};
+
+#[derive(Clone, PartialEq, Eq)]
+pub enum GeneralizedEigenvalue {
+ /// Finite generalized eigenvalue: `Finite(α/β, (α, β))`
+ Finite(T, (T, T)),
+
+ /// Indeterminate generalized eigenvalue: `Indeterminate((α, β))`
+ Indeterminate((T, T)),
+}
+
+impl std::fmt::Display for GeneralizedEigenvalue {
+ fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+ match self {
+ Self::Finite(e, (a, b)) => write!(f, "{e:.3e} ({a:.3e}/{b:.3e})"),
+ Self::Indeterminate((a, b)) => write!(f, "∞ ({a:.3e}/{b:.3e})"),
+ }
+ }
+}
+
+#[non_exhaustive]
+pub struct EigGeneralizedWork {
+ /// Problem size
+ pub n: i32,
+ /// Compute right eigenvectors or not
+ pub jobvr: JobEv,
+ /// Compute left eigenvectors or not
+ pub jobvl: JobEv,
+
+ /// Eigenvalues: alpha (numerators)
+ pub alpha: Vec>,
+ /// Eigenvalues: beta (denominators)
+ pub beta: Vec>,
+ /// Real part of alpha (eigenvalue numerators) used in real routines
+ pub alpha_re: Option>>,
+ /// Imaginary part of alpha (eigenvalue numerators) used in real routines
+ pub alpha_im: Option>>,
+ /// Real part of beta (eigenvalue denominators) used in real routines
+ pub beta_re: Option>>,
+ /// Imaginary part of beta (eigenvalue denominators) used in real routines
+ pub beta_im: Option>>,
+
+ /// Left eigenvectors
+ pub vc_l: Option>>,
+ /// Left eigenvectors used in real routines
+ pub vr_l: Option>>,
+ /// Right eigenvectors
+ pub vc_r: Option>>,
+ /// Right eigenvectors used in real routines
+ pub vr_r: Option>>,
+
+ /// Working memory
+ pub work: Vec>,
+ /// Working memory with `T::Real`
+ pub rwork: Option>>,
+}
+
+impl EigGeneralizedWork
+where
+ T: Scalar,
+ EigGeneralizedWork: EigGeneralizedWorkImpl,
+{
+ /// Create new working memory for eigenvalues compution.
+ pub fn new(calc_v: bool, l: MatrixLayout) -> Result {
+ EigGeneralizedWorkImpl::new(calc_v, l)
+ }
+
+ /// Compute eigenvalues and vectors on this working memory.
+ pub fn calc(&mut self, a: &mut [T], b: &mut [T]) -> Result> {
+ EigGeneralizedWorkImpl::calc(self, a, b)
+ }
+
+ /// Compute eigenvalues and vectors by consuming this working memory.
+ pub fn eval(self, a: &mut [T], b: &mut [T]) -> Result> {
+ EigGeneralizedWorkImpl::eval(self, a, b)
+ }
+}
+
+/// Owned result of eigenvalue problem by [EigGeneralizedWork::eval]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigGeneralizedOwned {
+ /// Eigenvalues
+ pub alpha: Vec,
+
+ pub beta: Vec,
+
+ /// Right eigenvectors
+ pub vr: Option>,
+
+ /// Left eigenvectors
+ pub vl: Option>,
+}
+
+/// Reference result of eigenvalue problem by [EigGeneralizedWork::calc]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigGeneralizedRef<'work, T: Scalar> {
+ /// Eigenvalues
+ pub alpha: &'work [T::Complex],
+
+ pub beta: &'work [T::Complex],
+
+ /// Right eigenvectors
+ pub vr: Option<&'work [T::Complex]>,
+
+ /// Left eigenvectors
+ pub vl: Option<&'work [T::Complex]>,
+}
+
+/// Helper trait for implementing [EigGeneralizedWork] methods
+pub trait EigGeneralizedWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(calc_v: bool, l: MatrixLayout) -> Result;
+ fn calc<'work>(
+ &'work mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result>;
+ fn eval(
+ self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result>;
+}
+
+macro_rules! impl_eig_generalized_work_c {
+ ($f:ty, $c:ty, $ev:path) => {
+ impl EigGeneralizedWorkImpl for EigGeneralizedWork<$c> {
+ type Elem = $c;
+
+ fn new(calc_v: bool, l: MatrixLayout) -> Result {
+ let (n, _) = l.size();
+ let (jobvl, jobvr) = if calc_v {
+ match l {
+ MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+ MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
+ }
+ } else {
+ (JobEv::None, JobEv::None)
+ };
+ let mut rwork = vec_uninit(8 * n as usize);
+
+ let mut alpha = vec_uninit(n as usize);
+ let mut beta = vec_uninit(n as usize);
+
+ let mut vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let mut vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
+
+ // calc work size
+ let mut info = 0;
+ let mut work_size = [<$c>::zero()];
+ unsafe {
+ $ev(
+ jobvl.as_ptr(),
+ jobvr.as_ptr(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut alpha),
+ AsPtr::as_mut_ptr(&mut beta),
+ AsPtr::as_mut_ptr(vc_l.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(vc_r.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ AsPtr::as_mut_ptr(&mut rwork),
+ &mut info,
+ )
+ };
+ info.as_lapack_result()?;
+
+ let lwork = work_size[0].to_usize().unwrap();
+ let work: Vec> = vec_uninit(lwork);
+ Ok(Self {
+ n,
+ jobvl,
+ jobvr,
+ alpha,
+ beta,
+ alpha_re: None,
+ alpha_im: None,
+ beta_re: None,
+ beta_im: None,
+ rwork: Some(rwork),
+ vc_l,
+ vc_r,
+ vr_l: None,
+ vr_r: None,
+ work,
+ })
+ }
+
+ fn calc<'work>(
+ &'work mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ unsafe {
+ $ev(
+ self.jobvl.as_ptr(),
+ self.jobvr.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(b),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.alpha),
+ AsPtr::as_mut_ptr(&mut self.beta),
+ AsPtr::as_mut_ptr(self.vc_l.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(self.vc_r.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
+ &mut info,
+ )
+ };
+ info.as_lapack_result()?;
+ // Hermite conjugate
+ if let Some(vl) = self.vc_l.as_mut() {
+ for value in vl {
+ let value = unsafe { value.assume_init_mut() };
+ value.im = -value.im;
+ }
+ }
+ Ok(EigGeneralizedRef {
+ alpha: unsafe { self.alpha.slice_assume_init_ref() },
+ beta: unsafe { self.beta.slice_assume_init_ref() },
+ vl: self
+ .vc_l
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ vr: self
+ .vc_r
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ })
+ }
+
+ fn eval(
+ mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let _eig_generalized_ref = self.calc(a, b)?;
+ Ok(EigGeneralizedOwned {
+ alpha: unsafe { self.alpha.assume_init() },
+ beta: unsafe { self.beta.assume_init() },
+ vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+ vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+ })
+ }
+ }
+
+ impl EigGeneralizedOwned<$c> {
+ pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec> {
+ self.alpha
+ .iter()
+ .zip(self.beta.iter())
+ .map(|(alpha, beta)| {
+ if let Some(thresh) = thresh_opt {
+ if beta.abs() < thresh {
+ GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+ } else {
+ GeneralizedEigenvalue::Finite(
+ alpha / beta,
+ (alpha.clone(), beta.clone()),
+ )
+ }
+ } else {
+ if beta.is_zero() {
+ GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+ } else {
+ GeneralizedEigenvalue::Finite(
+ alpha / beta,
+ (alpha.clone(), beta.clone()),
+ )
+ }
+ }
+ })
+ .collect::>()
+ }
+ }
+ };
+}
+
+impl_eig_generalized_work_c!(f32, c32, lapack_sys::cggev_);
+impl_eig_generalized_work_c!(f64, c64, lapack_sys::zggev_);
+
+macro_rules! impl_eig_generalized_work_r {
+ ($f:ty, $c:ty, $ev:path) => {
+ impl EigGeneralizedWorkImpl for EigGeneralizedWork<$f> {
+ type Elem = $f;
+
+ fn new(calc_v: bool, l: MatrixLayout) -> Result {
+ let (n, _) = l.size();
+ let (jobvl, jobvr) = if calc_v {
+ match l {
+ MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+ MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
+ }
+ } else {
+ (JobEv::None, JobEv::None)
+ };
+ let mut alpha_re = vec_uninit(n as usize);
+ let mut alpha_im = vec_uninit(n as usize);
+ let mut beta_re = vec_uninit(n as usize);
+ let mut vr_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let mut vr_r = jobvr.then(|| vec_uninit((n * n) as usize));
+ let vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+ let vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
+
+ // calc work size
+ let mut info = 0;
+ let mut work_size: [$f; 1] = [0.0];
+ unsafe {
+ $ev(
+ jobvl.as_ptr(),
+ jobvr.as_ptr(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut alpha_re),
+ AsPtr::as_mut_ptr(&mut alpha_im),
+ AsPtr::as_mut_ptr(&mut beta_re),
+ AsPtr::as_mut_ptr(vr_l.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(vr_r.as_deref_mut().unwrap_or(&mut [])),
+ &n,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ )
+ };
+ info.as_lapack_result()?;
+
+ // actual ev
+ let lwork = work_size[0].to_usize().unwrap();
+ let work = vec_uninit(lwork);
+
+ Ok(Self {
+ n,
+ jobvr,
+ jobvl,
+ alpha: vec_uninit(n as usize),
+ beta: vec_uninit(n as usize),
+ alpha_re: Some(alpha_re),
+ alpha_im: Some(alpha_im),
+ beta_re: Some(beta_re),
+ beta_im: None,
+ rwork: None,
+ vr_l,
+ vr_r,
+ vc_l,
+ vc_r,
+ work,
+ })
+ }
+
+ fn calc<'work>(
+ &'work mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ unsafe {
+ $ev(
+ self.jobvl.as_ptr(),
+ self.jobvr.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(b),
+ &self.n,
+ AsPtr::as_mut_ptr(self.alpha_re.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(self.alpha_im.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(self.beta_re.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(self.vr_l.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(self.vr_r.as_deref_mut().unwrap_or(&mut [])),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ )
+ };
+ info.as_lapack_result()?;
+
+ let alpha_re = self
+ .alpha_re
+ .as_ref()
+ .map(|e| unsafe { e.slice_assume_init_ref() })
+ .unwrap();
+ let alpha_im = self
+ .alpha_im
+ .as_ref()
+ .map(|e| unsafe { e.slice_assume_init_ref() })
+ .unwrap();
+ let beta_re = self
+ .beta_re
+ .as_ref()
+ .map(|e| unsafe { e.slice_assume_init_ref() })
+ .unwrap();
+ reconstruct_eigs_optional_im(alpha_re, Some(alpha_im), &mut self.alpha);
+ reconstruct_eigs_optional_im(beta_re, None, &mut self.beta);
+
+ if let Some(v) = self.vr_l.as_ref() {
+ let v = unsafe { v.slice_assume_init_ref() };
+ reconstruct_eigenvectors(true, alpha_im, v, self.vc_l.as_mut().unwrap());
+ }
+ if let Some(v) = self.vr_r.as_ref() {
+ let v = unsafe { v.slice_assume_init_ref() };
+ reconstruct_eigenvectors(false, alpha_im, v, self.vc_r.as_mut().unwrap());
+ }
+
+ Ok(EigGeneralizedRef {
+ alpha: unsafe { self.alpha.slice_assume_init_ref() },
+ beta: unsafe { self.beta.slice_assume_init_ref() },
+ vl: self
+ .vc_l
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ vr: self
+ .vc_r
+ .as_ref()
+ .map(|v| unsafe { v.slice_assume_init_ref() }),
+ })
+ }
+
+ fn eval(
+ mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let _eig_generalized_ref = self.calc(a, b)?;
+ Ok(EigGeneralizedOwned {
+ alpha: unsafe { self.alpha.assume_init() },
+ beta: unsafe { self.beta.assume_init() },
+ vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+ vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+ })
+ }
+ }
+
+ impl EigGeneralizedOwned<$f> {
+ pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec> {
+ self.alpha
+ .iter()
+ .zip(self.beta.iter())
+ .map(|(alpha, beta)| {
+ if let Some(thresh) = thresh_opt {
+ if beta.abs() < thresh {
+ GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+ } else {
+ GeneralizedEigenvalue::Finite(
+ alpha / beta,
+ (alpha.clone(), beta.clone()),
+ )
+ }
+ } else {
+ if beta.is_zero() {
+ GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+ } else {
+ GeneralizedEigenvalue::Finite(
+ alpha / beta,
+ (alpha.clone(), beta.clone()),
+ )
+ }
+ }
+ })
+ .collect::>()
+ }
+ }
+ };
+}
+impl_eig_generalized_work_r!(f32, c32, lapack_sys::sggev_);
+impl_eig_generalized_work_r!(f64, c64, lapack_sys::dggev_);
+
+/// Create complex eigenvalues from real and optional imaginary parts.
+fn reconstruct_eigs_optional_im(
+ re: &[T],
+ im_opt: Option<&[T]>,
+ eigs: &mut [MaybeUninit],
+) {
+ let n = eigs.len();
+ assert_eq!(re.len(), n);
+
+ if let Some(im) = im_opt {
+ assert_eq!(im.len(), n);
+ for i in 0..n {
+ eigs[i].write(T::complex(re[i], im[i]));
+ }
+ } else {
+ for i in 0..n {
+ eigs[i].write(T::complex(re[i], T::zero()));
+ }
+ }
+}
diff --git a/lax/src/eigh.rs b/lax/src/eigh.rs
index 46a3b131..bb3ca500 100644
--- a/lax/src/eigh.rs
+++ b/lax/src/eigh.rs
@@ -1,159 +1,190 @@
-//! Eigenvalue decomposition for Symmetric/Hermite matrices
+//! Eigenvalue problem for symmetric/Hermitian matricies
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32 | f64 | c32 | c64 |
+//! |:------|:------|:------|:------|
+//! | ssyev | dsyev | cheev | zheev |
use super::*;
use crate::{error::*, layout::MatrixLayout};
use cauchy::*;
use num_traits::{ToPrimitive, Zero};
-pub trait Eigh_: Scalar {
- /// Wraps `*syev` for real and `*heev` for complex
- fn eigh(
- calc_eigenvec: bool,
- layout: MatrixLayout,
- uplo: UPLO,
- a: &mut [Self],
- ) -> Result>;
+pub struct EighWork {
+ pub n: i32,
+ pub jobz: JobEv,
+ pub eigs: Vec>,
+ pub work: Vec>,
+ pub rwork: Option>>,
+}
- /// Wraps `*syegv` for real and `*heegv` for complex
- fn eigh_generalized(
- calc_eigenvec: bool,
- layout: MatrixLayout,
- uplo: UPLO,
- a: &mut [Self],
- b: &mut [Self],
- ) -> Result>;
+pub trait EighWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result;
+ fn calc(&mut self, uplo: UPLO, a: &mut [Self::Elem])
+ -> Result<&[::Real]>;
+ fn eval(self, uplo: UPLO, a: &mut [Self::Elem]) -> Result::Real>>;
}
-macro_rules! impl_eigh {
- (@real, $scalar:ty, $ev:path, $evg:path) => {
- impl_eigh!(@body, $scalar, $ev, $evg, );
- };
- (@complex, $scalar:ty, $ev:path, $evg:path) => {
- impl_eigh!(@body, $scalar, $ev, $evg, rwork);
- };
- (@body, $scalar:ty, $ev:path, $evg:path, $($rwork_ident:ident),*) => {
- impl Eigh_ for $scalar {
- fn eigh(
- calc_v: bool,
- layout: MatrixLayout,
- uplo: UPLO,
- mut a: &mut [Self],
- ) -> Result> {
+macro_rules! impl_eigh_work_c {
+ ($c:ty, $ev:path) => {
+ impl EighWorkImpl for EighWork<$c> {
+ type Elem = $c;
+
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result {
assert_eq!(layout.len(), layout.lda());
let n = layout.len();
- let jobz = if calc_v { b'V' } else { b'N' };
- let mut eigs = unsafe { vec_uninit(n as usize) };
-
- $(
- let mut $rwork_ident = unsafe { vec_uninit(3 * n as usize - 2 as usize) };
- )*
-
- // calc work size
+ let jobz = if calc_eigenvectors {
+ JobEv::All
+ } else {
+ JobEv::None
+ };
+ let mut eigs = vec_uninit(n as usize);
+ let mut rwork = vec_uninit(3 * n as usize - 2 as usize);
let mut info = 0;
- let mut work_size = [Self::zero()];
+ let mut work_size = [Self::Elem::zero()];
unsafe {
$ev(
- jobz,
- uplo as u8,
- n,
- &mut a,
- n,
- &mut eigs,
- &mut work_size,
- -1,
- $(&mut $rwork_ident,)*
+ jobz.as_ptr(),
+ UPLO::Upper.as_ptr(), // dummy, working memory is not affected by UPLO
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ AsPtr::as_mut_ptr(&mut rwork),
&mut info,
);
}
info.as_lapack_result()?;
-
- // actual ev
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
+ let work = vec_uninit(lwork);
+ Ok(EighWork {
+ n,
+ eigs,
+ jobz,
+ work,
+ rwork: Some(rwork),
+ })
+ }
+
+ fn calc(
+ &mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ ) -> Result<&[::Real]> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
unsafe {
$ev(
- jobz,
- uplo as u8,
- n,
- &mut a,
- n,
- &mut eigs,
- &mut work,
- lwork as i32,
- $(&mut $rwork_ident,)*
+ self.jobz.as_ptr(),
+ uplo.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.eigs),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
&mut info,
);
}
info.as_lapack_result()?;
- Ok(eigs)
+ Ok(unsafe { self.eigs.slice_assume_init_ref() })
}
- fn eigh_generalized(
- calc_v: bool,
- layout: MatrixLayout,
+ fn eval(
+ mut self,
uplo: UPLO,
- mut a: &mut [Self],
- mut b: &mut [Self],
- ) -> Result> {
- assert_eq!(layout.len(), layout.lda());
- let n = layout.len();
- let jobz = if calc_v { b'V' } else { b'N' };
- let mut eigs = unsafe { vec_uninit(n as usize) };
+ a: &mut [Self::Elem],
+ ) -> Result::Real>> {
+ let _eig = self.calc(uplo, a)?;
+ Ok(unsafe { self.eigs.assume_init() })
+ }
+ }
+ };
+}
+impl_eigh_work_c!(c64, lapack_sys::zheev_);
+impl_eigh_work_c!(c32, lapack_sys::cheev_);
- $(
- let mut $rwork_ident = unsafe { vec_uninit(3 * n as usize - 2) };
- )*
+macro_rules! impl_eigh_work_r {
+ ($f:ty, $ev:path) => {
+ impl EighWorkImpl for EighWork<$f> {
+ type Elem = $f;
- // calc work size
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result {
+ assert_eq!(layout.len(), layout.lda());
+ let n = layout.len();
+ let jobz = if calc_eigenvectors {
+ JobEv::All
+ } else {
+ JobEv::None
+ };
+ let mut eigs = vec_uninit(n as usize);
let mut info = 0;
- let mut work_size = [Self::zero()];
+ let mut work_size = [Self::Elem::zero()];
unsafe {
- $evg(
- &[1],
- jobz,
- uplo as u8,
- n,
- &mut a,
- n,
- &mut b,
- n,
- &mut eigs,
- &mut work_size,
- -1,
- $(&mut $rwork_ident,)*
+ $ev(
+ jobz.as_ptr(),
+ UPLO::Upper.as_ptr(), // dummy, working memory is not affected by UPLO
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
&mut info,
);
}
info.as_lapack_result()?;
-
- // actual evg
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
+ let work = vec_uninit(lwork);
+ Ok(EighWork {
+ n,
+ eigs,
+ jobz,
+ work,
+ rwork: None,
+ })
+ }
+
+ fn calc(
+ &mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ ) -> Result<&[::Real]> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
unsafe {
- $evg(
- &[1],
- jobz,
- uplo as u8,
- n,
- &mut a,
- n,
- &mut b,
- n,
- &mut eigs,
- &mut work,
- lwork as i32,
- $(&mut $rwork_ident,)*
+ $ev(
+ self.jobz.as_ptr(),
+ uplo.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.eigs),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
&mut info,
);
}
info.as_lapack_result()?;
- Ok(eigs)
+ Ok(unsafe { self.eigs.slice_assume_init_ref() })
+ }
+
+ fn eval(
+ mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ ) -> Result::Real>> {
+ let _eig = self.calc(uplo, a)?;
+ Ok(unsafe { self.eigs.assume_init() })
}
}
};
-} // impl_eigh!
-
-impl_eigh!(@real, f64, lapack::dsyev, lapack::dsygv);
-impl_eigh!(@real, f32, lapack::ssyev, lapack::ssygv);
-impl_eigh!(@complex, c64, lapack::zheev, lapack::zhegv);
-impl_eigh!(@complex, c32, lapack::cheev, lapack::chegv);
+}
+impl_eigh_work_r!(f64, lapack_sys::dsyev_);
+impl_eigh_work_r!(f32, lapack_sys::ssyev_);
diff --git a/lax/src/eigh_generalized.rs b/lax/src/eigh_generalized.rs
new file mode 100644
index 00000000..5d4d83ca
--- /dev/null
+++ b/lax/src/eigh_generalized.rs
@@ -0,0 +1,216 @@
+//! Generalized eigenvalue problem for symmetric/Hermitian matrices
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32 | f64 | c32 | c64 |
+//! |:------|:------|:------|:------|
+//! | ssygv | dsygv | chegv | zhegv |
+//!
+
+use super::*;
+use crate::{error::*, layout::MatrixLayout};
+use cauchy::*;
+use num_traits::{ToPrimitive, Zero};
+
+pub struct EighGeneralizedWork {
+ pub n: i32,
+ pub jobz: JobEv,
+ pub eigs: Vec>,
+ pub work: Vec>,
+ pub rwork: Option>>,
+}
+
+pub trait EighGeneralizedWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result;
+ fn calc(
+ &mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result<&[::Real]>;
+ fn eval(
+ self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result::Real>>;
+}
+
+macro_rules! impl_eigh_generalized_work_c {
+ ($c:ty, $gv:path) => {
+ impl EighGeneralizedWorkImpl for EighGeneralizedWork<$c> {
+ type Elem = $c;
+
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result {
+ assert_eq!(layout.len(), layout.lda());
+ let n = layout.len();
+ let jobz = if calc_eigenvectors {
+ JobEv::All
+ } else {
+ JobEv::None
+ };
+ let mut eigs = vec_uninit(n as usize);
+ let mut rwork = vec_uninit(3 * n as usize - 2 as usize);
+ let mut info = 0;
+ let mut work_size = [Self::Elem::zero()];
+ unsafe {
+ $gv(
+ &1, // ITYPE A*x = (lambda)*B*x
+ jobz.as_ptr(),
+ UPLO::Upper.as_ptr(), // dummy, working memory is not affected by UPLO
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ AsPtr::as_mut_ptr(&mut rwork),
+ &mut info,
+ );
+ }
+ info.as_lapack_result()?;
+ let lwork = work_size[0].to_usize().unwrap();
+ let work = vec_uninit(lwork);
+ Ok(EighGeneralizedWork {
+ n,
+ eigs,
+ jobz,
+ work,
+ rwork: Some(rwork),
+ })
+ }
+
+ fn calc(
+ &mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result<&[::Real]> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ unsafe {
+ $gv(
+ &1, // ITYPE A*x = (lambda)*B*x
+ self.jobz.as_ptr(),
+ uplo.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(b),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.eigs),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
+ &mut info,
+ );
+ }
+ info.as_lapack_result()?;
+ Ok(unsafe { self.eigs.slice_assume_init_ref() })
+ }
+
+ fn eval(
+ mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result::Real>> {
+ let _eig = self.calc(uplo, a, b)?;
+ Ok(unsafe { self.eigs.assume_init() })
+ }
+ }
+ };
+}
+impl_eigh_generalized_work_c!(c64, lapack_sys::zhegv_);
+impl_eigh_generalized_work_c!(c32, lapack_sys::chegv_);
+
+macro_rules! impl_eigh_generalized_work_r {
+ ($f:ty, $gv:path) => {
+ impl EighGeneralizedWorkImpl for EighGeneralizedWork<$f> {
+ type Elem = $f;
+
+ fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result {
+ assert_eq!(layout.len(), layout.lda());
+ let n = layout.len();
+ let jobz = if calc_eigenvectors {
+ JobEv::All
+ } else {
+ JobEv::None
+ };
+ let mut eigs = vec_uninit(n as usize);
+ let mut info = 0;
+ let mut work_size = [Self::Elem::zero()];
+ unsafe {
+ $gv(
+ &1, // ITYPE A*x = (lambda)*B*x
+ jobz.as_ptr(),
+ UPLO::Upper.as_ptr(), // dummy, working memory is not affected by UPLO
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ std::ptr::null_mut(),
+ &n,
+ AsPtr::as_mut_ptr(&mut eigs),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ );
+ }
+ info.as_lapack_result()?;
+ let lwork = work_size[0].to_usize().unwrap();
+ let work = vec_uninit(lwork);
+ Ok(EighGeneralizedWork {
+ n,
+ eigs,
+ jobz,
+ work,
+ rwork: None,
+ })
+ }
+
+ fn calc(
+ &mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result<&[::Real]> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ unsafe {
+ $gv(
+ &1, // ITYPE A*x = (lambda)*B*x
+ self.jobz.as_ptr(),
+ uplo.as_ptr(),
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.n,
+ AsPtr::as_mut_ptr(b),
+ &self.n,
+ AsPtr::as_mut_ptr(&mut self.eigs),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ );
+ }
+ info.as_lapack_result()?;
+ Ok(unsafe { self.eigs.slice_assume_init_ref() })
+ }
+
+ fn eval(
+ mut self,
+ uplo: UPLO,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result::Real>> {
+ let _eig = self.calc(uplo, a, b)?;
+ Ok(unsafe { self.eigs.assume_init() })
+ }
+ }
+ };
+}
+impl_eigh_generalized_work_r!(f64, lapack_sys::dsygv_);
+impl_eigh_generalized_work_r!(f32, lapack_sys::ssygv_);
diff --git a/lax/src/error.rs b/lax/src/error.rs
index fb4b9838..e1c314ee 100644
--- a/lax/src/error.rs
+++ b/lax/src/error.rs
@@ -11,7 +11,7 @@ pub enum Error {
LapackInvalidValue { return_code: i32 },
#[error(
- "Comutational failure in LAPACK subroutine: return_code = {}",
+ "Computational failure in LAPACK subroutine: return_code = {}",
return_code
)]
LapackComputationalFailure { return_code: i32 },
diff --git a/lax/src/flags.rs b/lax/src/flags.rs
new file mode 100644
index 00000000..f9dea20d
--- /dev/null
+++ b/lax/src/flags.rs
@@ -0,0 +1,140 @@
+//! Charactor flags, e.g. `'T'`, used in LAPACK API
+use core::ffi::c_char;
+
+/// Upper/Lower specification for seveal usages
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum UPLO {
+ Upper = b'U',
+ Lower = b'L',
+}
+
+impl UPLO {
+ pub fn t(self) -> Self {
+ match self {
+ UPLO::Upper => UPLO::Lower,
+ UPLO::Lower => UPLO::Upper,
+ }
+ }
+
+ /// To use Fortran LAPACK API in lapack-sys crate
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const UPLO as *const c_char
+ }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum Transpose {
+ No = b'N',
+ Transpose = b'T',
+ Hermite = b'C',
+}
+
+impl Transpose {
+ /// To use Fortran LAPACK API in lapack-sys crate
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const Transpose as *const c_char
+ }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum NormType {
+ One = b'O',
+ Infinity = b'I',
+ Frobenius = b'F',
+}
+
+impl NormType {
+ pub fn transpose(self) -> Self {
+ match self {
+ NormType::One => NormType::Infinity,
+ NormType::Infinity => NormType::One,
+ NormType::Frobenius => NormType::Frobenius,
+ }
+ }
+
+ /// To use Fortran LAPACK API in lapack-sys crate
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const NormType as *const c_char
+ }
+}
+
+/// Flag for calculating eigenvectors or not
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum JobEv {
+ /// Calculate eigenvectors in addition to eigenvalues
+ All = b'V',
+ /// Do not calculate eigenvectors. Only calculate eigenvalues.
+ None = b'N',
+}
+
+impl JobEv {
+ pub fn is_calc(&self) -> bool {
+ match self {
+ JobEv::All => true,
+ JobEv::None => false,
+ }
+ }
+
+ pub fn then T>(&self, f: F) -> Option {
+ if self.is_calc() {
+ Some(f())
+ } else {
+ None
+ }
+ }
+
+ /// To use Fortran LAPACK API in lapack-sys crate
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const JobEv as *const c_char
+ }
+}
+
+/// Specifies how many singular vectors are computed
+///
+/// For an input matrix $A$ of shape $m \times n$,
+/// the following are computed on the singular value decomposition $A = U\Sigma V^T$:
+#[cfg_attr(doc, katexit::katexit)]
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum JobSvd {
+ /// All $m$ columns of $U$, and/or all $n$ rows of $V^T$.
+ All = b'A',
+ /// The first $\min(m, n)$ columns of $U$ and/or the first $\min(m, n)$ rows of $V^T$.
+ Some = b'S',
+ /// No columns of $U$ and/or rows of $V^T$.
+ None = b'N',
+}
+
+impl JobSvd {
+ pub fn from_bool(calc_uv: bool) -> Self {
+ if calc_uv {
+ JobSvd::All
+ } else {
+ JobSvd::None
+ }
+ }
+
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const JobSvd as *const c_char
+ }
+}
+
+/// Specify whether input triangular matrix is unit or not
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[repr(u8)]
+pub enum Diag {
+ /// Unit triangular matrix, i.e. all diagonal elements of the matrix are `1`
+ Unit = b'U',
+ /// Non-unit triangular matrix. Its diagonal elements may be different from `1`
+ NonUnit = b'N',
+}
+
+impl Diag {
+ pub fn as_ptr(&self) -> *const c_char {
+ self as *const Diag as *const c_char
+ }
+}
diff --git a/lax/src/layout.rs b/lax/src/layout.rs
index e7ab1da4..28b35122 100644
--- a/lax/src/layout.rs
+++ b/lax/src/layout.rs
@@ -37,7 +37,7 @@
//! This `S` for a matrix `A` is called "leading dimension of the array A" in LAPACK document, and denoted by `lda`.
//!
-use cauchy::Scalar;
+use super::*;
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum MatrixLayout {
@@ -153,7 +153,7 @@ impl MatrixLayout {
/// ------
/// - If size of `a` and `layout` size mismatch
///
-pub fn square_transpose(layout: MatrixLayout, a: &mut [T]) {
+pub fn square_transpose(layout: MatrixLayout, a: &mut [T]) {
let (m, n) = layout.size();
let n = n as usize;
let m = m as usize;
@@ -162,23 +162,78 @@ pub fn square_transpose(layout: MatrixLayout, a: &mut [T]) {
for j in (i + 1)..n {
let a_ij = a[i * n + j];
let a_ji = a[j * m + i];
- a[i * n + j] = a_ji.conj();
- a[j * m + i] = a_ij.conj();
+ a[i * n + j] = a_ji;
+ a[j * m + i] = a_ij;
}
}
}
/// Out-place transpose for general matrix
///
-/// Inplace transpose of non-square matrices is hard.
-/// See also: https://en.wikipedia.org/wiki/In-place_matrix_transposition
+/// Examples
+/// ---------
+///
+/// ```rust
+/// # use lax::layout::*;
+/// let layout = MatrixLayout::C { row: 2, lda: 3 };
+/// let a = vec![1., 2., 3., 4., 5., 6.];
+/// let (l, b) = transpose(layout, &a);
+/// assert_eq!(l, MatrixLayout::F { col: 3, lda: 2 });
+/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
+/// ```
+///
+/// ```rust
+/// # use lax::layout::*;
+/// let layout = MatrixLayout::F { col: 2, lda: 3 };
+/// let a = vec![1., 2., 3., 4., 5., 6.];
+/// let (l, b) = transpose(layout, &a);
+/// assert_eq!(l, MatrixLayout::C { row: 3, lda: 2 });
+/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
+/// ```
+///
+/// Panics
+/// ------
+/// - If input array size and `layout` size mismatch
+///
+pub fn transpose(layout: MatrixLayout, input: &[T]) -> (MatrixLayout, Vec) {
+ let (m, n) = layout.size();
+ let transposed = layout.resized(n, m).t();
+ let m = m as usize;
+ let n = n as usize;
+ assert_eq!(input.len(), m * n);
+
+ let mut out: Vec> = vec_uninit(m * n);
+
+ match layout {
+ MatrixLayout::C { .. } => {
+ for i in 0..m {
+ for j in 0..n {
+ out[j * m + i].write(input[i * n + j]);
+ }
+ }
+ }
+ MatrixLayout::F { .. } => {
+ for i in 0..m {
+ for j in 0..n {
+ out[i * n + j].write(input[j * m + i]);
+ }
+ }
+ }
+ }
+ (transposed, unsafe { out.assume_init() })
+}
+
+/// Out-place transpose for general matrix
+///
+/// Examples
+/// ---------
///
/// ```rust
/// # use lax::layout::*;
/// let layout = MatrixLayout::C { row: 2, lda: 3 };
/// let a = vec![1., 2., 3., 4., 5., 6.];
/// let mut b = vec![0.0; a.len()];
-/// let l = transpose(layout, &a, &mut b);
+/// let l = transpose_over(layout, &a, &mut b);
/// assert_eq!(l, MatrixLayout::F { col: 3, lda: 2 });
/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
/// ```
@@ -188,16 +243,16 @@ pub fn square_transpose(layout: MatrixLayout, a: &mut [T]) {
/// let layout = MatrixLayout::F { col: 2, lda: 3 };
/// let a = vec![1., 2., 3., 4., 5., 6.];
/// let mut b = vec![0.0; a.len()];
-/// let l = transpose(layout, &a, &mut b);
+/// let l = transpose_over(layout, &a, &mut b);
/// assert_eq!(l, MatrixLayout::C { row: 3, lda: 2 });
/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
/// ```
///
/// Panics
/// ------
-/// - If size of `a` and `layout` size mismatch
+/// - If input array sizes and `layout` size mismatch
///
-pub fn transpose(layout: MatrixLayout, from: &[T], to: &mut [T]) -> MatrixLayout {
+pub fn transpose_over(layout: MatrixLayout, from: &[T], to: &mut [T]) -> MatrixLayout {
let (m, n) = layout.size();
let transposed = layout.resized(n, m).t();
let m = m as usize;
diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs
index fc378aa6..d0bb7def 100644
--- a/lax/src/least_squares.rs
+++ b/lax/src/least_squares.rs
@@ -5,154 +5,328 @@ use cauchy::*;
use num_traits::{ToPrimitive, Zero};
/// Result of LeastSquares
-pub struct LeastSquaresOutput {
+pub struct LeastSquaresOwned {
/// singular values
pub singular_values: Vec,
/// The rank of the input matrix A
pub rank: i32,
}
-/// Wraps `*gelsd`
-pub trait LeastSquaresSvdDivideConquer_: Scalar {
- fn least_squares(
- a_layout: MatrixLayout,
- a: &mut [Self],
- b: &mut [Self],
- ) -> Result>;
-
- fn least_squares_nrhs(
- a_layout: MatrixLayout,
- a: &mut [Self],
- b_layout: MatrixLayout,
- b: &mut [Self],
- ) -> Result>;
+/// Result of LeastSquares
+pub struct LeastSquaresRef<'work, A: Scalar> {
+ /// singular values
+ pub singular_values: &'work [A::Real],
+ /// The rank of the input matrix A
+ pub rank: i32,
}
-macro_rules! impl_least_squares {
- (@real, $scalar:ty, $gelsd:path) => {
- impl_least_squares!(@body, $scalar, $gelsd, );
- };
- (@complex, $scalar:ty, $gelsd:path) => {
- impl_least_squares!(@body, $scalar, $gelsd, rwork);
- };
+pub struct LeastSquaresWork {
+ pub a_layout: MatrixLayout,
+ pub b_layout: MatrixLayout,
+ pub singular_values: Vec>,
+ pub work: Vec>,
+ pub iwork: Vec>,
+ pub rwork: Option>>,
+}
- (@body, $scalar:ty, $gelsd:path, $($rwork:ident),*) => {
- impl LeastSquaresSvdDivideConquer_ for $scalar {
- fn least_squares(
- l: MatrixLayout,
- a: &mut [Self],
- b: &mut [Self],
- ) -> Result> {
- let b_layout = l.resized(b.len() as i32, 1);
- Self::least_squares_nrhs(l, a, b_layout, b)
- }
+pub trait LeastSquaresWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(a_layout: MatrixLayout, b_layout: MatrixLayout) -> Result;
+ fn calc(
+ &mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result>;
+ fn eval(
+ self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result>;
+}
+
+macro_rules! impl_least_squares_work_c {
+ ($c:ty, $lsd:path) => {
+ impl LeastSquaresWorkImpl for LeastSquaresWork<$c> {
+ type Elem = $c;
- fn least_squares_nrhs(
- a_layout: MatrixLayout,
- a: &mut [Self],
- b_layout: MatrixLayout,
- b: &mut [Self],
- ) -> Result> {
- // Minimize |b - Ax|_2
- //
- // where
- // A : (m, n)
- // b : (max(m, n), nrhs) // `b` has to store `x` on exit
- // x : (n, nrhs)
+ fn new(a_layout: MatrixLayout, b_layout: MatrixLayout) -> Result {
let (m, n) = a_layout.size();
let (m_, nrhs) = b_layout.size();
let k = m.min(n);
assert!(m_ >= m);
+ let rcond = -1.;
+ let mut singular_values = vec_uninit(k as usize);
+ let mut rank: i32 = 0;
+
+ // eval work size
+ let mut info = 0;
+ let mut work_size = [Self::Elem::zero()];
+ let mut iwork_size = [0];
+ let mut rwork = [::Real::zero()];
+ unsafe {
+ $lsd(
+ &m,
+ &n,
+ &nrhs,
+ std::ptr::null_mut(),
+ &m,
+ std::ptr::null_mut(),
+ &m_,
+ AsPtr::as_mut_ptr(&mut singular_values),
+ &rcond,
+ &mut rank,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ AsPtr::as_mut_ptr(&mut rwork),
+ iwork_size.as_mut_ptr(),
+ &mut info,
+ )
+ };
+ info.as_lapack_result()?;
+
+ let lwork = work_size[0].to_usize().unwrap();
+ let liwork = iwork_size[0].to_usize().unwrap();
+ let lrwork = rwork[0].to_usize().unwrap();
+
+ let work = vec_uninit(lwork);
+ let iwork = vec_uninit(liwork);
+ let rwork = vec_uninit(lrwork);
+
+ Ok(LeastSquaresWork {
+ a_layout,
+ b_layout,
+ work,
+ iwork,
+ rwork: Some(rwork),
+ singular_values,
+ })
+ }
+
+ fn calc(
+ &mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let (m, n) = self.a_layout.size();
+ let (m_, nrhs) = self.b_layout.size();
+ assert!(m_ >= m);
+
+ let lwork = self.work.len().to_i32().unwrap();
+
// Transpose if a is C-continuous
let mut a_t = None;
- let a_layout = match a_layout {
+ let _ = match self.a_layout {
MatrixLayout::C { .. } => {
- a_t = Some(unsafe { vec_uninit( a.len()) });
- transpose(a_layout, a, a_t.as_mut().unwrap())
+ let (layout, t) = transpose(self.a_layout, a);
+ a_t = Some(t);
+ layout
}
- MatrixLayout::F { .. } => a_layout,
+ MatrixLayout::F { .. } => self.a_layout,
};
// Transpose if b is C-continuous
let mut b_t = None;
- let b_layout = match b_layout {
+ let b_layout = match self.b_layout {
MatrixLayout::C { .. } => {
- b_t = Some(unsafe { vec_uninit( b.len()) });
- transpose(b_layout, b, b_t.as_mut().unwrap())
+ let (layout, t) = transpose(self.b_layout, b);
+ b_t = Some(t);
+ layout
}
- MatrixLayout::F { .. } => b_layout,
+ MatrixLayout::F { .. } => self.b_layout,
};
- let rcond: Self::Real = -1.;
- let mut singular_values: Vec = unsafe { vec_uninit( k as usize) };
+ let rcond: ::Real = -1.;
+ let mut rank: i32 = 0;
+
+ let mut info = 0;
+ unsafe {
+ $lsd(
+ &m,
+ &n,
+ &nrhs,
+ AsPtr::as_mut_ptr(a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a)),
+ &m,
+ AsPtr::as_mut_ptr(b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b)),
+ &m_,
+ AsPtr::as_mut_ptr(&mut self.singular_values),
+ &rcond,
+ &mut rank,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
+ AsPtr::as_mut_ptr(&mut self.iwork),
+ &mut info,
+ );
+ }
+ info.as_lapack_result()?;
+
+ let singular_values = unsafe { self.singular_values.slice_assume_init_ref() };
+
+ // Skip a_t -> a transpose because A has been destroyed
+ // Re-transpose b
+ if let Some(b_t) = b_t {
+ transpose_over(b_layout, &b_t, b);
+ }
+
+ Ok(LeastSquaresRef {
+ singular_values,
+ rank,
+ })
+ }
+
+ fn eval(
+ mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let LeastSquaresRef { rank, .. } = self.calc(a, b)?;
+ let singular_values = unsafe { self.singular_values.assume_init() };
+ Ok(LeastSquaresOwned {
+ singular_values,
+ rank,
+ })
+ }
+ }
+ };
+}
+impl_least_squares_work_c!(c64, lapack_sys::zgelsd_);
+impl_least_squares_work_c!(c32, lapack_sys::cgelsd_);
+
+macro_rules! impl_least_squares_work_r {
+ ($c:ty, $lsd:path) => {
+ impl LeastSquaresWorkImpl for LeastSquaresWork<$c> {
+ type Elem = $c;
+
+ fn new(a_layout: MatrixLayout, b_layout: MatrixLayout) -> Result {
+ let (m, n) = a_layout.size();
+ let (m_, nrhs) = b_layout.size();
+ let k = m.min(n);
+ assert!(m_ >= m);
+
+ let rcond = -1.;
+ let mut singular_values = vec_uninit(k as usize);
let mut rank: i32 = 0;
// eval work size
let mut info = 0;
- let mut work_size = [Self::zero()];
+ let mut work_size = [Self::Elem::zero()];
let mut iwork_size = [0];
- $(
- let mut $rwork = [Self::Real::zero()];
- )*
unsafe {
- $gelsd(
- m,
- n,
- nrhs,
- a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a),
- a_layout.lda(),
- b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b),
- b_layout.lda(),
- &mut singular_values,
- rcond,
+ $lsd(
+ &m,
+ &n,
+ &nrhs,
+ std::ptr::null_mut(),
+ &m,
+ std::ptr::null_mut(),
+ &m_,
+ AsPtr::as_mut_ptr(&mut singular_values),
+ &rcond,
&mut rank,
- &mut work_size,
- -1,
- $(&mut $rwork,)*
- &mut iwork_size,
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ iwork_size.as_mut_ptr(),
&mut info,
)
};
info.as_lapack_result()?;
- // calc
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit( lwork) };
let liwork = iwork_size[0].to_usize().unwrap();
- let mut iwork = unsafe { vec_uninit( liwork) };
- $(
- let lrwork = $rwork[0].to_usize().unwrap();
- let mut $rwork = unsafe { vec_uninit( lrwork) };
- )*
+
+ let work = vec_uninit(lwork);
+ let iwork = vec_uninit(liwork);
+
+ Ok(LeastSquaresWork {
+ a_layout,
+ b_layout,
+ work,
+ iwork,
+ rwork: None,
+ singular_values,
+ })
+ }
+
+ fn calc(
+ &mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let (m, n) = self.a_layout.size();
+ let (m_, nrhs) = self.b_layout.size();
+ assert!(m_ >= m);
+
+ let lwork = self.work.len().to_i32().unwrap();
+
+ // Transpose if a is C-continuous
+ let mut a_t = None;
+ let _ = match self.a_layout {
+ MatrixLayout::C { .. } => {
+ let (layout, t) = transpose(self.a_layout, a);
+ a_t = Some(t);
+ layout
+ }
+ MatrixLayout::F { .. } => self.a_layout,
+ };
+
+ // Transpose if b is C-continuous
+ let mut b_t = None;
+ let b_layout = match self.b_layout {
+ MatrixLayout::C { .. } => {
+ let (layout, t) = transpose(self.b_layout, b);
+ b_t = Some(t);
+ layout
+ }
+ MatrixLayout::F { .. } => self.b_layout,
+ };
+
+ let rcond: ::Real = -1.;
+ let mut rank: i32 = 0;
+
+ let mut info = 0;
unsafe {
- $gelsd(
- m,
- n,
- nrhs,
- a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a),
- a_layout.lda(),
- b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b),
- b_layout.lda(),
- &mut singular_values,
- rcond,
+ $lsd(
+ &m,
+ &n,
+ &nrhs,
+ AsPtr::as_mut_ptr(a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a)),
+ &m,
+ AsPtr::as_mut_ptr(b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b)),
+ &m_,
+ AsPtr::as_mut_ptr(&mut self.singular_values),
+ &rcond,
&mut rank,
- &mut work,
- lwork as i32,
- $(&mut $rwork,)*
- &mut iwork,
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ AsPtr::as_mut_ptr(&mut self.iwork),
&mut info,
);
}
info.as_lapack_result()?;
+ let singular_values = unsafe { self.singular_values.slice_assume_init_ref() };
+
// Skip a_t -> a transpose because A has been destroyed
// Re-transpose b
if let Some(b_t) = b_t {
- transpose(b_layout, &b_t, b);
+ transpose_over(b_layout, &b_t, b);
}
- Ok(LeastSquaresOutput {
+ Ok(LeastSquaresRef {
+ singular_values,
+ rank,
+ })
+ }
+
+ fn eval(
+ mut self,
+ a: &mut [Self::Elem],
+ b: &mut [Self::Elem],
+ ) -> Result> {
+ let LeastSquaresRef { rank, .. } = self.calc(a, b)?;
+ let singular_values = unsafe { self.singular_values.assume_init() };
+ Ok(LeastSquaresOwned {
singular_values,
rank,
})
@@ -160,8 +334,5 @@ macro_rules! impl_least_squares {
}
};
}
-
-impl_least_squares!(@real, f64, lapack::dgelsd);
-impl_least_squares!(@real, f32, lapack::sgelsd);
-impl_least_squares!(@complex, c64, lapack::zgelsd);
-impl_least_squares!(@complex, c32, lapack::cgelsd);
+impl_least_squares_work_r!(f64, lapack_sys::dgelsd_);
+impl_least_squares_work_r!(f32, lapack_sys::sgelsd_);
diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 41c15237..680ff0db 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -1,63 +1,80 @@
-//! Linear Algebra eXtension (LAX)
-//! ===============================
+//! Safe Rust wrapper for LAPACK without external dependency.
//!
-//! ndarray-free safe Rust wrapper for LAPACK FFI
+//! [Lapack] trait
+//! ----------------
//!
-//! Linear equation, Inverse matrix, Condition number
-//! --------------------------------------------------
+//! This crates provides LAPACK wrapper as a traits.
+//! For example, LU decomposition of general matrices is provided like:
//!
-//! As the property of $A$, several types of triangular factorization are used:
+//! ```ignore
+//! pub trait Lapack {
+//! fn lu(l: MatrixLayout, a: &mut [Self]) -> Result;
+//! }
+//! ```
//!
-//! - LU-decomposition for general matrix
-//! - $PA = LU$, where $L$ is lower matrix, $U$ is upper matrix, and $P$ is permutation matrix
-//! - Bunch-Kaufman diagonal pivoting method for nonpositive-definite Hermitian matrix
-//! - $A = U D U^\dagger$, where $U$ is upper matrix,
-//! $D$ is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
+//! see [Lapack] for detail.
+//! This trait is implemented for [f32], [f64], [c32] which is an alias to `num::Complex`,
+//! and [c64] which is an alias to `num::Complex`.
+//! You can use it like `f64::lu`:
//!
-//! | matrix type | Triangler factorization (TRF) | Solve (TRS) | Inverse matrix (TRI) | Reciprocal condition number (CON) |
-//! |:--------------------------------|:------------------------------|:------------|:---------------------|:----------------------------------|
-//! | General (GE) | [lu] | [solve] | [inv] | [rcond] |
-//! | Symmetric (SY) / Hermitian (HE) | [bk] | [solveh] | [invh] | - |
+//! ```
+//! use lax::{Lapack, layout::MatrixLayout, Transpose};
//!
-//! [lu]: solve/trait.Solve_.html#tymethod.lu
-//! [solve]: solve/trait.Solve_.html#tymethod.solve
-//! [inv]: solve/trait.Solve_.html#tymethod.inv
-//! [rcond]: solve/trait.Solve_.html#tymethod.rcond
+//! let mut a = vec![
+//! 1.0, 2.0,
+//! 3.0, 4.0
+//! ];
+//! let mut b = vec![1.0, 2.0];
+//! let layout = MatrixLayout::C { row: 2, lda: 2 };
+//! let pivot = f64::lu(layout, &mut a).unwrap();
+//! f64::solve(layout, Transpose::No, &a, &pivot, &mut b).unwrap();
+//! ```
//!
-//! [bk]: solveh/trait.Solveh_.html#tymethod.bk
-//! [solveh]: solveh/trait.Solveh_.html#tymethod.solveh
-//! [invh]: solveh/trait.Solveh_.html#tymethod.invh
+//! When you want to write generic algorithm for real and complex matrices,
+//! this trait can be used as a trait bound:
//!
-//! Eigenvalue Problem
-//! -------------------
+//! ```
+//! use lax::{Lapack, layout::MatrixLayout, Transpose};
//!
-//! Solve eigenvalue problem for a matrix $A$
+//! fn solve_at_once(layout: MatrixLayout, a: &mut [T], b: &mut [T]) -> Result<(), lax::error::Error> {
+//! let pivot = T::lu(layout, a)?;
+//! T::solve(layout, Transpose::No, a, &pivot, b)?;
+//! Ok(())
+//! }
+//! ```
//!
-//! $$ Av_i = \lambda_i v_i $$
+//! There are several similar traits as described below to keep development easy.
+//! They are merged into a single trait, [Lapack].
//!
-//! or generalized eigenvalue problem
+//! Linear equation, Inverse matrix, Condition number
+//! --------------------------------------------------
+//!
+//! According to the property input metrix, several types of triangular decomposition are used:
//!
-//! $$ Av_i = \lambda_i B v_i $$
+//! - [solve] module provides methods for LU-decomposition for general matrix.
+//! - [solveh] module provides methods for Bunch-Kaufman diagonal pivoting method for symmetric/Hermitian indefinite matrix.
+//! - [cholesky] module provides methods for Cholesky decomposition for symmetric/Hermitian positive dinite matrix.
//!
-//! | matrix type | Eigenvalue (EV) | Generalized Eigenvalue Problem (EG) |
-//! |:--------------------------------|:----------------|:------------------------------------|
-//! | General (GE) |[eig] | - |
-//! | Symmetric (SY) / Hermitian (HE) |[eigh] |[eigh_generalized] |
+//! Eigenvalue Problem
+//! -------------------
//!
-//! [eig]: eig/trait.Eig_.html#tymethod.eig
-//! [eigh]: eigh/trait.Eigh_.html#tymethod.eigh
-//! [eigh_generalized]: eigh/trait.Eigh_.html#tymethod.eigh_generalized
+//! According to the property input metrix,
+//! there are several types of eigenvalue problem API
//!
-//! Singular Value Decomposition (SVD), Least square problem
-//! ----------------------------------------------------------
+//! - [eig] module for eigenvalue problem for general matrix.
+//! - [eig_generalized] module for generalized eigenvalue problem for general matrix.
+//! - [eigh] module for eigenvalue problem for symmetric/Hermitian matrix.
+//! - [eigh_generalized] module for generalized eigenvalue problem for symmetric/Hermitian matrix.
//!
-//! | matrix type | Singular Value Decomposition (SVD) | SVD with divided-and-conquer (SDD) | Least square problem (LSD) |
-//! |:-------------|:-----------------------------------|:-----------------------------------|:---------------------------|
-//! | General (GE) | [svd] | [svddc] | [least_squares] |
+//! Singular Value Decomposition
+//! -----------------------------
//!
-//! [svd]: svd/trait.SVD_.html#tymethod.svd
-//! [svddc]: svddck/trait.SVDDC_.html#tymethod.svddc
-//! [least_squares]: least_squares/trait.LeastSquaresSvdDivideConquer_.html#tymethod.least_squares
+//! - [svd] module for singular value decomposition (SVD) for general matrix
+//! - [svddc] module for singular value decomposition (SVD) with divided-and-conquer algorithm for general matrix
+//! - [least_squares] module for solving least square problem using SVD
+//!
+
+#![deny(rustdoc::broken_intra_doc_links, rustdoc::private_intra_doc_links)]
#[cfg(any(feature = "intel-mkl-system", feature = "intel-mkl-static"))]
extern crate intel_mkl_src as _src;
@@ -68,115 +85,482 @@ extern crate openblas_src as _src;
#[cfg(any(feature = "netlib-system", feature = "netlib-static"))]
extern crate netlib_src as _src;
+pub mod alloc;
+pub mod cholesky;
+pub mod eig;
+pub mod eig_generalized;
+pub mod eigh;
+pub mod eigh_generalized;
pub mod error;
+pub mod flags;
pub mod layout;
+pub mod least_squares;
+pub mod opnorm;
+pub mod qr;
+pub mod rcond;
+pub mod solve;
+pub mod solveh;
+pub mod svd;
+pub mod svddc;
+pub mod triangular;
+pub mod tridiagonal;
-mod cholesky;
-mod eig;
-mod eigh;
-mod least_squares;
-mod opnorm;
-mod qr;
-mod rcond;
-mod solve;
-mod solveh;
-mod svd;
-mod svddc;
-mod triangular;
-mod tridiagonal;
-
-pub use self::cholesky::*;
-pub use self::eig::*;
-pub use self::eigh::*;
-pub use self::least_squares::*;
-pub use self::opnorm::*;
-pub use self::qr::*;
-pub use self::rcond::*;
-pub use self::solve::*;
-pub use self::solveh::*;
-pub use self::svd::*;
-pub use self::svddc::*;
-pub use self::triangular::*;
-pub use self::tridiagonal::*;
+pub use crate::eig_generalized::GeneralizedEigenvalue;
+pub use self::flags::*;
+pub use self::least_squares::LeastSquaresOwned;
+pub use self::svd::{SvdOwned, SvdRef};
+pub use self::tridiagonal::{LUFactorizedTridiagonal, Tridiagonal};
+
+use self::{alloc::*, error::*, layout::*};
use cauchy::*;
+use std::mem::MaybeUninit;
pub type Pivot = Vec;
+#[cfg_attr(doc, katexit::katexit)]
/// Trait for primitive types which implements LAPACK subroutines
-pub trait Lapack:
- OperatorNorm_
- + QR_
- + SVD_
- + SVDDC_
- + Solve_
- + Solveh_
- + Cholesky_
- + Eig_
- + Eigh_
- + Triangular_
- + Tridiagonal_
- + Rcond_
- + LeastSquaresSvdDivideConquer_
-{
-}
+pub trait Lapack: Scalar {
+ /// Compute right eigenvalue and eigenvectors for a general matrix
+ fn eig(
+ calc_v: bool,
+ l: MatrixLayout,
+ a: &mut [Self],
+ ) -> Result<(Vec, Vec)>;
-impl Lapack for f32 {}
-impl Lapack for f64 {}
-impl Lapack for c32 {}
-impl Lapack for c64 {}
-
-/// Upper/Lower specification for seveal usages
-#[derive(Debug, Clone, Copy)]
-#[repr(u8)]
-pub enum UPLO {
- Upper = b'U',
- Lower = b'L',
-}
+ /// Compute right eigenvalue and eigenvectors for a general matrix
+ fn eig_generalized(
+ calc_v: bool,
+ l: MatrixLayout,
+ a: &mut [Self],
+ b: &mut [Self],
+ thresh_opt: Option,
+ ) -> Result<(
+ Vec>,
+ Vec,
+ )>;
-impl UPLO {
- pub fn t(self) -> Self {
- match self {
- UPLO::Upper => UPLO::Lower,
- UPLO::Lower => UPLO::Upper,
- }
- }
-}
+ /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
+ fn eigh(
+ calc_eigenvec: bool,
+ layout: MatrixLayout,
+ uplo: UPLO,
+ a: &mut [Self],
+ ) -> Result>;
-#[derive(Debug, Clone, Copy)]
-#[repr(u8)]
-pub enum Transpose {
- No = b'N',
- Transpose = b'T',
- Hermite = b'C',
-}
+ /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
+ fn eigh_generalized(
+ calc_eigenvec: bool,
+ layout: MatrixLayout,
+ uplo: UPLO,
+ a: &mut [Self],
+ b: &mut [Self],
+ ) -> Result>;
-#[derive(Debug, Clone, Copy)]
-#[repr(u8)]
-pub enum NormType {
- One = b'O',
- Infinity = b'I',
- Frobenius = b'F',
-}
+ /// Execute Householder reflection as the first step of QR-decomposition
+ ///
+ /// For C-continuous array,
+ /// this will call LQ-decomposition of the transposed matrix $ A^T = LQ^T $
+ fn householder(l: MatrixLayout, a: &mut [Self]) -> Result>;
-impl NormType {
- pub fn transpose(self) -> Self {
- match self {
- NormType::One => NormType::Infinity,
- NormType::Infinity => NormType::One,
- NormType::Frobenius => NormType::Frobenius,
- }
- }
+ /// Reconstruct Q-matrix from Householder-reflectors
+ fn q(l: MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>;
+
+ /// Execute QR-decomposition at once
+ fn qr(l: MatrixLayout, a: &mut [Self]) -> Result>;
+
+ /// Compute singular-value decomposition (SVD)
+ fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result>;
+
+ /// Compute singular value decomposition (SVD) with divide-and-conquer algorithm
+ fn svddc(layout: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result>;
+
+ /// Compute a vector $x$ which minimizes Euclidian norm $\| Ax - b\|$
+ /// for a given matrix $A$ and a vector $b$.
+ fn least_squares(
+ a_layout: MatrixLayout,
+ a: &mut [Self],
+ b: &mut [Self],
+ ) -> Result>;
+
+ /// Solve least square problems $\argmin_X \| AX - B\|$
+ fn least_squares_nrhs(
+ a_layout: MatrixLayout,
+ a: &mut [Self],
+ b_layout: MatrixLayout,
+ b: &mut [Self],
+ ) -> Result>;
+
+ /// Computes the LU decomposition of a general $m \times n$ matrix
+ /// with partial pivoting with row interchanges.
+ ///
+ /// For a given matrix $A$, LU decomposition is described as $A = PLU$ where:
+ ///
+ /// - $L$ is lower matrix
+ /// - $U$ is upper matrix
+ /// - $P$ is permutation matrix represented by [Pivot]
+ ///
+ /// This is designed as two step computation according to LAPACK API:
+ ///
+ /// 1. Factorize input matrix $A$ into $L$, $U$, and $P$.
+ /// 2. Solve linear equation $Ax = b$ by [Lapack::solve]
+ /// or compute inverse matrix $A^{-1}$ by [Lapack::inv] using the output of LU decomposition.
+ ///
+ /// Output
+ /// -------
+ /// - $U$ and $L$ are stored in `a` after LU decomposition has succeeded.
+ /// - $P$ is returned as [Pivot]
+ ///
+ /// Error
+ /// ------
+ /// - if the matrix is singular
+ /// - On this case, `return_code` in [Error::LapackComputationalFailure] means
+ /// `return_code`-th diagonal element of $U$ becomes zero.
+ ///
+ fn lu(l: MatrixLayout, a: &mut [Self]) -> Result;
+
+ /// Compute inverse matrix $A^{-1}$ from the output of LU-decomposition
+ fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()>;
+
+ /// Solve linear equations $Ax = b$ using the output of LU-decomposition
+ fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
+
+ /// Factorize symmetric/Hermitian matrix using Bunch-Kaufman diagonal pivoting method
+ ///
+ /// For a given symmetric matrix $A$,
+ /// this method factorizes $A = U^T D U$ or $A = L D L^T$ where
+ ///
+ /// - $U$ (or $L$) are is a product of permutation and unit upper (lower) triangular matrices
+ /// - $D$ is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
+ ///
+ /// This takes two-step approach based in LAPACK:
+ ///
+ /// 1. Factorize given matrix $A$ into upper ($U$) or lower ($L$) form with diagonal matrix $D$
+ /// 2. Then solve linear equation $Ax = b$, and/or calculate inverse matrix $A^{-1}$
+ ///
+ fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result;
+
+ /// Compute inverse matrix $A^{-1}$ using the result of [Lapack::bk]
+ fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()>;
+
+ /// Solve symmetric/Hermitian linear equation $Ax = b$ using the result of [Lapack::bk]
+ fn solveh(l: MatrixLayout, uplo: UPLO, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()>;
+
+ /// Solve symmetric/Hermitian positive-definite linear equations using Cholesky decomposition
+ ///
+ /// For a given positive definite matrix $A$,
+ /// Cholesky decomposition is described as $A = U^T U$ or $A = LL^T$ where
+ ///
+ /// - $L$ is lower matrix
+ /// - $U$ is upper matrix
+ ///
+ /// This is designed as two step computation according to LAPACK API
+ ///
+ /// 1. Factorize input matrix $A$ into $L$ or $U$
+ /// 2. Solve linear equation $Ax = b$ by [Lapack::solve_cholesky]
+ /// or compute inverse matrix $A^{-1}$ by [Lapack::inv_cholesky]
+ ///
+ fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+
+ /// Compute inverse matrix $A^{-1}$ using $U$ or $L$ calculated by [Lapack::cholesky]
+ fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+
+ /// Solve linear equation $Ax = b$ using $U$ or $L$ calculated by [Lapack::cholesky]
+ fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
+
+ /// Estimates the the reciprocal of the condition number of the matrix in 1-norm.
+ ///
+ /// `anorm` should be the 1-norm of the matrix `a`.
+ fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result;
+
+ /// Compute norm of matrices
+ ///
+ /// For a $n \times m$ matrix
+ /// $$
+ /// A = \begin{pmatrix}
+ /// a_{11} & \cdots & a_{1m} \\\\
+ /// \vdots & \ddots & \vdots \\\\
+ /// a_{n1} & \cdots & a_{nm}
+ /// \end{pmatrix}
+ /// $$
+ /// LAPACK can compute three types of norms:
+ ///
+ /// - Operator norm based on 1-norm for its domain linear space:
+ /// $$
+ /// \Vert A \Vert_1 = \sup_{\Vert x \Vert_1 = 1} \Vert Ax \Vert_1
+ /// = \max_{1 \le j \le m } \sum_{i=1}^n |a_{ij}|
+ /// $$
+ /// where
+ /// $\Vert x\Vert_1 = \sum_{j=1}^m |x_j|$
+ /// is 1-norm for a vector $x$.
+ ///
+ /// - Operator norm based on $\infty$-norm for its domain linear space:
+ /// $$
+ /// \Vert A \Vert_\infty = \sup_{\Vert x \Vert_\infty = 1} \Vert Ax \Vert_\infty
+ /// = \max_{1 \le i \le n } \sum_{j=1}^m |a_{ij}|
+ /// $$
+ /// where
+ /// $\Vert x\Vert_\infty = \max_{j=1}^m |x_j|$
+ /// is $\infty$-norm for a vector $x$.
+ ///
+ /// - Frobenious norm
+ /// $$
+ /// \Vert A \Vert_F = \sqrt{\mathrm{Tr} \left(AA^\dagger\right)} = \sqrt{\sum_{i=1}^n \sum_{j=1}^m |a_{ij}|^2}
+ /// $$
+ ///
+ fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
+
+ fn solve_triangular(
+ al: MatrixLayout,
+ bl: MatrixLayout,
+ uplo: UPLO,
+ d: Diag,
+ a: &[Self],
+ b: &mut [Self],
+ ) -> Result<()>;
+
+ /// Computes the LU factorization of a tridiagonal `m x n` matrix `a` using
+ /// partial pivoting with row interchanges.
+ fn lu_tridiagonal(a: Tridiagonal) -> Result>;
+
+ fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result;
+
+ fn solve_tridiagonal(
+ lu: &LUFactorizedTridiagonal,
+ bl: MatrixLayout,
+ t: Transpose,
+ b: &mut [Self],
+ ) -> Result<()>;
}
-/// Create a vector without initialization
-///
-/// Safety
-/// ------
-/// - Memory is not initialized. Do not read the memory before write.
-///
-unsafe fn vec_uninit(n: usize) -> Vec {
- let mut v = Vec::with_capacity(n);
- v.set_len(n);
- v
+macro_rules! impl_lapack {
+ ($s:ty) => {
+ impl Lapack for $s {
+ fn eig(
+ calc_v: bool,
+ l: MatrixLayout,
+ a: &mut [Self],
+ ) -> Result<(Vec, Vec)> {
+ use eig::*;
+ let work = EigWork::<$s>::new(calc_v, l)?;
+ let EigOwned { eigs, vr, vl } = work.eval(a)?;
+ Ok((eigs, vr.or(vl).unwrap_or_default()))
+ }
+
+ fn eig_generalized(
+ calc_v: bool,
+ l: MatrixLayout,
+ a: &mut [Self],
+ b: &mut [Self],
+ thresh_opt: Option,
+ ) -> Result<(
+ Vec>,
+ Vec,
+ )> {
+ use eig_generalized::*;
+ let work = EigGeneralizedWork::<$s>::new(calc_v, l)?;
+ let eig_generalized_owned = work.eval(a, b)?;
+ let eigs = eig_generalized_owned.calc_eigs(thresh_opt);
+ let vr = eig_generalized_owned.vr;
+ let vl = eig_generalized_owned.vl;
+ Ok((eigs, vr.or(vl).unwrap_or_default()))
+ }
+
+ fn eigh(
+ calc_eigenvec: bool,
+ layout: MatrixLayout,
+ uplo: UPLO,
+ a: &mut [Self],
+ ) -> Result> {
+ use eigh::*;
+ let work = EighWork::<$s>::new(calc_eigenvec, layout)?;
+ work.eval(uplo, a)
+ }
+
+ fn eigh_generalized(
+ calc_eigenvec: bool,
+ layout: MatrixLayout,
+ uplo: UPLO,
+ a: &mut [Self],
+ b: &mut [Self],
+ ) -> Result> {
+ use eigh_generalized::*;
+ let work = EighGeneralizedWork::<$s>::new(calc_eigenvec, layout)?;
+ work.eval(uplo, a, b)
+ }
+
+ fn householder(l: MatrixLayout, a: &mut [Self]) -> Result> {
+ use qr::*;
+ let work = HouseholderWork::<$s>::new(l)?;
+ work.eval(a)
+ }
+
+ fn q(l: MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()> {
+ use qr::*;
+ let mut work = QWork::<$s>::new(l)?;
+ work.calc(a, tau)?;
+ Ok(())
+ }
+
+ fn qr(l: MatrixLayout, a: &mut [Self]) -> Result> {
+ let tau = Self::householder(l, a)?;
+ let r = Vec::from(&*a);
+ Self::q(l, a, &tau)?;
+ Ok(r)
+ }
+
+ fn svd(
+ l: MatrixLayout,
+ calc_u: bool,
+ calc_vt: bool,
+ a: &mut [Self],
+ ) -> Result> {
+ use svd::*;
+ let work = SvdWork::<$s>::new(l, calc_u, calc_vt)?;
+ work.eval(a)
+ }
+
+ fn svddc(layout: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result> {
+ use svddc::*;
+ let work = SvdDcWork::<$s>::new(layout, jobz)?;
+ work.eval(a)
+ }
+
+ fn least_squares(
+ l: MatrixLayout,
+ a: &mut [Self],
+ b: &mut [Self],
+ ) -> Result> {
+ let b_layout = l.resized(b.len() as i32, 1);
+ Self::least_squares_nrhs(l, a, b_layout, b)
+ }
+
+ fn least_squares_nrhs(
+ a_layout: MatrixLayout,
+ a: &mut [Self],
+ b_layout: MatrixLayout,
+ b: &mut [Self],
+ ) -> Result> {
+ use least_squares::*;
+ let work = LeastSquaresWork::<$s>::new(a_layout, b_layout)?;
+ work.eval(a, b)
+ }
+
+ fn lu(l: MatrixLayout, a: &mut [Self]) -> Result {
+ use solve::*;
+ LuImpl::lu(l, a)
+ }
+
+ fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()> {
+ use solve::*;
+ let mut work = InvWork::<$s>::new(l)?;
+ work.calc(a, p)?;
+ Ok(())
+ }
+
+ fn solve(
+ l: MatrixLayout,
+ t: Transpose,
+ a: &[Self],
+ p: &Pivot,
+ b: &mut [Self],
+ ) -> Result<()> {
+ use solve::*;
+ SolveImpl::solve(l, t, a, p, b)
+ }
+
+ fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result {
+ use solveh::*;
+ let work = BkWork::<$s>::new(l)?;
+ work.eval(uplo, a)
+ }
+
+ fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
+ use solveh::*;
+ let mut work = InvhWork::<$s>::new(l)?;
+ work.calc(uplo, a, ipiv)
+ }
+
+ fn solveh(
+ l: MatrixLayout,
+ uplo: UPLO,
+ a: &[Self],
+ ipiv: &Pivot,
+ b: &mut [Self],
+ ) -> Result<()> {
+ use solveh::*;
+ SolvehImpl::solveh(l, uplo, a, ipiv, b)
+ }
+
+ fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
+ use cholesky::*;
+ CholeskyImpl::cholesky(l, uplo, a)
+ }
+
+ fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
+ use cholesky::*;
+ InvCholeskyImpl::inv_cholesky(l, uplo, a)
+ }
+
+ fn solve_cholesky(
+ l: MatrixLayout,
+ uplo: UPLO,
+ a: &[Self],
+ b: &mut [Self],
+ ) -> Result<()> {
+ use cholesky::*;
+ SolveCholeskyImpl::solve_cholesky(l, uplo, a, b)
+ }
+
+ fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result {
+ use rcond::*;
+ let mut work = RcondWork::<$s>::new(l);
+ work.calc(a, anorm)
+ }
+
+ fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
+ use opnorm::*;
+ let mut work = OperatorNormWork::<$s>::new(t, l);
+ work.calc(a)
+ }
+
+ fn solve_triangular(
+ al: MatrixLayout,
+ bl: MatrixLayout,
+ uplo: UPLO,
+ d: Diag,
+ a: &[Self],
+ b: &mut [Self],
+ ) -> Result<()> {
+ use triangular::*;
+ SolveTriangularImpl::solve_triangular(al, bl, uplo, d, a, b)
+ }
+
+ fn lu_tridiagonal(a: Tridiagonal) -> Result> {
+ use tridiagonal::*;
+ let work = LuTridiagonalWork::<$s>::new(a.l);
+ work.eval(a)
+ }
+
+ fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal) -> Result {
+ use tridiagonal::*;
+ let mut work = RcondTridiagonalWork::<$s>::new(lu.a.l);
+ work.calc(lu)
+ }
+
+ fn solve_tridiagonal(
+ lu: &LUFactorizedTridiagonal,
+ bl: MatrixLayout,
+ t: Transpose,
+ b: &mut [Self],
+ ) -> Result<()> {
+ use tridiagonal::*;
+ SolveTridiagonalImpl::solve_tridiagonal(lu, bl, t, b)
+ }
+ }
+ };
}
+impl_lapack!(c64);
+impl_lapack!(c32);
+impl_lapack!(f64);
+impl_lapack!(f32);
diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs
index dd84f441..1789f385 100644
--- a/lax/src/opnorm.rs
+++ b/lax/src/opnorm.rs
@@ -1,35 +1,58 @@
-//! Operator norms of matrices
+//! Operator norm
-use super::NormType;
+use super::{AsPtr, NormType};
use crate::{layout::MatrixLayout, *};
use cauchy::*;
-pub trait OperatorNorm_: Scalar {
- fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
+pub struct OperatorNormWork {
+ pub ty: NormType,
+ pub layout: MatrixLayout,
+ pub work: Vec>,
}
-macro_rules! impl_opnorm {
- ($scalar:ty, $lange:path) => {
- impl OperatorNorm_ for $scalar {
- fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
- let m = l.lda();
- let n = l.len();
- let t = match l {
- MatrixLayout::F { .. } => t,
- MatrixLayout::C { .. } => t.transpose(),
+pub trait OperatorNormWorkImpl {
+ type Elem: Scalar;
+ fn new(t: NormType, l: MatrixLayout) -> Self;
+ fn calc(&mut self, a: &[Self::Elem]) -> ::Real;
+}
+
+macro_rules! impl_operator_norm {
+ ($s:ty, $lange:path) => {
+ impl OperatorNormWorkImpl for OperatorNormWork<$s> {
+ type Elem = $s;
+
+ fn new(ty: NormType, layout: MatrixLayout) -> Self {
+ let m = layout.lda();
+ let work = match (ty, layout) {
+ (NormType::Infinity, MatrixLayout::F { .. })
+ | (NormType::One, MatrixLayout::C { .. }) => vec_uninit(m as usize),
+ _ => Vec::new(),
};
- let mut work = if matches!(t, NormType::Infinity) {
- unsafe { vec_uninit(m as usize) }
- } else {
- Vec::new()
+ OperatorNormWork { ty, layout, work }
+ }
+
+ fn calc(&mut self, a: &[Self::Elem]) -> ::Real {
+ let m = self.layout.lda();
+ let n = self.layout.len();
+ let t = match self.layout {
+ MatrixLayout::F { .. } => self.ty,
+ MatrixLayout::C { .. } => self.ty.transpose(),
};
- unsafe { $lange(t as u8, m, n, a, m, &mut work) }
+ unsafe {
+ $lange(
+ t.as_ptr(),
+ &m,
+ &n,
+ AsPtr::as_ptr(a),
+ &m,
+ AsPtr::as_mut_ptr(&mut self.work),
+ )
+ }
}
}
};
-} // impl_opnorm!
-
-impl_opnorm!(f64, lapack::dlange);
-impl_opnorm!(f32, lapack::slange);
-impl_opnorm!(c64, lapack::zlange);
-impl_opnorm!(c32, lapack::clange);
+}
+impl_operator_norm!(c64, lapack_sys::zlange_);
+impl_operator_norm!(c32, lapack_sys::clange_);
+impl_operator_norm!(f64, lapack_sys::dlange_);
+impl_operator_norm!(f32, lapack_sys::slange_);
diff --git a/lax/src/qr.rs b/lax/src/qr.rs
index 6460b8b9..f37bd579 100644
--- a/lax/src/qr.rs
+++ b/lax/src/qr.rs
@@ -4,152 +4,213 @@ use crate::{error::*, layout::MatrixLayout, *};
use cauchy::*;
use num_traits::{ToPrimitive, Zero};
-pub trait QR_: Sized {
- /// Execute Householder reflection as the first step of QR-decomposition
- ///
- /// For C-continuous array,
- /// this will call LQ-decomposition of the transposed matrix $ A^T = LQ^T $
- fn householder(l: MatrixLayout, a: &mut [Self]) -> Result>;
-
- /// Reconstruct Q-matrix from Householder-reflectors
- fn q(l: MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>;
+pub struct HouseholderWork {
+ pub m: i32,
+ pub n: i32,
+ pub layout: MatrixLayout,
+ pub tau: Vec>,
+ pub work: Vec>,
+}
- /// Execute QR-decomposition at once
- fn qr(l: MatrixLayout, a: &mut [Self]) -> Result>;
+pub trait HouseholderWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(l: MatrixLayout) -> Result;
+ fn calc(&mut self, a: &mut [Self::Elem]) -> Result<&[Self::Elem]>;
+ fn eval(self, a: &mut [Self::Elem]) -> Result>;
}
-macro_rules! impl_qr {
- ($scalar:ty, $qrf:path, $lqf:path, $gqr:path, $glq:path) => {
- impl QR_ for $scalar {
- fn householder(l: MatrixLayout, mut a: &mut [Self]) -> Result> {
- let m = l.lda();
- let n = l.len();
- let k = m.min(n);
- let mut tau = unsafe { vec_uninit(k as usize) };
+macro_rules! impl_householder_work {
+ ($s:ty, $qrf:path, $lqf: path) => {
+ impl HouseholderWorkImpl for HouseholderWork<$s> {
+ type Elem = $s;
- // eval work size
+ fn new(layout: MatrixLayout) -> Result {
+ let m = layout.lda();
+ let n = layout.len();
+ let k = m.min(n);
+ let mut tau = vec_uninit(k as usize);
let mut info = 0;
- let mut work_size = [Self::zero()];
- unsafe {
- match l {
- MatrixLayout::F { .. } => {
- $qrf(m, n, &mut a, m, &mut tau, &mut work_size, -1, &mut info);
- }
- MatrixLayout::C { .. } => {
- $lqf(m, n, &mut a, m, &mut tau, &mut work_size, -1, &mut info);
- }
- }
+ let mut work_size = [Self::Elem::zero()];
+ match layout {
+ MatrixLayout::F { .. } => unsafe {
+ $qrf(
+ &m,
+ &n,
+ std::ptr::null_mut(),
+ &m,
+ AsPtr::as_mut_ptr(&mut tau),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ )
+ },
+ MatrixLayout::C { .. } => unsafe {
+ $lqf(
+ &m,
+ &n,
+ std::ptr::null_mut(),
+ &m,
+ AsPtr::as_mut_ptr(&mut tau),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ )
+ },
}
info.as_lapack_result()?;
-
- // calc
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
- unsafe {
- match l {
- MatrixLayout::F { .. } => {
- $qrf(
- m,
- n,
- &mut a,
- m,
- &mut tau,
- &mut work,
- lwork as i32,
- &mut info,
- );
- }
- MatrixLayout::C { .. } => {
- $lqf(
- m,
- n,
- &mut a,
- m,
- &mut tau,
- &mut work,
- lwork as i32,
- &mut info,
- );
- }
- }
+ let work = vec_uninit(lwork);
+ Ok(HouseholderWork {
+ n,
+ m,
+ layout,
+ tau,
+ work,
+ })
+ }
+
+ fn calc(&mut self, a: &mut [Self::Elem]) -> Result<&[Self::Elem]> {
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ match self.layout {
+ MatrixLayout::F { .. } => unsafe {
+ $qrf(
+ &self.m,
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.m,
+ AsPtr::as_mut_ptr(&mut self.tau),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ );
+ },
+ MatrixLayout::C { .. } => unsafe {
+ $lqf(
+ &self.m,
+ &self.n,
+ AsPtr::as_mut_ptr(a),
+ &self.m,
+ AsPtr::as_mut_ptr(&mut self.tau),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ );
+ },
}
info.as_lapack_result()?;
+ Ok(unsafe { self.tau.slice_assume_init_ref() })
+ }
- Ok(tau)
+ fn eval(mut self, a: &mut [Self::Elem]) -> Result> {
+ let _eig = self.calc(a)?;
+ Ok(unsafe { self.tau.assume_init() })
}
+ }
+ };
+}
+impl_householder_work!(c64, lapack_sys::zgeqrf_, lapack_sys::zgelqf_);
+impl_householder_work!(c32, lapack_sys::cgeqrf_, lapack_sys::cgelqf_);
+impl_householder_work!(f64, lapack_sys::dgeqrf_, lapack_sys::dgelqf_);
+impl_householder_work!(f32, lapack_sys::sgeqrf_, lapack_sys::sgelqf_);
- fn q(l: MatrixLayout, mut a: &mut [Self], tau: &[Self]) -> Result<()> {
- let m = l.lda();
- let n = l.len();
- let k = m.min(n);
- assert_eq!(tau.len(), k as usize);
+pub struct QWork {
+ pub layout: MatrixLayout,
+ pub work: Vec>,
+}
- // eval work size
- let mut info = 0;
- let mut work_size = [Self::zero()];
- unsafe {
- match l {
- MatrixLayout::F { .. } => {
- $gqr(m, k, k, &mut a, m, &tau, &mut work_size, -1, &mut info)
- }
- MatrixLayout::C { .. } => {
- $glq(k, n, k, &mut a, m, &tau, &mut work_size, -1, &mut info)
- }
- }
- };
+pub trait QWorkImpl: Sized {
+ type Elem: Scalar;
+ fn new(layout: MatrixLayout) -> Result;
+ fn calc(&mut self, a: &mut [Self::Elem], tau: &[Self::Elem]) -> Result<()>;
+}
- // calc
+macro_rules! impl_q_work {
+ ($s:ty, $gqr:path, $glq:path) => {
+ impl QWorkImpl for QWork<$s> {
+ type Elem = $s;
+
+ fn new(layout: MatrixLayout) -> Result {
+ let m = layout.lda();
+ let n = layout.len();
+ let k = m.min(n);
+ let mut info = 0;
+ let mut work_size = [Self::Elem::zero()];
+ match layout {
+ MatrixLayout::F { .. } => unsafe {
+ $gqr(
+ &m,
+ &k,
+ &k,
+ std::ptr::null_mut(),
+ &m,
+ std::ptr::null_mut(),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ )
+ },
+ MatrixLayout::C { .. } => unsafe {
+ $glq(
+ &k,
+ &n,
+ &k,
+ std::ptr::null_mut(),
+ &m,
+ std::ptr::null_mut(),
+ AsPtr::as_mut_ptr(&mut work_size),
+ &(-1),
+ &mut info,
+ )
+ },
+ }
let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
- unsafe {
- match l {
- MatrixLayout::F { .. } => {
- $gqr(m, k, k, &mut a, m, &tau, &mut work, lwork as i32, &mut info)
- }
- MatrixLayout::C { .. } => {
- $glq(k, n, k, &mut a, m, &tau, &mut work, lwork as i32, &mut info)
- }
- }
+ let work = vec_uninit(lwork);
+ Ok(QWork { layout, work })
+ }
+
+ fn calc(&mut self, a: &mut [Self::Elem], tau: &[Self::Elem]) -> Result<()> {
+ let m = self.layout.lda();
+ let n = self.layout.len();
+ let k = m.min(n);
+ let lwork = self.work.len().to_i32().unwrap();
+ let mut info = 0;
+ match self.layout {
+ MatrixLayout::F { .. } => unsafe {
+ $gqr(
+ &m,
+ &k,
+ &k,
+ AsPtr::as_mut_ptr(a),
+ &m,
+ AsPtr::as_ptr(&tau),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ )
+ },
+ MatrixLayout::C { .. } => unsafe {
+ $glq(
+ &k,
+ &n,
+ &k,
+ AsPtr::as_mut_ptr(a),
+ &m,
+ AsPtr::as_ptr(&tau),
+ AsPtr::as_mut_ptr(&mut self.work),
+ &lwork,
+ &mut info,
+ )
+ },
}
info.as_lapack_result()?;
Ok(())
}
-
- fn qr(l: MatrixLayout, a: &mut [Self]) -> Result> {
- let tau = Self::householder(l, a)?;
- let r = Vec::from(&*a);
- Self::q(l, a, &tau)?;
- Ok(r)
- }
}
};
-} // endmacro
+}
-impl_qr!(
- f64,
- lapack::dgeqrf,
- lapack::dgelqf,
- lapack::dorgqr,
- lapack::dorglq
-);
-impl_qr!(
- f32,
- lapack::sgeqrf,
- lapack::sgelqf,
- lapack::sorgqr,
- lapack::sorglq
-);
-impl_qr!(
- c64,
- lapack::zgeqrf,
- lapack::zgelqf,
- lapack::zungqr,
- lapack::zunglq
-);
-impl_qr!(
- c32,
- lapack::cgeqrf,
- lapack::cgelqf,
- lapack::cungqr,
- lapack::cunglq
-);
+impl_q_work!(c64, lapack_sys::zungqr_, lapack_sys::zunglq_);
+impl_q_work!(c32, lapack_sys::cungqr_, lapack_sys::cunglq_);
+impl_q_work!(f64, lapack_sys::dorgqr_, lapack_sys::dorglq_);
+impl_q_work!(f32, lapack_sys::sorgqr_, lapack_sys::sorglq_);
diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs
index 91d7458c..4d4a4c92 100644
--- a/lax/src/rcond.rs
+++ b/lax/src/rcond.rs
@@ -1,85 +1,124 @@
+//! Reciprocal conditional number
+
use crate::{error::*, layout::MatrixLayout, *};
use cauchy::*;
use num_traits::Zero;
-pub trait Rcond_: Scalar + Sized {
- /// Estimates the the reciprocal of the condition number of the matrix in 1-norm.
- ///
- /// `anorm` should be the 1-norm of the matrix `a`.
- fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result;
+pub struct RcondWork {
+ pub layout: MatrixLayout,
+ pub work: Vec>,
+ pub rwork: Option>>,
+ pub iwork: Option>>,
}
-macro_rules! impl_rcond_real {
- ($scalar:ty, $gecon:path) => {
- impl Rcond_ for $scalar {
- fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result {
- let (n, _) = l.size();
- let mut rcond = Self::Real::zero();
- let mut info = 0;
+pub trait RcondWorkImpl {
+ type Elem: Scalar;
+ fn new(l: MatrixLayout) -> Self;
+ fn calc(
+ &mut self,
+ a: &[Self::Elem],
+ anorm: ::Real,
+ ) -> Result<::Real>;
+}
+
+macro_rules! impl_rcond_work_c {
+ ($c:ty, $con:path) => {
+ impl RcondWorkImpl for RcondWork<$c> {
+ type Elem = $c;
- let mut work = unsafe { vec_uninit(4 * n as usize) };
- let mut iwork = unsafe { vec_uninit(n as usize) };
- let norm_type = match l {
+ fn new(layout: MatrixLayout) -> Self {
+ let (n, _) = layout.size();
+ let work = vec_uninit(2 * n as usize);
+ let rwork = vec_uninit(2 * n as usize);
+ RcondWork {
+ layout,
+ work,
+ rwork: Some(rwork),
+ iwork: None,
+ }
+ }
+
+ fn calc(
+ &mut self,
+ a: &[Self::Elem],
+ anorm: ::Real,
+ ) -> Result<::Real> {
+ let (n, _) = self.layout.size();
+ let mut rcond = ::Real::zero();
+ let mut info = 0;
+ let norm_type = match self.layout {
MatrixLayout::C { .. } => NormType::Infinity,
MatrixLayout::F { .. } => NormType::One,
- } as u8;
+ };
unsafe {
- $gecon(
- norm_type,
- n,
- a,
- l.lda(),
- anorm,
+ $con(
+ norm_type.as_ptr(),
+ &n,
+ AsPtr::as_ptr(a),
+ &self.layout.lda(),
+ &anorm,
&mut rcond,
- &mut work,
- &mut iwork,
+ AsPtr::as_mut_ptr(&mut self.work),
+ AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
&mut info,
)
};
info.as_lapack_result()?;
-
Ok(rcond)
}
}
};
}
+impl_rcond_work_c!(c64, lapack_sys::zgecon_);
+impl_rcond_work_c!(c32, lapack_sys::cgecon_);
-impl_rcond_real!(f32, lapack::sgecon);
-impl_rcond_real!(f64, lapack::dgecon);
+macro_rules! impl_rcond_work_r {
+ ($r:ty, $con:path) => {
+ impl RcondWorkImpl for RcondWork<$r> {
+ type Elem = $r;
+
+ fn new(layout: MatrixLayout) -> Self {
+ let (n, _) = layout.size();
+ let work = vec_uninit(4 * n as usize);
+ let iwork = vec_uninit(n as usize);
+ RcondWork {
+ layout,
+ work,
+ rwork: None,
+ iwork: Some(iwork),
+ }
+ }
-macro_rules! impl_rcond_complex {
- ($scalar:ty, $gecon:path) => {
- impl Rcond_ for $scalar {
- fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result {
- let (n, _) = l.size();
- let mut rcond = Self::Real::zero();
+ fn calc(
+ &mut self,
+ a: &[Self::Elem],
+ anorm: ::Real,
+ ) -> Result<::Real> {
+ let (n, _) = self.layout.size();
+ let mut rcond = ::Real::zero();
let mut info = 0;
- let mut work = unsafe { vec_uninit(2 * n as usize) };
- let mut rwork = unsafe { vec_uninit(2 * n as usize) };
- let norm_type = match l {
+ let norm_type = match self.layout {
MatrixLayout::C { .. } => NormType::Infinity,
MatrixLayout::F { .. } => NormType::One,
- } as u8;
+ };
unsafe {
- $gecon(
- norm_type,
- n,
- a,
- l.lda(),
- anorm,
+ $con(
+ norm_type.as_ptr(),
+ &n,
+ AsPtr::as_ptr(a),
+ &self.layout.lda(),
+ &anorm,
&mut rcond,
- &mut work,
- &mut rwork,
+ AsPtr::as_mut_ptr(&mut self.work),
+ AsPtr::as_mut_ptr(self.iwork.as_mut().unwrap()),
&mut info,
)
};
info.as_lapack_result()?;
-
Ok(rcond)
}
}
};
}
-
-impl_rcond_complex!(c32, lapack::cgecon);
-impl_rcond_complex!(c64, lapack::zgecon);
+impl_rcond_work_r!(f64, lapack_sys::dgecon_);
+impl_rcond_work_r!(f32, lapack_sys::sgecon_);
diff --git a/lax/src/solve.rs b/lax/src/solve.rs
index 39498a04..63f69983 100644
--- a/lax/src/solve.rs
+++ b/lax/src/solve.rs
@@ -1,30 +1,25 @@
-//! Solve linear problem using LU decomposition
+//! Solve linear equations using LU-decomposition
use crate::{error::*, layout::MatrixLayout, *};
use cauchy::*;
use num_traits::{ToPrimitive, Zero};
-pub trait Solve_: Scalar + Sized {
- /// Computes the LU factorization of a general `m x n` matrix `a` using
- /// partial pivoting with row interchanges.
- ///
- /// $ PA = LU $
- ///
- /// Error
- /// ------
- /// - `LapackComputationalFailure { return_code }` when the matrix is singular
- /// - Division by zero will occur if it is used to solve a system of equations
- /// because `U[(return_code-1, return_code-1)]` is exactly zero.
+/// Helper trait to abstract `*getrf` LAPACK routines for implementing [Lapack::lu]
+///
+/// LAPACK correspondance
+/// ----------------------
+///
+/// | f32 | f64 | c32 | c64 |
+/// |:-------|:-------|:-------|:-------|
+/// | sgetrf | dgetrf | cgetrf | zgetrf |
+///
+pub trait LuImpl: Scalar {
fn lu(l: MatrixLayout, a: &mut [Self]) -> Result;
-
- fn inv(l: MatrixLayout, a: &mut [Self], p: &Pivot) -> Result<()>;
-
- fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
}
-macro_rules! impl_solve {
- ($scalar:ty, $getrf:path, $getri:path, $getrs:path) => {
- impl Solve_ for $scalar {
+macro_rules! impl_lu {
+ ($scalar:ty, $getrf:path) => {
+ impl LuImpl for $scalar {
fn lu(l: MatrixLayout, a: &mut [Self]) -> Result {
let (row, col) = l.size();
assert_eq!(a.len() as i32, row * col);
@@ -33,41 +28,71 @@ macro_rules! impl_solve {
return Ok(Vec::new());
}
let k = ::std::cmp::min(row, col);
- let mut ipiv = unsafe { vec_uninit(k as usize) };
- let mut info = 0;
- unsafe { $getrf(l.lda(), l.len(), a, l.lda(), &mut ipiv, &mut info) };
- info.as_lapack_result()?;
- Ok(ipiv)
- }
-
- fn inv(l: MatrixLayout, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
- let (n, _) = l.size();
-
- // calc work size
+ let mut ipiv = vec_uninit(k as usize);
let mut info = 0;
- let mut work_size = [Self::zero()];
- unsafe { $getri(n, a, l.lda(), ipiv, &mut work_size, -1, &mut info) };
- info.as_lapack_result()?;
-
- // actual
- let lwork = work_size[0].to_usize().unwrap();
- let mut work = unsafe { vec_uninit(lwork) };
unsafe {
- $getri(
- l.len(),
- a,
- l.lda(),
- ipiv,
- &mut work,
- lwork as i32,
+ $getrf(
+ &l.lda(),
+ &l.len(),
+ AsPtr::as_mut_ptr(a),
+ &l.lda(),
+ AsPtr::as_mut_ptr(&mut ipiv),
&mut info,
)
};
info.as_lapack_result()?;
-
- Ok(())
+ let ipiv = unsafe { ipiv.assume_init() };
+ Ok(ipiv)
}
+ }
+ };
+}
+
+impl_lu!(c64, lapack_sys::zgetrf_);
+impl_lu!(c32, lapack_sys::cgetrf_);
+impl_lu!(f64, lapack_sys::dgetrf_);
+impl_lu!(f32, lapack_sys::sgetrf_);
+
+#[cfg_attr(doc, katexit::katexit)]
+/// Helper trait to abstract `*getrs` LAPACK routines for implementing [Lapack::solve]
+///
+/// If the array has C layout, then it needs to be handled
+/// specially, since LAPACK expects a Fortran-layout array.
+/// Reinterpreting a C layout array as Fortran layout is
+/// equivalent to transposing it. So, we can handle the "no
+/// transpose" and "transpose" cases by swapping to "transpose"
+/// or "no transpose", respectively. For the "Hermite" case, we
+/// can take advantage of the following:
+///
+/// $$
+/// \begin{align*}
+/// A^H x &= b \\\\
+/// \Leftrightarrow \overline{A^T} x &= b \\\\
+/// \Leftrightarrow \overline{\overline{A^T} x} &= \overline{b} \\\\
+/// \Leftrightarrow \overline{\overline{A^T}} \overline{x} &= \overline{b} \\\\
+/// \Leftrightarrow A^T \overline{x} &= \overline{b}
+/// \end{align*}
+/// $$
+///
+/// So, we can handle this case by switching to "no transpose"
+/// (which is equivalent to transposing the array since it will
+/// be reinterpreted as Fortran layout) and applying the
+/// elementwise conjugate to `x` and `b`.
+///
+pub trait SolveImpl: Scalar {
+ /// LAPACK correspondance
+ /// ----------------------
+ ///
+ /// | f32 | f64 | c32 | c64 |
+ /// |:-------|:-------|:-------|:-------|
+ /// | sgetrs | dgetrs | cgetrs | zgetrs |
+ ///
+ fn solve(l: MatrixLayout, t: Transpose, a: &[Self], p: &Pivot, b: &mut [Self]) -> Result<()>;
+}
+macro_rules! impl_solve {
+ ($scalar:ty, $getrs:path) => {
+ impl SolveImpl for $scalar {
fn solve(
l: MatrixLayout,
t: Transpose,
@@ -75,18 +100,41 @@ macro_rules! impl_solve {
ipiv: &Pivot,
b: &mut [Self],
) -> Result<()> {
- let t = match l {
+ let (t, conj) = match l {
MatrixLayout::C { .. } => match t {
- Transpose::No => Transpose::Transpose,
- Transpose::Transpose | Transpose::Hermite => Transpose::No,
+ Transpose::No => (Transpose::Transpose, false),
+ Transpose::Transpose => (Transpose::No, false),
+ Transpose::Hermite => (Transpose::No, true),
},
- _ => t,
+ MatrixLayout::F { .. } => (t, false),
};
let (n, _) = l.size();
let nrhs = 1;
let ldb = l.lda();
let mut info = 0;
- unsafe { $getrs(t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb, &mut info) };
+ if conj {
+ for b_elem in &mut *b {
+ *b_elem = b_elem.conj();
+ }
+ }
+ unsafe {
+ $getrs(
+ t.as_ptr(),
+ &n,
+ &nrhs,
+ AsPtr::as_ptr(a),
+ &l.lda(),
+ ipiv.as_ptr(),
+ AsPtr::as_mut_ptr(b),
+ &ldb,
+ &mut info,
+ )
+ };
+ if conj {
+ for b_elem in &mut *b {
+ *b_elem = b_elem.conj();
+ }
+ }
info.as_lapack_result()?;
Ok(())
}
@@ -94,7 +142,83 @@ macro_rules! impl_solve {
};
} // impl_solve!
-impl_solve!(f64, lapack::dgetrf, lapack::dgetri, lapack::dgetrs);
-impl_solve!(f32, lapack::sgetrf, lapack::sgetri, lapack::sgetrs);
-impl_solve!(c64, lapack::zgetrf, lapack::zgetri, lapack::zgetrs);
-impl_solve!(c32, lapack::cgetrf, lapack::cgetri, lapack::cgetrs);
+impl_solve!(f64, lapack_sys::dgetrs_);
+impl_solve!(f32, lapack_sys::sgetrs_);
+impl_solve!(c64, lapack_sys::zgetrs_);
+impl_solve!(c32, lapack_sys::cgetrs_);
+
+/// Working memory for computing inverse matrix
+pub struct InvWork {
+ pub layout: MatrixLayout,
+ pub work: Vec