BLAS bindings

Stream: lean4

[Tomas Skrivan]

I have created C BLAS bindings for Lean, LeanBLAS. Right now I have just level 1 for double precision.

Also I started writing specification in Spec directory, for level one I'm still missing iamax, rotg, rotmg and rot.

If anyone is interested in contributing, here is a list of things to do:

- add level one for single precision float
- add level one for complex numbers
- finish specification for level one
- add bindings for level 2 and 3
- write specification for level 2 and 3
- set up build for cuBLAS
- set up build for clBLAS

Copilot is really good at writing this kind of code. I would be interested if anyone would manage to oneshot all of this with some AI code generator.

[Eric Wieser]

Should ComplexFloat be in Lean core, since C99 has it natively?

[Eric Wieser]

Maybe the better question is "while it's not in Lean core, can we use implemented_by to compile it to C complex double?".

[Tomas Skrivan]

Good question, I'm planing to implement single precision float as

opaque Float32 := UInt32

Unfortunately, current compiler boxes types when they are marked opaque, so for some time it would have to be def Float32 := UInt32 and just not exploit defeq.

Similarly complex float can be done using UInt64 but for complex double I would need UInt128.

But it is not so bad as values in ComplexFloat are no boxed as they would be in Float × Float.

Hopefully, compiling ComplexFloat down to complex double will be possible as user code in the new complier.

[Henrik Böving]

Tomas Skrivan said:

Good question, I'm planing to implement single precision float as

opaque Float32 := UInt32

Unfortunately, current compiler boxes types when they are marked opaque, so for some time it would have to be def Float32 := UInt32 and just not exploit defeq.

Similarly complex float can be done using UInt64 but for complex double I would need UInt128.

But it is not so bad as values in ComplexFloat are no boxed as they would be in Float × Float.

Hopefully, compiling ComplexFloat down to complex double will be possible as user code in the new complier.

Leo added built-in support for Float32 this week

[Henrik Böving]

Eric Wieser said:

Maybe the better question is "while it's not in Lean core, can we use implemented_by to compile it to C complex double?".

That is not something that is technically possible with implemented_by, the set of types that Lean accepts as unboxed is currently not open to extension and hard controlled by the compiler.

[Eric Wieser]

Relatedly, perhaps mathlib's Complex should be generalized to Complex Real or similar (perhaps in batteries), so that you can get Complex Float64 for "free"

[Eric Wieser]

I don't know if defaulting the type argument to Real is well-behaved, but at least mathlib pretty much exclusively uses the notation anyway.

[Eric Wieser]

Plus then we can merge Complex and GaussianInt

[Tomas Skrivan]

Also you could use it for vector space complexification.

[Edward van de Meent]

Eric Wieser said:

Plus then we can merge Complex and GaussianInt

Then presumably we want to generalize docs#Zsqrtd to any type rather than just Int, and possibly upstream?

[Kevin Buzzard]

We should do it the Bourbaki way really, and have $R[X]/(X^2-aX-b)$ instead of just $R[X]/(X^2-d)$. This is the general free rank 2 $R$-algebra, you can't complete the square if 2 isn't invertible in $R$.

[Kevin Buzzard]

Maybe we have special notation or whatever for the case $a=0$ but we should really have the general thing too.

[Edward van de Meent]

that makes it noncomputable though

[Edward van de Meent]

(at least, i think?)

[Edward van de Meent]

there is docs#AdjoinRoot

[Edward van de Meent]

anyhow, the point is not to do fancy math things for Complex numbers, but to get a nice memory layout, and i don't believe going through Polynomial does that.

[Edward van de Meent]

Kevin Buzzard said:

We should do it the Bourbaki way really, and have $R[X]/(X^2-aX-b)$ instead of just $R[X]/(X^2-d)$.

As a small note, right now neither Complex nor GaussianInt is defined this way.

[Kim Morrison]

I think one can have both memory layout and mathematical generality. We'll define (apologies, on mobile)

structure Complex' (R : Type) (b : R := by exact -1) (a := by exact 0) where
  re : R
  im : R

The parameters a and b are then not used until it's time to define multiplication.

[Andrew Yang]

No the thing Kevin suggesting is not going through Polynomial but just generalize docs#Zsqrtd from pairs $n + mX$ such that $X^2 = d$ to $X^2=aX+b$.

[Kim Morrison]

The by exact is a trick to allow default values without typeclasses in the signature.

[Edward van de Meent]

i'm trying to define instances for the suggested structure... Just to be sure, is there a sensible instance of Star when a is nonzero? because simply negating the imaginary component is not it...

[Edward van de Meent]

i.e. does not yield a StarRing assuming the original ring is a CommRing

[Edward van de Meent]

(maybe at this point the conversation should be moved to either #mathlib4 or #maths?)

[Kim Morrison]

You need 2 to be invertible, so you can send X to the other root.

Doesn't X \mapsto -X + a work?

[Edward van de Meent]

i've run into a small issue, and i don't know what's causing it:

import Mathlib
set_option autoImplicit false

structure Complex' (R : Type*) (b : R := by exact -1) (a : R := by exact 0) where
  re : R
  im : R

namespace Complex'

def of {R : Type*} {b a : R} [Zero R] (r:R) : Complex' R b a where
  re := r
  im := 0

-- not relevant afaict
instance {R : Type*} {b a : R} : Mul (Complex' R b a) := sorry

-- issue in this declaration
theorem exists_coprime_of_gcd_pos {d:ℤ} {x : Complex' ℤ d} :
    ∃ (y : Complex' ℤ d), x = (of (x.re.gcd x.im : ℤ) : Complex' ℤ d) * y := by -- could not synthesize default value for parameter 'a' using tactics
  sorry

-- similar issue in this declaration, now at a different position
theorem exists_coprime_of_gcd_pos {d:ℤ} {x : Complex' ℤ d} : -- could not synthesize default value for parameter 'a' using tactics
    ∃ y, x = (of (x.re.gcd x.im : ℤ)) * y := by
  sorry


end Complex'

the error occurs at the uses of Complex' in the marked lines.
Fixing is easy enough, just provide the parameter explicitly... still, i'd love to know what triggers this issue...

[Kim Morrison]

I'm not sure. I wonder if there is some obstruction to running autoparams inside type ascriptions??

[Tomas Skrivan]

Does anyone have experience with BLAS on Windows? I'm really do not know how can I make the library easily portable :( Should I include source of OpenBLAS and just compile it on user machine? Should I download binaries?

cc @Mac Malone what is the recommended way for Lean project to depend on C library?

On Linux or Mac you just run sudo apt-get install libblas-dev or brew install openblas but on Windows you can't do anything simple like that.

[Tomas Skrivan]

Also I would like to provide optional dependency on CUDA/cuBLAS. What would be the recommended way of adding an ENABLE_CUDA switch? And how should downstream projects use this switch?