quil-lang / magicl Goto Github PK

Magicl loads various things through cffi, which itself marshals objects via libffi.dylib (on macOS). When run inside of my emacs/slime setup, this runs (mysterioulsy) smoothly, but when I save-lisp-and-die and reinitialize the Lisp image, sbcl's foreign library reinitialization does not remember the full pathname where libffi.dylib is located, and macOS's default LD_LIBRARY_PATH is empty. Setting LD_LIBRARY_PATH=/usr/lib before reopening the image fixes the issue. How can this be smoothed?

Here's an example run of a saved lisp image and its failure, so that onlookers can at least see where reinitialization happens.

debugger invoked on a SIMPLE-ERROR in thread
#<THREAD "main thread" RUNNING {10035A6A63}>:
  Error opening shared object "libffi.dylib":
  dlopen(libffi.dylib, 10): image not found.

Type HELP for debugger help, or (SB-EXT:EXIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
  0: [CONTINUE       ] Skip this shared object and continue.
  1: [RETRY          ] Retry loading this shared object.
  2: [CHANGE-PATHNAME] Specify a different pathname to load the shared object from.

(SB-SYS:DLOPEN-OR-LOSE #S(SB-ALIEN::SHARED-OBJECT :PATHNAME #P"libffi.dylib" :NAMESTRING "libffi.dylib" :HANDLE NIL :DONT-SAVE NIL))
0] backtrace

Backtrace for: #<SB-THREAD:THREAD "main thread" RUNNING {10035A6A63}>
0: (SB-SYS:DLOPEN-OR-LOSE #S(SB-ALIEN::SHARED-OBJECT :PATHNAME #P"libffi.dylib" :NAMESTRING "libffi.dylib" :HANDLE NIL :DONT-SAVE NIL))
1: (SB-ALIEN::TRY-REOPEN-SHARED-OBJECT #S(SB-ALIEN::SHARED-OBJECT :PATHNAME #P"libffi.dylib" :NAMESTRING "libffi.dylib" :HANDLE NIL :DONT-SAVE NIL))
2: (SB-SYS:REOPEN-SHARED-OBJECTS)
3: (SB-IMPL::FOREIGN-REINIT)
4: (SB-IMPL::REINIT)
5: ((FLET "WITHOUT-INTERRUPTS-BODY-10" :IN SAVE-LISP-AND-DIE))
6: ((LABELS SB-IMPL::RESTART-LISP :IN SAVE-LISP-AND-DIE))

0]

The OpenBLAS implementation seems to be all-round better (at least on macos), owing to concurrency and hardware-specific optimizations.

A (perhaps too) simple comparison:

With system (macos) default BLAS

(time (progn (reduce #'magicl:multiply-complex-matrices (mapcar (lambda (_) (quil::random-special-unitary 1000)) (list 1 2 3 4))) 1))

Evaluation took:
  53.799 seconds of real time
  53.642728 seconds of total run time (53.109868 user, 0.532860 system)
  [ Run times consist of 0.618 seconds GC time, and 53.025 seconds non-GC time. ]
  99.71% CPU
  145,903,290,806 processor cycles
  28 page faults
  5,055,004,624 bytes consed
  
1

With OpenBLAS

(time (progn (reduce #'magicl:multiply-complex-matrices (mapcar (lambda (_) (quil::random-special-unitary 1000)) (list 1 2 3 4))) 1))

Evaluation took:
  6.386 seconds of real time
  8.379127 seconds of total run time (7.713601 user, 0.665526 system)
  [ Run times consist of 0.482 seconds GC time, and 7.898 seconds non-GC time. ]
  131.21% CPU
  17,310,196,679 processor cycles
  5,078,847,024 bytes consed
  
1

MAGICL.LAPACK-CFFI:%ZUNCSD is broken, forcing the high-level library to define a working version

Request: define t and h as synonyms to (respectively) denote the transpose and Hermitian conjugate.

As a matter of fact, the numpy.matrix class has an H attribute that does this.

Originally posted by @jmbr in https://github.com/_render_node/MDExOlB1bGxSZXF1ZXN0MzE5MDI4MTYz/timeline/more_items

As @jmbr pointed out, these could be implemented with read tables as #t(matrix) and #h(matrix).

MAGICL currently has no matrix sum function—pretty goofy situation for a matrix library.

It looks like it might be possible to speed up the from-list function:

MAGICL> (let* ((n 50000)
                     (xs (loop :for i :below n :collect (coerce i 'double-float))))
                (sb-ext:gc :full t)
                (time (magicl:from-list xs (list n) :type 'double-float))
                (sb-ext:gc :full t)
                (time (magicl:from-array (coerce xs '(array double-float (*))) (list n)))
                (values))
Evaluation took:
  1.544 seconds of real time
  1.544000 seconds of total run time (1.544000 user, 0.000000 system)
  100.00% CPU
  4,188,179,784 processor cycles
  3,578,512 bytes consed
  
Evaluation took:
  0.000 seconds of real time
  0.004000 seconds of total run time (0.004000 user, 0.000000 system)
  100.00% CPU
  1,477,558 processor cycles
  400,016 bytes consed

We emit some really simple safety checks. Using the POLICY-COND library, we can check if SPEED > SAFETY, and only emit checks if that's not true.

Implement generic versions of some of the algorithms described here: https://arxiv.org/abs/physics/9806030 .

Just a suggestion (again, this PR covers a lot of ground, so I don't view this as a requirement or anything), it would be nifty to have an additional constructor which allows for one to generate a tensor by specifying index variables, corresponding dimensions, and a generating expression for the actual tensor entries. For example, I whipped together this macro (I think there are similar things for other multidimensional array libraries):

(defmacro generate-tensor (type index-exprs &body body-expr)
  "Generate an expression which evaluates to a MAGICL TENSOR of the given TYPE.

Here, INDEX-EXPRS is an expression of the form
   ((i0 d0) (i1 d1) ... (iN dN)).
Here D0, D1, ..., DN denote the shape of the resulting tensor. The variables
I0, I1, ... IN will be bound in BODY-EXPR, and with the value of BODY-EXPR
determining the value of the (I0, I1, ..., IN) entry of the tensor. "
  (alexandria:with-gensyms (shape tensor)
    (let ((vars (mapcar #'first index-exprs)))
      (labels ((generate-loops (i)
                 (if (= i (length vars))
                     `(setf (magicl:tref ,tensor ,@vars)
                            (coerce (progn ,@body-expr) ,type))
                     `(dotimes (,(elt vars i)
                                (elt ,shape ,i))
                        ,(generate-loops (1+ i))))))
        `(let* ((,shape (list ,@(mapcar #'second index-exprs)))
                (,tensor (magicl:const 0.0 ,shape :type ,type)))
           ,(generate-loops 0)
           ,tensor)))))

so that we can do something like

> (generate-tensor 'single-float ((i 3) (j 2)) 
     (+ i j))

#<MAGICL:MATRIX/SINGLE-FLOAT (3x2):
   0.000     1.000
   1.000     2.000
   2.000     3.000>

Originally posted by @kilimanjaro in https://github.com/_render_node/MDExOlB1bGxSZXF1ZXN0MzE5MDI4MTYz/timeline/more_items

This is a request for advice from the current maintainers; I plan on doing the necessary work to actually get this running on ABCL.

With abcl-1.7.0, I've been able to extend the implemention's creation of arrays specialized on commonly used byte types to use memory allocated outside of the JVM's heap via system interfaces like malloc().

To do this, ABCL has extended cl:make-array to take an additional keyword vis. :nio-buffer to use a generic contract for byte buffers for which I can use JNA to use an implementation that is backed by malloc()d memory

The version of CFFI distributed with the latest (2020-06-10) Quicklisp now contains an implementation of cffi-sys:make-shareable-byte-vector which utilizes this extension if present in the ABCL runtime.

As I understand things magicl relies on the ability of a given implementation to "pin" an already allocated array to a given memory location that is then used by the LAPACK et. al. libraries.

Unfortunately for ABCL, I can't "pin" an object that is already allocated on the JVM heap. The JVM GC needs to be able to move things around at will, unconditionally.

Therefore I conclude that I need to somehow "wrap" invocations of the cl:make-array used by magicl to get it to work on ABCL.

Possible paths:

1. Macroize cl:make-array. CONS: seems hairy; need to check all implementations

2. Require magicl consumers to only use arrays created cffi-sys:make-shareable-byte-vector. CONS: too much code would break to be practical as an initial implementation.

Please correct my misunderstandings. I've never used magicl: I just wanna run quilc programs on ABCL…

Use FIASCO for the framework.

All the error messages since #63 have been somewhat cryptic. This will likely also require a PR to stylewarning/policy-cond to allow for custom assertion messages in with-expectations.

For the wish list: generic support for matrix-vector multiplication, e.g.

(magicl:@ (magicl:from-list '(0d0 1d0 2d0 3d0) '(2 2)) 
          (magicl:from-list '(4d0 5d0) '(2)))

which currently gives an error, but works if the rightmost VECTOR is instead a MATRIX with shape `(2 1).

Originally posted by @kilimanjaro in https://github.com/rigetti/magicl/diffs/12

Just a thought: another nice bit of syntactic sugar is allowing at least one of the indices in the target shape to have unknown dimension, which is uniquely determined by the others when considering the size of the input tensor. For example, supposing that the shape argument can have one optional wildcard '*, then (reshape tensor '(* 1)) would be the same as (reshape tensor (list (size tensor) 1).

Originally posted by @kilimanjaro in https://github.com/_render_node/MDExOlB1bGxSZXF1ZXN0MzE5MDI4MTYz/timeline/more_items

When you do something like SAVE-LISP-AND-DIE, pointers to foreign data are not correctly saved.

It is annoying that re-generating the files causes 80k+ line commits.

1. Put them into a special directory.
2. Provide default files to compile if they haven't been user generated.

An offset argument should be added to FROM-DIAG and EYE to allow for setting a diagonal band not along the main diagonal.

From @jmbr

...the caller should be able to set which diagonal to fill (via a diagonal index as MATLAB and NumPy do). This is useful when creating tri-diagonal or banded matrices (e.g., when doing finite difference approximations of second order derivatives).

As LLA is now an abandonedware, I am looking into Magicl as a potential back-end of numcl. I would like to be briefed about the library in order to make myself 100% sure about the decision, especially the state of its high-level interface.

• I hope there is not much overlap between numcl and magicl yet --- numcl is still focusing on the high-level interface, and I believe magicl is still more on the low-level BLAS/LAPACK side. I only see einsum and couple of other things in the src/high-level/ directory, and hope not many libraries are depending on it yet.
• Could you brief me about the datastructure mainly passed around by magicl functions. numcl works on the natural CL simple 1D arrays in order to make future migration to C API easy.
• Handling of element-types. I know LLA handles element-type nicely, dispatches to the right functions and returns correct arrays. Is it handled by magicl already?
• Requirement to the pointer pinning and word alignment. Some Lisp allows to pin a pointer so that GC does not move them for some CFFI applications. Does magicl have any particular requirement regarding the raw pointers?

A pet peeve of mine in matrix packages is that they often do not properly handle matrices where one dimension is / both dimensions are set to zero. Since MAGICL already tracks the matrix dimension separately from the flat array of matrix entries, it should be relatively easy to include support for this edge case here.

Default thresholds should be x * epsilon where 1e1 <= x <= 1e3.

In many cases a sparse implementation of tensors will improve performance. Dense can remain the default implementation used but all abstract-tensor calls should be made implementation-agnostic.

This would likely require subclasses from abstract-tensor: sparse-tensor and dense-tensor.

The project is a bit of a mess and can be cleaned up. I propose something like this:

magicl/
    README, LICENSE, etc.
    .asd files
    src/
        high-level/
            high-level.lisp
            random.lisp
            ...
        bindings/
            computer-generated bindings files
        generate-interface/
            generate-interface.lisp
    test/
        test files

For some reason I'm not able to install magicl. Here's what I get in SBCL:

• (ql:quickload :magicl)

debugger invoked on a QUICKLISP-CLIENT:SYSTEM-NOT-FOUND in thread
#<THREAD "main thread" RUNNING {10005505B3}>:
System "magicl" not found

Type HELP for debugger help, or (SB-EXT:EXIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
0: [CONTINUE] Try again
1: [ABORT ] Give up on "magicl"
2: Exit debugger, returning to top level.

((LABELS QUICKLISP-CLIENT::RECURSE :IN QUICKLISP-CLIENT::COMPUTE-LOAD-STRATEGY) "magicl")
source: (CERROR "Try again" 'SYSTEM-NOT-FOUND :NAME NAME)
0]

What am I missing?

LAPACK often takes a parameter named info that is mutated during the call to a LAPACK routine to reflect for example whether the function was able to perform its calculation, and, if not, to then provide information on how to debug the issue. For example, the function dsegv uses this info parameter to signal whether a solution to some system of equations could be found (info = 0), and if not, info > 0 points to the relevent zero-valued diagonal.

In Lisp-world, we do provide the info parameter but it is passed-by-value and so cannot be mutated, leaving us without the diagnostic information to test whether e.g. MAGICL:SOLVE succeeded. The C bindings for LAPACK get around this same problem by instead returning the info value as the function return value and so there is no need to box the number as a mutable type. I don't know how that change would be made to the CL bindings that are generated, but I suspect that would be the best approach.

Beyond that, there are of course problem-specific ways of testing whether a LAPACK computation succeeded. For the example of solving linear equations, we could test the determinant or the condition number, or perhaps look for zeros on the diagonal. But that is probably duplicating effort / complicating things that the LAPACK routines will perform as a matter of course.

This was originally planned for #63 but needed to be pushed back.

Add support to create a "view" of a vector, allowing for free slicing and manipulation of sub-tensors. A suffix similar to the ! for mutation should be added to make it obvious when calls return a view.

This can be achieved with a "stride" field or similar on the tensor type or possibly defined at the abstract-tensor level.

The readme says to load magicl with (ql:quickload :magicl), but magicl is not available as a public quicklisp project. After adding it as a local project, the quickload operation was successful.

When running
CL-USER> (magicl:from-array (make-array '(3) :initial-contents '(3d0 3d0 3d0)) '(3))
I get an error of
no compatible tensor constructor for type T

This is because FROM-ARRAY by default specializes the tensor on the ARRAY-ELEMENT-TYPE of the input array when TYPE is not specified. This should be changed to inspect the first element in the case that ARRAY-ELEMENT-TYPE is T.

Right now, high-level.lisp deals with these relatively raw BLAS-inspired objects called matrix. Most of high-level.lisp just works with complex double floats. I think going the (old) MATLISP route of making nice object-oriented shims over all this business would be good.

Expokit doesn't like to exponentiate the zero matrix. This is pretty silly, especially since we can just check if a matrix is identically zero and return an identity matrix of the same size.

Expokit is enough of a nuisance to install and use, so it may be wise to separate it off into its own system, magicl.expokit which can be loaded optionally.

Loading magicl.expokit can still install the matrix functions into the expokit package.

The current implementation acts up when non-square matrices are passed in. These methods should allow for non-square matrices

There are no entries in the form cffi:define-foreign-library for libmkl_rt.dylib in :darwin, only for libmkl_rt.so on :unix.

Generate these all in one go, instead of spliced throughout the file.

Suggested by asjackson. https://developer.nvidia.com/cublas

This will allow einsum to be super efficient

The lists of matrix and vector classes should be stored in a central registry for introspection right at the time when they are declared (i.e., within defmatrix and defvector).

Originally posted by @jmbr in #92

No one has tested MAGICL with MKL in its different incarnations. Make sure everything works, or stop advertising support for this backend.

magicl faided in CCL64, OSX 10.14
CCL shows:

Error: Subprocess #<UIOP/LAUNCH-PROGRAM::PROCESS-INFO #x3020056B4CFD>
with command ("gfortran" "-fPIC" "-c" "/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f" "-o" "/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.o")
exited with error code 71
While executing: UIOP/RUN-PROGRAM::%CHECK-RESULT, in process Listener(4).
Type cmd-/ to continue, cmd-. to abort, cmd-\ for a list of available restarts.
If continued: IGNORE-ERROR-STATUS
Type :? for other options.

So I try this

gfortran -fPIC -c /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f -o /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.o

and I get

gfortran -fPIC -c /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f -o /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.o
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:294:72:

      call DNCHBV(mx,sgn*t_step,wsp(ih),mh,wsp(iexph),wsp(ifree+mx))
                                                                    1

Warning: Type mismatch in argument 'wsp' at (1); passed REAL(8) to COMPLEX(8) [-Wargument-mismatch] /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1580:15:

   do ip = 1,ndeg
                                                                    2
      theta(ndeg+ip) = CONJG( theta(ip) )
           1

Warning: Array reference at (1) out of bounds (8 > 7) in loop beginning at (2) /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1580:15:

   do ip = 1,ndeg
                                                                    2
      theta(ndeg+ip) = CONJG( theta(ip) )
           1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2) /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1581:15: /Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1579:72:

   do ip = 1,ndeg
                                                                    2

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1581:15:

      alpha(ndeg+ip) = CONJG( alpha(ip) )
           1

Warning: Array reference at (1) out of bounds (8 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1581:15:

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1579:72:

   do ip = 1,ndeg
                                                                    2

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1581:15:

      alpha(ndeg+ip) = CONJG( alpha(ip) )
           1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1591:15:

   do ip = 1,2*ndeg
                                                                    2
      alpha(ip) = 0.5d0*alpha(ip)
           1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1591:33:

   do ip = 1,2*ndeg
                                                                    2
      alpha(ip) = 0.5d0*alpha(ip)
                             1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1598:60:

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1590:72:

   do ip = 1,2*ndeg
                                                                    2

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1598:60:

         wsp(ih+(j-1)*m+j-1) = wsp(ih+(j-1)*m+j-1)-theta(ip)
                                                        1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1625:32:

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1590:72:

   do ip = 1,2*ndeg
                                                                    2

/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1625:32:

         y(j) = y(j) + alpha(ip)*wsp(iy+j-1)
                            1

Warning: Array reference at (1) out of bounds (14 > 7) in loop beginning at (2)
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:1901:72:

      call DNCHBV(mx,sgn*t_step,wsp(ih),mh,wsp(iexph),wsp(ifree+mx))
                                                                    1

Warning: Type mismatch in argument 'wsp' at (1); passed REAL(8) to COMPLEX(8) [-Wargument-mismatch]
/Users/user1/quicklisp/local-projects/magicl/transcendental/expokit.f:2283:72:

      call DNCHBV(mx,sgn*t_step,wsp(ih),mh,wsp(iexph),wsp(ifree+mx))
                                                                    1

Warning: Type mismatch in argument 'wsp' at (1); passed REAL(8) to COMPLEX(8) [-Wargument-mismatch]

any advice?

A mult method on these specialized matrices would speed up matrix arithmetic.

e.g.

(magicl:from-array (make-array '(3 3 3) :initial-element 0d0 :element-type '(complex double-float)) '(3 3 3))

The number of dimensions not equal to rank of array.
   [Condition of type SIMPLE-ERROR]

There are two sorts of promotion I have in mind:

• "broadcasting" tensors of one shape to a larger shape (e.g. to multiply a matrix by a vector)
• coercing tensor entries so that we can do stuff like multiply a MATRIX/DOUBLE-FLOAT by a MATRIX/COMPLEX-DOUBLE-FLOAT

The implementation will require a bit of strategy, but if done right it can be both convenient and avoid allocating new objects.

Current bindings produce defuns for which various LAPACK output values (e.g. mutated INFO variables) are lost to a Lisp caller. E.g.

(COMMON-LISP:DEFUN %DGETRF (M N A LDA IPIV INFO)
  (COMMON-LISP:DECLARE (COMMON-LISP:INLINE %%DGETRF)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) M)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) N)
                       (COMMON-LISP:TYPE
                        (COMMON-LISP:SIMPLE-ARRAY COMMON-LISP:DOUBLE-FLOAT) A)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) LDA)
                       (COMMON-LISP:TYPE
                        (COMMON-LISP:SIMPLE-ARRAY (COMMON-LISP:SIGNED-BYTE 32)
                         (COMMON-LISP:*))
                        IPIV)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) INFO))
  (CFFI:WITH-FOREIGN-OBJECTS ((M-REF6702 ':INT32) (N-REF6703 ':INT32)
                              (LDA-REF6705 ':INT32) (INFO-REF6707 ':INT32))
    (COMMON-LISP:SETF (CFFI:MEM-REF M-REF6702 :INT32) M)
    (COMMON-LISP:SETF (CFFI:MEM-REF N-REF6703 :INT32) N)
    (COMMON-LISP:SETF (CFFI:MEM-REF LDA-REF6705 :INT32) LDA)
    (COMMON-LISP:SETF (CFFI:MEM-REF INFO-REF6707 :INT32) INFO)
    (MAGICL.CFFI-TYPES:WITH-ARRAY-POINTERS ((A-REF6704 A) (IPIV-REF6706 IPIV))
      (%%DGETRF M-REF6702 N-REF6703 A-REF6704 LDA-REF6705 IPIV-REF6706
                INFO-REF6707))))

(from https://github.com/rigetti/magicl/blob/feature/magicl-ng/src/bindings/lapack02-cffi.lisp#L3725).
Mutations to INFO within the defun body will not be visible to callers.

It seems that any workaround is going to require rebuilding the lapack bindings.

• One option is to change them to explicitly take boxed values by default (say, a length 1 vector), and then users need to manage boxing/unboxing around the call site. (barf)
• Another option is to change the bindings so that the functions return the cffi:mem-ref values corresponding to the defun arguments. This is kind of heavyweight but shouldn't break existing code. Then callers can just check the last of the return values (corresponding to INFO).

There are probably other approaches.

Originally posted by @kilimanjaro in #63

The outcome of

(magicl.blas-cffi::%ddot 1 (make-array 1 :element-type 'double-float :initial-element 1.0d0) 1 (make-array 1 :element-type 'double-float :initial-element -1.0d0) 1)

is -1.0d0, as expected. However, the outcome of

(magicl.blas-cffi::%zdotu 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(1.0d0 0.0d0)) 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(-1.0d0 0.0d0)) 1)

is #C(0.0d0 0.0d0). The resulting zero value is consistent no matter which arrays are passed as input.

Something like (from-list '(1.0 2.0 3.0)) should suffice to instantiate a vector with 3 components. Compare with the current state of affairs: (from-list '(1.0 2.0 3.0) '(3)).

Nested lists with consistent sizes should also be transformed automatically (MAGICLly?) into tensors of the appropriate size.

The generated bindings file employ lots of CFFI magic with defcfun and the like. And the files are huge! lapack-cffi.lisp is 108497 lines long.

We ought to think of ways to improve compilation speed and memory usage, and (2) allow for a more streamlined binding set.