Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generic scalar implementation #7

Open
wants to merge 12 commits into
base: master
Choose a base branch
from

Conversation

hmunozb
Copy link

@hmunozb hmunozb commented Mar 29, 2020

Enhances the 1-norm estimator to be generic over (f32, f64, c32, c64). Required to continue with SuperFluffy/rust-expm#3.

@SuperFluffy
Copy link
Owner

Hi, thanks for the PR and sorry it took me so long to respond. With the general craziness right now I am a bit tight on time. I promise I'll have a more thorough look in the next few weeks.

But a general question beside that: what do you base the extension to the complex numbers on? I admit, I don't have the details of normest1 and expm in mind, but are you certain that you can just extend the implementation like that?

@hmunozb
Copy link
Author

hmunozb commented Apr 14, 2020

The Higham and Tisseur paper had this to say about the complex case:

Everything in this section remains valid for complex matrices provided that sign(A) is redefined as the matrix (a_ij/|a_ij|) (and sign(0) = 1) and transposes are replaced by conjugate transposes. The matrix S is now complex with elements of unit modulus and we are much less likely to find parallel columns of S from one iteration to the next or within the current S. Therefore for complex matrices we omit the tests for parallel columns. However, we take the same, real, starting matrix.

as well as some additional justification for generalizing Z = A^T S to Z = A^* S in the complex case. So it was perfectly valid to just drop in an appropriately generic type, as the complex signum is precisely the complex phase as required, and replace transposition with hermitian conjugation. At the very least, this gives the advantage of making it straightforward to check that f32/f64 behavior is unchanged.

That said, this PR as-is right now is probably not without some sticking points. For instance, I opted not to omit the parallel column test in the complex case due to the special type check that it needed. I'm also regretting not marking the generic BLAS wrapper as unsafe, so I might do that and put back the unsafe blocks here in a sec.

Copy link

@emmatyping emmatyping left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks great! Just a couple of thoughts

lapacke = "0.2"
ndarray = "0.12"
cblas = "0.4"
lapacke = "0.5"

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd prefer to not require lapacke as well as lapack (it makes building on macOS much more difficult). Would you consider switching to using lapack? (It is mostly the same)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lapacke was the original requirement, so I just stuck to it. Upstream, maybe I can make lapacke and lapack options to the lapack-traits package, but it's probably up to @SuperFluffy which wrapper is used here.

I know linkage on macOS is a pain, but I got it to work out with a homebrew installed openblas backend. You just have to set the rustflags to point to the lib directories for openblas and GCC (for libgfortran). You can make your .cargo/config include something like this

[target.x86_64-apple-darwin]
rustflags = [  "-L", "native=/usr/local/opt/gcc/lib/gcc/11/",
               "-L", "native=/usr/local/opt/openblas/lib"
            ]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I unfortunately have requirements to statically link all my BLAS dependencies, as I know nothing of my target system other than "macOS", so it would complicate many things to link to homebrew.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see the problem. I've started an upstream issue for this: hmunozb/lapack-traits#1

src/lib.rs Outdated Show resolved Hide resolved
src/lib.rs Outdated Show resolved Hide resolved
{
let mut max_norm = 0.0;
//todo:

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the TODO here for?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

¯\_(ツ)_/¯

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants