Is it possible to support complex dtypes on Windows?

[eriknw]

Is there any workaround we could do to make this work? Can we write any custom C code to serve as a bridge between Python and C GraphBLAS?

If anybody needs complex dtypes on Windows or has any information for how we might get it to work, please reply here!

[DrTimothyAldenDavis]

SuiteSparse:GraphBLAS should work with complex data types on Windows. Is this the python-to-C connection that doesn't work?

[eriknw]

Correct, we currently strip out all complex stuff for Windows. I wish we could get this working, since complex works with SuiteSparse:GraphBLAS on all OSes.

I don't think we need anything from you; we'll need to figure it out, or live with this shortcoming. Similarly, I wonder whether wrapping SuiteSparse:GraphBLAS with Cython would work with complex dtypes...

[DrTimothyAldenDavis]

It might be possible to change how I handle my 2 complex types. I could use my own struct, as in

typedef struct { double real ; double imag ; } GxB_FC64_t ;

I already do many of the complex operators myself (see for example https://github.com/DrTimothyAldenDavis/GraphBLAS/blob/0b79097dfba798c158af244eb7323a55adf9cf30/Source/GB_math.h#L25 and following).

[DrTimothyAldenDavis]

For example, I do my own complex division because the compiler-provided versions are mixed in quality. I want to match how MATLAB does it, and they are more careful about Inf and NaN behavior. So I use this:

https://github.com/DrTimothyAldenDavis/GraphBLAS/blob/0b79097dfba798c158af244eb7323a55adf9cf30/Source/GB_math.h#L336

and revised in my master branch as:

https://github.com/DrTimothyAldenDavis/GraphBLAS/blob/031fd110ba9e8c431ee5bb31f5609fc1aee00efa/Source/GB_math.h#L497

If I had an #ifdef'd option to use my own complex type definition, then the GxB_CMPLX(r,i) macro would be

#define GxB_CMPLX(r,i) (GxB_FC64_t) { .real=(r) ; .imag=(i) ; }

because this works:

#include <stdio.h>
#include <stddef.h>
typedef struct { double real ; double imag ; } mycomplex ;
int main (void)
{
    double pi = 3.14 ;
    double one = 1.0 ;
    mycomplex t = (mycomplex) { .real=pi, .imag=one } ;
    printf ("(%g,%g)\n", t.real, t.imag) ;
    #undef creal
    #define creal(x) ((x).real)
    #undef cimag
    #define cimag(x) ((x).imag)
    printf ("(%g,%g)\n", creal (t), cimag (t)) ;
}

If this was an #ifdef'd option, then GxB_FC32_t and GxB_FC64_t would both be typedef'd structs. Can you pass in scalars and arrays of that type?

[DrTimothyAldenDavis]

Downside to this: I would need to write my own complex trig and hyperbolic functions: cacos, cacosf, casin, casinf, catan, catanf, ..., cacosh, ... cpow, cexp, clog, and so many more. That would be very hard to do. I would need to find a set of definitions of these functions that I could use and incorporate into GraphBLAS.

[DrTimothyAldenDavis]

What if instead you pretend on your side that GxB_FC64_t has this typedef:

typedef struct { double real ; double imag ; } GxB_FC64_t ;

and then whereever you need to pass in a C scalar of that type, pass in that struct. Or a C array of that type -- just pass in a pointer to an array of those types.

[eriknw]

What if instead you pretend on your side that GxB_FC64_t has this typedef:

typedef struct { double real ; double imag ; } GxB_FC64_t ;

and then whereever you need to pass in a C scalar of that type, pass in that struct. Or a C array of that type -- just pass in a pointer to an array of those types.

Yeah, we should be able to do that for Windows easily enough. If we tried this, do you think it should work with the latest SuiteSparse:GraphBLAS, or would you need to do something?

[DrTimothyAldenDavis]

I would need to make a 3-line change to GraphBLAS.h.in, and the change would have no effect on anything inside SS:GrB, nor any effect on any other users of GraphBLAS. Here's what I would do:

The 3 lines would be:

#ifndef GXB_COMPLEX_DEFN
#define GXB_COMPLEX_DEFN

... here I place lines 137 to 203 of GraphBLAS.h
but I would remove the definition of GB_restrict (already done in the master branch).

#endif

Then in your C wrapper, you do this:

#define GXB_COMPLEX_DEFN
typedef struct { double real ; double imag ; } GxB_FC64_t ;
typedef struct { float real ; float imag ; } GxB_FC32_t ;
#define GxB_CMPLX(r,i) (GxB_FC64_t) { .real=(r) ; .imag=(i) ; }
#define GxB_CMPLXF(r,i) (GxB_FC32_t) { .real=(r) ; .imag=(i) ; }
#include "GraphBLAS.h"

Then your notion of the two complex types would differ from mine, but they would be compatible at the bit level. Your C wrapper would think they are structs of 2 floats or 2 doubles. All of the GraphBLAS prototypes that you see on your side would use that definition. All of they prototypes GraphBLAS would see internally would be the built-in complex type in Windows, since you would still compile SS:GrB the same, without this preface of #define GXB_COMPLEX_DEFN.

This would be a minor change to GraphBLAS so I could do it as a v7.5.0 version, released any time you like. I'd bump the minor version from v7.4.x to v7.5.0 because it's a new user-visible feature. It would have no effect at all on any other GraphBLAS users, except to add the #define preprocessing token GXB_COMPLEX_DEFN.

Would that work on your side?

[alugowski]

In case this is helpful: here is how pybind11 does individual complex values: https://github.com/pybind/pybind11/blob/master/include/pybind11/complex.h
And their nparray code: https://github.com/pybind/pybind11/blob/master/include/pybind11/numpy.h
I don't see anything windows specific, nor anything weird about complex values.

I can confirm that pybind11 works for me with complex numpy arrays being read and written by a C++ extension, including on Windows.

[alugowski]

Oh, I think I see what's causing the problem. The MSVC command that compiles the Python extensions does not specify a C standard. The MSVC default is something ancient, which does not have complex value support.

I'll submit a PR shortly that will include a /std:c11 switch. Last few days I've gotten cibuildwheel working so that maybe we can get Python GraphBLAS working normally on all platforms. I've added that switch to my changes.

Tim, GraphBLAS/CMakeLists.txt does set the C standard, but in a confusing way. The line that does so is this one:

set_property ( TARGET graphblas PROPERTY C_STANDARD 11 )

(line 293). The confusing bits:

• there is also C_STANDARD_REQUIRED 11 which is not a correct use, since that is an ON/OFF value, not a number. It just says that the value of C_STANDARD is required, and should be set to "ON" not "11". See https://cmake.org/cmake/help/latest/prop_tgt/C_STANDARD_REQUIRED.html.
• The big compiler-specific set of IFs list C-standard switches such as -std=c11, which is redundant given the C_STANDARD property is set. Of course, all except MSVC :)

Also FYI, the preferred way to check for Clang is:

if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")

Notice the MATCHES, not STREQUAL. That way you handle both regular Clang and AppleClang, the default (and often only) compiler on macOS.

[alugowski]

The cibuildwheel PR with /std:c11 is ready: #66

I'll let someone familiar with the complex value issue see if that fixes it or not.

[DrTimothyAldenDavis]

Thanks for the comments. I'll work on fixing those issues in the GraphBLAS/CMakeLists.txt file. It was the first cmake I wrote and you can see it's a bit crufty in a few places -- in particular, the large if-else block for handing different compilers. Some of that can be cleaned up.

[alugowski]

This took far longer than I hoped, but I have some progress.

The /std:c11 flag isn't enough. MSVC still chokes on float _Complex as if it can't parse a two-token type. I don't know if that's a general MSVC bug (googling suggests this issue is common, but MSVC docs imply it should work) or if cffi and distutils interfere in some way.

However, using _Dcomplex instead of double _Complex does work. The problem is that cffi won't let you. So if you insist on cffi then you need hacks.

• There might be a "cleaner" hack to somehow convince cffi to accept _Dcomplex. Maybe there's some way to declare a type? Maybe declare a double_complex to something neutral then add a #define double_complex _Dcomplex in source.c? Someone familiar with cffi might have ideas.
• The dirtier hack is to just emit the code with ffibuilder.emit_c_code('_graphblas.c') then:
  code = code.replace("float _Complex", "_Fcomplex").replace("double _Complex", "_Dcomplex")
  and load that as a regular non-cffi extension.

I hacked an ugly version of approach 2 ( alugowski@dd70141 ) and:

Run coverage run --branch -m pytest
============================= test session starts =============================
platform win32 -- Python 3.9.16, pytest-7.2.1, pluggy-1.0.0
Using --randomly-seed=1119680806
rootdir: D:\a\python-suitesparse-graphblas\python-suitesparse-graphblas, configfile: setup.cfg, testpaths: suitesparse_graphblas/tests
plugins: randomly-3.12.0
collected 4 items / 1 skipped

suitesparse_graphblas\tests\test_scalar.py .                             [ 25%]
suitesparse_graphblas\tests\test_doctest.py .                            [ 50%]
suitesparse_graphblas\tests\test_exceptions.py .                         [ 75%]
suitesparse_graphblas\tests\test_package.py .                            [100%]

======================== 4 passed, 1 skipped in 0.34s =========================

Then I noticed those tests have nearly zero coverage so barely confirm anything, so I re-enabled test_io.py which at least touches complex values and got this:

=================================== ERRORS ====================================
___________ ERROR collecting suitesparse_graphblas/tests/test_io.py ___________
suitesparse_graphblas\tests\test_io.py:27: in <module>
    from suitesparse_graphblas.io import binary  # isort:skip
suitesparse_graphblas\io\binary.py:14: in <module>
    stdlib = stdffi.dlopen(find_library("c"))
C:\Miniconda\envs\suitesparse-graphblas\lib\site-packages\cffi\api.py:150: in dlopen
    lib, function_cache = _make_ffi_library(self, name, flags)
C:\Miniconda\envs\suitesparse-graphblas\lib\site-packages\cffi\api.py:832: in _make_ffi_library
    backendlib = _load_backend_lib(backend, libname, flags)
C:\Miniconda\envs\suitesparse-graphblas\lib\site-packages\cffi\api.py:821: in _load_backend_lib
    raise OSError("dlopen(None) cannot work on Windows for Python 3 "
E   OSError: dlopen(None) cannot work on Windows for Python 3 (see http://bugs.python.org/issue23606)
=========================== short test summary info ===========================
ERROR suitesparse_graphblas/tests/test_io.py - OSError: dlopen(None) cannot work on Windows for Python 3 (see http://bugs.python.org/issue23606)
!!!!!!!!!!!!!!!!!!! Interrupted: 1 error during collection !!!!!!!!!!!!!!!!!!!!
============================== 1 error in 0.55s ===============================

See https://github.com/alugowski/python-suitesparse-graphblas/actions/runs/4140578467/jobs/7159365826

That looks like a different issue, though, so maybe this is a success?

[alugowski]

The cffi docs say code from emit_c_code() is intended to be distributed; it's not specific to the cffi that created it. So it's feasible to convert build.py to a generator akin to create_headers.py. Then ship a _graphblas.c and a _graphblas_msvc.c and build suitesparse_graphblas._graphblas as a regular extension (choosing which .c to use based on platform).

Ugly, but it'll work.

[alugowski]

The cffi docs say code from emit_c_code() is intended to be distributed; it's not specific to the cffi that created it. So it's feasible to convert build.py to a generator akin to create_headers.py. Then ship a _graphblas.c and a _graphblas_msvc.c and build suitesparse_graphblas._graphblas as a regular extension (choosing which .c to use based on platform).

Ugly, but it'll work.

Or just generate those files at setup time instead of checking them in. That would be nicer; those generated .c files have 81k lines.

[eriknw]

Interesting. This does sound like progress. So on Windows we would emit _graphblas.c, modify it, then compile it to _graphblas.dll? That doesn't sound too bad 🤞, but I want this to happen automatically when building. These files get created (then deleted) during the build process anyway.

I feel like I've seen that dlopen(None) error before. I'll let you know if I remember/find anything.

We can test things more thoroughly in python-graphblas. In our CI, we have jobs that install python-suitesparse-graphblas from the latest main branch, and we can easily point to a specific commit or PR. We test installing from wheels and conda-forge too.

[eriknw]

suitesparse_graphblas\tests\test_scalar.py . [ 25%]

Oh, hey, this test ran and passed! This test is called test_complex. So, looks like it worked!

so I re-enabled test_io.py

Okay, I'm more confident now that the exception you're seeing on Windows is expected and is why we skip those tests. This is a limitation of Windows we can live with.

[alugowski]

Sounds like this approach is acceptable in principle, so I'll clean it up and make a PR.