# Lapack **URL:** \ **Category:** General\ **Created:** [December 15, 2023, 8:17am UTC](https://swi-prolog.discourse.group/t/lapack/7070 "2023-12-15T08:17:42Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![mgondan1](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/mgondan1/32/803_2.png) [@mgondan1](https://swi-prolog.discourse.group/u/mgondan1)\ **Post date:** [December 15, 2023, 8:17am UTC](https://swi-prolog.discourse.group/t/lapack/7070/1 "2023-12-15T08:17:42Z")
Suppose I wanted to add a lapack binding for swipl. I don’t think such a thing exists, or does it? What would be a good strategy? I guess I would put it on top of either clapack or lapack++ (the c++ binding) to avoid having to interface with fortran. I am not sure, but I guess lapack++ would fit my needs even if some functions are missing. The matrix objects would be swipl “blobs” (?), with some accessor functions for row or column vectors etc. Of course, there already exists a binding to R, which I could extend for my needs. ---
**Author:** ![jan](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/jan/32/4_2.png) [@jan](https://swi-prolog.discourse.group/u/jan)\ **Post date:** [December 15, 2023, 9:17am UTC](https://swi-prolog.discourse.group/t/lapack/7070/2 "2023-12-15T09:17:09Z")
Roughly, I think there are three routes: - As you suggest, use it through R - Use it through the new Janus Python interface - Add a C(++) interface. The latter is surely the most elaborate. Indeed, you represent external objects as blobs, either by packing a pointer or the data itself. In the end this route will give the best performance. The Janus approach could be the most attractive. Probably almost no work and pretty ok performance. ---
**Author:** ![mgondan1](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/mgondan1/32/803_2.png) [@mgondan1](https://swi-prolog.discourse.group/u/mgondan1)\ **Post date:** [December 15, 2023, 9:53am UTC](https://swi-prolog.discourse.group/t/lapack/7070/3 "2023-12-15T09:53:19Z")
The annoying thing with R is that the type of the return values is so unpredictable. I have no experience with Janus, though, and I don’t really know python. I guess I’ll try a little proof-of-concept with C++ and see how much work it is in terms of code lines. ---
**Author:** ![jan](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/jan/32/4_2.png) [@jan](https://swi-prolog.discourse.group/u/jan)\ **Post date:** [December 15, 2023, 10:21am UTC](https://swi-prolog.discourse.group/t/lapack/7070/4 "2023-12-15T10:21:43Z")
> [@mgondan1](#): > > I have no experience with Janus, though, and I don’t really know python. You might be surprised how easy it is if there is a good Python binding. I’m not sure about that. ---
**Author:** ![emiruz](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/emiruz/32/4495_2.png) [@emiruz](https://swi-prolog.discourse.group/u/emiruz)\ **Post date:** [December 15, 2023, 12:57pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/5 "2023-12-15T12:57:38Z")
If you just need matrix ops, numpy is more or less the same as R matrix ops. All the usual slicing and indexing is supported. You can use numpy through Janus. ---
**Author:** ![Boris](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/boris/32/7486_2.png) [@Boris](https://swi-prolog.discourse.group/u/Boris)\ **Post date:** [December 15, 2023, 12:58pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/6 "2023-12-15T12:58:13Z")
I am almost certain LAPACK and BLAS are only transitive dependencies of Numpy or Scipy or similar Python libraries ---
**Author:** ![jan](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/jan/32/4_2.png) [@jan](https://swi-prolog.discourse.group/u/jan)\ **Post date:** [December 15, 2023, 1:25pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/8 "2023-12-15T13:25:30Z")
> [@emiruz](#): > > You can use numpy through Janus. I’m sure you can 🙂 Do you have to have nice examples to share, for example in your blog? ---
**Author:** ![emiruz](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/emiruz/32/4495_2.png) [@emiruz](https://swi-prolog.discourse.group/u/emiruz)\ **Post date:** [December 15, 2023, 1:39pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/9 "2023-12-15T13:39:16Z")
I’ve nothing handy but I’ll rustle something up for you in the next few days ☺ ---
**Author:** ![peter.ludemann](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/peter.ludemann/32/48_2.png) [@peter.ludemann](https://swi-prolog.discourse.group/u/peter.ludemann)\ **Post date:** [December 15, 2023, 8:16pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/10 "2023-12-15T20:16:08Z")
> [@jan](#): > > Add a C(++) interface IMHO, the C++ interface (`SWI-cpp2.h`) is a lot easier than the C interface (the C interface requires correct code to do a lot of error checking; this is done automatically by the C++ interface). There’s also a new `PlBlob` class that simplifies implementing blobs. But there’s also a bit of a learning curve, and the documentation could use some improvement. ☹ > [@emiruz](#): > > You can use numpy through Janus. There’s also scipy, which seems to have an interface to LAPACK. ---
**Author:** ![mgondan1](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/mgondan1/32/803_2.png) [@mgondan1](https://swi-prolog.discourse.group/u/mgondan1)\ **Post date:** [December 15, 2023, 9:28pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/11 "2023-12-15T21:28:07Z")
I tried the PlBlob sample code, ```prolog #include PREDICATE(add, 3) { return A3.unify_integer(A1.as_long() + A2.as_long()); } struct MyBlob; static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob"); struct MyBlob : public PlBlob { std::unique_ptr connection; std::string name_; // Used for error terms ``` etc., but I get an error when compiling it (ubuntu): ```prolog matthias@rechner:~/plblas$ ~/swipl-devel/build/src/swipl-ld -c plblas.cpp plblas.cpp:10:53: error: expected primary-expression before ‘,’ token 10 | static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob"); | ^ plblas.cpp:10:28: error: ‘PL_BLOB_DEFINITION’ was not declared in this scope; did you mean ‘PL_BLOB_VERSION’? 10 | static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob"); | ^ ~~~~~~~~~~~~~~~~~ | PL_BLOB_VERSION plblas.cpp:13:1: error: expected class-name before ‘{’ token 13 | { std::unique_ptr connection; | ^ ``` Am I missing some include file? ---
**Author:** ![peter.ludemann](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/peter.ludemann/32/48_2.png) [@peter.ludemann](https://swi-prolog.discourse.group/u/peter.ludemann)\ **Post date:** [December 15, 2023, 10:50pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/12 "2023-12-15T22:50:51Z")
> [@mgondan1](#): > > tried the PlBlob sample code, I did a copy&paste from `test_cpp.cpp` and it compiled fine. Perhaps you’re missing the definition of `MyConnection`? (It appears that the documentation is slightly behind the test code; I’ll update the documentation) Here’s what worked for me, using the command `swipl-ld`: > **sample blob code** > > ```cpp > #include > > PREDICATE(add, 3) > { return A3.unify_integer(A1.as_long() + A2.as_long()); > } > > struct MyConnection > { std::string name; > > explicit MyConnection() { } > explicit MyConnection(const std::string& _name) > : name(_name) { } > > ~MyConnection() { } > > bool open() > { if ( name == "FAIL" ) // for testing error handling > return false; > return true; > } > > bool close() noexcept > { if ( name == "FAIL_close" ) // for testing error handling > return false; > return true; > } > }; > > // The following code is taken from > // pl2cpp2.doc \subsubsection{Sample PlBlob code} > // with some minor changes for testing > > struct MyBlob; > > static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob"); > > struct MyBlob : public PlBlob > { std::unique_ptr connection; > std::string name_; // Used for error terms > > explicit MyBlob() > : PlBlob(&my_blob) { } > > explicit MyBlob(const std::string& connection_name) > : PlBlob(&my_blob), > connection(std::make_unique(connection_name)), > name_(connection_name) > { if ( !connection->open() ) > throw MyBlobError("my_blob_open_error"); > } > > PL_BLOB_SIZE > > ~MyBlob() noexcept > { if ( !close() ) > Sdprintf("***ERROR: Close MyBlob failed: %s\n", name_.c_str()); // Can't use PL_warning() > } > > bool close() noexcept > { if ( !connection ) > return true; > bool rc = connection->close(); > connection.reset(); // Can be omitted, leaving deletion to ~MyBlob() > return rc; > } > > PlException MyBlobError(const char* error) const > { return PlGeneralError(PlCompound(error, PlTermv(symbol_term()))); > } > > int compare_fields(const PlBlob* _b_data) const override > { auto b_data = static_cast(_b_data); // See note about cast > return name_.compare(b_data->name_); > } > > bool write_fields(IOSTREAM *s, int flags) const override > { PlStream strm(s); > > strm.printf(",name=%s", name_.c_str()); > if ( !connection ) > strm.printf(",closed"); > return true; > } > }; > > // %! create_my_blob(+Name: atom, -MyBlob) is semidet. > PREDICATE(create_my_blob, 2) > { // Allocating the blob uses std::unique_ptr so that it'll be > // deleted if an error happens - the auto-deletion is disabled by > // ref.release() before returning success. > > auto ref = std::unique_ptr(new MyBlob(A1.as_atom().as_string())); > return A2.unify_blob(&ref); > } > > // %! close_my_blob(+MyBlob) is det. > // % Close the connection, silently succeeding if is already > // % closed; throw an exception if something goes wrong. > PREDICATE(close_my_blob, 1) > { auto ref = PlBlobV::cast_ex(A1, my_blob); > if ( !ref->close() ) > throw ref->MyBlobError("my_blob_close_error"); > return true; > } > > int main(void) { return 0; } > > ``` ---
**Author:** ![mgondan1](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/mgondan1/32/803_2.png) [@mgondan1](https://swi-prolog.discourse.group/u/mgondan1)\ **Post date:** [December 16, 2023, 7:58am UTC](https://swi-prolog.discourse.group/t/lapack/7070/13 "2023-12-16T07:58:45Z")
Problem solved. I think PL\_BLOB\_DEFINITION appeared lately, it is not yet in 9.0.4. For the records: I compile the current swipl-devel (in subfolder build, say). I don’t install it yet. A working swipl is found in swipl-devel/build/src/swipl, also swipl-ld. The latter does not seem to adjust the include path properly. Edit: it uses the /usr/bin/swipl found in the path. Probably a feature, not necessarily a bug. ```prolog matthias@rechner:~/plblas$ ~/swipl-devel/build/src/swipl-ld -c plblas.cpp -v eval `swipl --dump-runtime-variables` PLBASE="/usr/lib/swi-prolog" <<<<<<<<< PLARCH="x86_64-linux" PLSOEXT="so" PLLIBDIR="/usr/lib/swi-prolog/lib/x86_64-linux" PLLIB="-lswipl" PLTHREADS="yes" /usr/bin/c++ -c -D_REENTRANT -D __SWI_PROLOG__ -D __SWI_EMBEDDED__ ***-I/usr/lib/swi-prolog/include*** -o plblas.o plblas.cpp ``` ---
**Author:** ![peter.ludemann](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/peter.ludemann/32/48_2.png) [@peter.ludemann](https://swi-prolog.discourse.group/u/peter.ludemann)\ **Post date:** [December 16, 2023, 10:30am UTC](https://swi-prolog.discourse.group/t/lapack/7070/14 "2023-12-16T10:30:02Z")
You need a fairly recent “development” version. What do you get from this? `swipl --version` ---
**Author:** ![emiruz](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/emiruz/32/4495_2.png) [@emiruz](https://swi-prolog.discourse.group/u/emiruz)\ **Post date:** [December 16, 2023, 9:25pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/15 "2023-12-16T21:25:00Z")
Assuming you have a python3 environment, get numpy by typing `pip3 install numpy` at the command line, or `python3 -m pip install numpy` or perhaps `py -m pip install numpy` if under Windows. Here is some code: ```prolog % Doesn't try to convert complex Python objects to Prolog. raw(Cmd,Out) :- py_call(Cmd,Out,[py_object(true)]). % Tries to convert Python objects to Prolog. cmd(Cmd,Out) :- py_call(Cmd,Out). % Multiply two vectors. example1 :- X = [1,2,3,4], Y = [2,3,4,5], raw(numpy:array(X),V1), raw(numpy:array(Y),V2), raw(numpy:multiply(V1,V2),V3), cmd(numpy:ndarray:tolist(V3), Out), writeln(Out). % Square a matrix. example2 :- X = [[1,2,3], [3,5,3], [0,2,1]], raw(numpy:array(X),M), raw(numpy:matmul(X,X),Square), cmd(numpy:ndarray:tolist(Square),Out), writeln(Out). ``` It works as advertised by I’d note that getting something wrong (e.g. a typo) in py\_call results in the REPL going into an infinite loop, but I’m sure that will be ironed out in due course. There is an API reference for the vast array of function calls available to you here: [https://numpy.org/doc/stable/reference/index.html](https://numpy.org/doc/stable/reference/index.html) If you want to put together a specific example I may be able to help! NB: _The NumPy linear algebra functions rely on BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms. Those libraries may be provided by NumPy itself using C versions of a subset of their reference implementations but, when possible, highly optimized libraries that take advantage of specialized processor functionality are preferred. Examples of such libraries are [OpenBLAS](https://www.openblas.net/), MKL ™, and ATLAS. Because those libraries are multithreaded and processor dependent, environmental variables and external packages such as [threadpoolctl](https://github.com/joblib/threadpoolctl) may be needed to control the number of threads or specify the processor architecture._ - [Linear algebra (numpy.linalg) — NumPy v1.26 Manual](https://numpy.org/doc/stable/reference/routines.linalg.html) ---
**Author:** ![mgondan1](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/mgondan1/32/803_2.png) [@mgondan1](https://swi-prolog.discourse.group/u/mgondan1)\ **Post date:** [December 17, 2023, 8:18am UTC](https://swi-prolog.discourse.group/t/lapack/7070/16 "2023-12-17T08:18:20Z")
Impressive, indeed. I added a use\_module-directive, and removed a line in Example 2 that was copied from Example 1 (?). ```prolog use_module(library(janus)). % Does not try to convert complex Python objects to Prolog. raw(Cmd,Out) :- py_call(Cmd,Out,[py_object(true)]). % Tries to convert Python objects to Prolog. cmd(Cmd,Out) :- py_call(Cmd,Out). % Square a matrix. example2 :- X = [[1,2,3], [3,5,3], [0,2,1]], raw(numpy:transpose(X), XT), raw(numpy:matmul(XT, X),Square), cmd(numpy:ndarray:tolist(Square),Out), writeln(Out). ``` Can I get rid of the namespace numpy somehow? And is there a way to have two functions in the same expression? E.g. matmul(t(X), X) ---
**Author:** ![anon419083](https://avatars.discourse-cdn.com/v4/letter/a/53a042/32.png) [@anon419083](https://swi-prolog.discourse.group/u/anon419083)\ **Post date:** [December 17, 2023, 8:55am UTC](https://swi-prolog.discourse.group/t/lapack/7070/17 "2023-12-17T08:55:06Z")
> [@mgondan1](#): > > ```prolog > % Does not try to convert complex Python objects to Prolog. > raw(Cmd,Out) :- > py_call(Cmd,Out,[py_object(true)]). > % Tries to convert Python objects to Prolog. > cmd(Cmd,Out) :- py_call(Cmd,Out). > > ``` I had seen this syntax elsewhere using the janus interface, I think I can understand better what it means. The fact that “it tries” to convert Python objects means that the conversion may not succeed? And what one does get as an output in such cases? An error message? Thanks PS: I mean: the call that doesn’t try to convert should always succeed, right? ---
**Author:** ![jan](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/jan/32/4_2.png) [@jan](https://swi-prolog.discourse.group/u/jan)\ **Post date:** [December 17, 2023, 9:32am UTC](https://swi-prolog.discourse.group/t/lapack/7070/18 "2023-12-17T09:32:25Z")
> [@emiruz](#): > > Here is some code Thanks! > [@emiruz](#): > > It works as advertised by I’d note that getting something wrong (e.g. a typo) in py\_call results in the REPL going into an infinite loop, but I’m sure that will be ironed out in due course. Sure, I tried a few things though, but all works nicely. Any example? > [@mgondan1](#): > > Can I get rid of the namespace numpy somehow? And is there a way to have two functions in the same expression? E.g. matmul(t(X), X) Well, rather than raw/2 and cmd/2, you could define numpy/2, etc. You can pass the results of a Python function call into an argument using eval(Expr), e.g. ```prolog ex1 :- X = [1,2,3,4], Y = [2,3,4,5], py_call(numpy:ndarray:tolist(eval(numpy:multiply(eval(numpy:array(X)), eval(numpy:array(Y))))), Out), writeln(Out). ``` If you want to use this for serious Prolog based linear algebra you probably want to define a DSL that concisely expresses your intend and interpret or compile this to janus calls. Maybe you can implement [APL](https://en.wikipedia.org/wiki/APL_(programming_language)) 🙂 > [@anon419083](#): > > The fact that “it tries” to convert Python objects means that the conversion may not succeed? And what one does get as an output in such cases? An error message? Conversion from Python to Prolog always succeeds (except for resource errors or Python exceptions in case an iterator is translated into a list). Without py\_object(true), a number of Python objects are translated to natural counterparts in Prolog. Others produce, as with this option, a Python object reference. ---
**Author:** ![emiruz](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/emiruz/32/4495_2.png) [@emiruz](https://swi-prolog.discourse.group/u/emiruz)\ **Post date:** [December 17, 2023, 2:17pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/19 "2023-12-17T14:17:06Z")
> [@jan](#): > > Sure, I tried a few things though, but all works nicely. Any example? I’m running version 9.1.20-9-g4344213cf under emacs. Starts like this: ```prolog ?- py_call(zzz,_). % Interative session; added `.` to Python `sys.path` ERROR: Python 'ModuleNotFoundError': ERROR: No module named 'zzz' ERROR: In: ERROR: [12] janus:py_call(zzz,_156) ERROR: [11] toplevel_call('') at /home/me/lib/swipl/boot/toplevel.pl:1317 ^ Call: (6) forall(prolog:repl_loop_hook(end, _28), true) ? ``` I get an option regarding what to do, I pick “a” for abort. It then goes into an infinite loop, each segment of which looks like this: ```prolog ERROR: Arguments are not sufficiently instantiated ERROR: In: ERROR: [9] _1322>=0 ERROR: [8] prolog:repl_loop_hook(begin,_1356) at /home/me/lib/swipl/library/janus.pl:869 ERROR: [7] forall(prolog:repl_loop_hook(begin,_1394),'$toplevel':true) at /home/me/lib/swipl/boot/apply.pl:52 ERROR: [6] ERROR: [5] sig_atomic('$toplevel':forall(...,true)) ERROR: [4] setup_call_cleanup('$toplevel':forall(...,true),'$toplevel':'$query_loop'(0),'$toplevel':forall(...,true)) at /home/me/lib/swipl/boot/init.pl:679 ``` ---
**Author:** ![emiruz](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/emiruz/32/4495_2.png) [@emiruz](https://swi-prolog.discourse.group/u/emiruz)\ **Post date:** [December 17, 2023, 2:20pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/20 "2023-12-17T14:20:58Z")
> [@jan](#): > > Conversion from Python to Prolog always succeeds (except for resource errors or Python exceptions in case an iterator is translated into a list). Without py\_object(true), a number of Python objects are translated to natural counterparts in Prolog. Others produce, as with this option, a Python object reference. Here is an example: ```prolog ?- py_call(numpy:array([1,2,3]),X). % Interative session; added `.` to Python `sys.path` X = [(0x7f6963261950), (0x7f6963261d90), (0x7f69632618f0)]. ``` To me this seems like a partial conversion, since the values inside of the list are not Prolog types. Meanwhile, X can’t be put back into any numpy function which expect an ndarray structures because its been converted to a list, so its practical to receive references instead and then convert at the end. ---
**Author:** ![jan](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/jan/32/4_2.png) [@jan](https://swi-prolog.discourse.group/u/jan)\ **Post date:** [December 17, 2023, 4:02pm UTC](https://swi-prolog.discourse.group/t/lapack/7070/21 "2023-12-17T16:02:00Z")
> [@emiruz](#): > > ```prolog > ERROR: Arguments are not sufficiently instantiated > ERROR: In: > ERROR: [9] _1322>=0 > > ``` That has already been fixed shortly after it was introduced. Just bad luck 🙂 Upgrade and this should be fixed. > [@emiruz](#): > > `X = [(0x7f6963261950), (0x7f6963261d90), (0x7f69632618f0)].` I saw this. It is a rather unfortunate result of the default translation. The numbers are instances of `numpy.int64`, which subclasses eventually from `numpy.number`, but this class hierarchy seems unrelated to Pythons own `int` class. As a result, Janus does not translate to integers, which would be nice. As is, the result is rather useless: it is neither useful for Python nor for the Prolog programmer. Should we allow adding a method, e.g., ` __janus__ `() to any class that would then be used to obtain a value that is subsequently translated into a Prolog term? And, possibly, if ` __janus__ ()` returns `None`, return an object reference? That would allow us extend libraries to behave the way we want. Hmmm. It seems `numpy.int64` is an _immutable_ type and we cannot dynamically add methods to this. Some advice? [Next page](https://swi-prolog.discourse.group/t/lapack/7070.md?page=2)