CLISP - an ANSI Common Lisp

Hi again,

I sent plenty of comments in the other e-mail, so here's the summary =
and proposed solutions (all untested).

(defun solve (matrice y)
  "R=E9solution MATRICE * x =3D Y par lapack"
  (let ((n (array-dimension matrice 0)))
    (with-foreign-object (mat-alien `(c-array double-float ,(* n n)) =
matrice)
      (with-c-var (x-alien `(c-array double-float ,n) y)
        (with-foreign-object (ipiv-alien `(c-array int ,n))
	  (let ((info (dgesv
		       n
		       1 ; nrhs
		       mat-alien
		       n ; lda
		       ipiv-alien
		       (c-var-address x-alien)
		       n))) ; ldb
	    (if (=3D info 0) x-alien info)))))))
Just an example, I would actually not recommend mixing c-var and
foreign-object. Just stick to one, consistently.

To be used with
(ffi:def-call-out dgesv
  (:name "dgesv_")
  (:return-type nil)
  (:language :stdc)
  (:library "liblapack.so")
  (:arguments (n (c-ptr int))
	      (nrhs (c-ptr int))
	      (a c-pointer) ; conceptually :in-out (c-ptr double-float ...)
	      (lda (c-ptr int))
	      (ipiv c-pointer) ; conceptually :out (c-ptr double-float n)
	      (b c-pointer) ; conceptually :in-out (c-ptr double-float ...)
	      (ldb (c-ptr int))
	      (info (c-ptr int) :out)))


One may recognize that this, again, is an instance of passing variable =
length arrays to/from functions. Now that the foreign function =
constructor is easily available in clisp, we can do this:

(defun memofn (n) ; actual memoization not shown
  (parse-c-type
   `(c-function
     (:return-type nil)
     (:arguments (n (c-ptr int))
		 (nrhs (c-ptr int))
		 ; a is strictly :in-out, but we throw away the result here
		 (a (c-ptr (c-array double-float ,n ,n)) :in)
		 (lda (c-ptr int))
		 (ipiv (c-ptr (c-array double-float ,n)) :out)
		 ; b is strictly :in-out, but we throw away the result here
		 (b (c-ptr (c-array double-float ,n)))
		 (ldb (c-ptr int))
		 (info (c-ptr int) :out))
     (:language :stdc))))

(defun solve (matrice y)
  ;; This code is bad, because the result is dereferenced even in case =
of
  ;; error. In general, doing so could result in run-time conversion =
errors.
  (multiple-value-bind (ipiv info)
      (funcall
       (foreign-function #'dgesv (memofn (length y)))
       n 1 matrice n y n)
    (if (=3D info 0) ipiv)))


In order to not dereference :out variables upon failures, we have to do =
what clisp does internally by hand: stack-allocated variables are =
manipulated explicitly then, with-foreign resurfaces:

(defun memofn (n) ; actual memoization not shown
  (parse-c-type
   `(c-function
     (:return-type nil)
     (:arguments (n (c-ptr int))
		 (nrhs (c-ptr int))
		 ; a is strictly :in-out, but we throw away the result here
		 (a (c-ptr (c-array double-float ,n ,n)) :in)
		 (lda (c-ptr int))
		 (ipiv c-pointer) ; conceptually :out (c-ptr double-float n)
		 ; b is strictly :in-out, but we throw away the result here
		 (b (c-ptr (c-array double-float ,n)))
		 (ldb (c-ptr int))
		 (info (c-ptr int) :out))
     (:language :stdc))))
; this differs from the above for ipiv.

(defun solve (matrice y)
  (let ((n (length y)))
    (with-foreign-object (ipiv `(c-array double-float ,n))
      (let ((info (funcall (foreign-function
			    #'dgesv (memofn (length y)))
			   n 1 matrice n ipiv y n)))
	(if (=3D info 0) (foreign-value ipiv))))))

I expect this memoized function declaration stuff to beat any custom =
nesting of with-foreign-object any time -- if memoized. Please report.

One could also memoize the result of foreign-function, but that's less =
likely to work across saved images.

BTW, UFFI is very far from such conciseness. But then, it's just a =
matter of piling up macros.


Remembers me that some month ago, I talked about a macro to provide =
exactly such variable-length function definition ability, with :guards =
(for :output checks) and bells and whistles. Any volunteer?

Regards,
	J=F6rg H=F6hle

"Hoehle, Joerg-Cyril" <Joerg-Cyril.Hoehle@...> writes:

> Hi again,
>
> I sent plenty of comments in the other e-mail, so here's the summary and =
proposed solutions (all untested).
>
> (defun solve (matrice y)
>   "R=E9solution MATRICE * x =3D Y par lapack"
>   (let ((n (array-dimension matrice 0)))
>     (with-foreign-object (mat-alien `(c-array double-float ,(* n n)) matr=
ice)

It does not work, due to rank differences. The conversion works if we
write:
(with-foreign-object (mat-alien `(c-array double-float (,n ,n)) matrice)
Here mat-alien is a 2-dimensional C array, that means a pointer to a
pointer to a double. But 2-D arrays in Fortran are always 1-D behind
the scene. So this won't work either.

>       (with-c-var (x-alien `(c-array double-float ,n) y)
>         (with-foreign-object (ipiv-alien `(c-array int ,n))
> 	  (let ((info (dgesv
> 		       n
> 		       1 ; nrhs
> 		       mat-alien
> 		       n ; lda
> 		       ipiv-alien
> 		       (c-var-address x-alien)
> 		       n))) ; ldb

Nice, you do not need to pass the :out argument ! After verification,
yes, it is in the FFI documentation, but I did not believe the FFI was
so smart.

> 	    (if (=3D info 0) x-alien info)))))))

Actually, in most of the Lapack subroutines, an non-null info value is
not an error, and the returned values contain meaningful information
(reduction to triangular form, some eigenvalues computed...).=20
If there is an error, Lapack quits (hence my question about how to
intercept the Fortran STOP statement).

[...]

> One may recognize that this, again, is an instance of passing variable le=
ngth arrays to/from functions. Now that the foreign function constructor is=
 easily available in clisp, we can do this:
>
> (defun memofn (n) ; actual memoization not shown
>   (parse-c-type
>    `(c-function
>      (:return-type nil)
>      (:arguments (n (c-ptr int))
> 		 (nrhs (c-ptr int))
> 		 ; a is strictly :in-out, but we throw away the result here
> 		 (a (c-ptr (c-array double-float ,n ,n)) :in)

You mean
 (a (c-ptr (c-array double-float (,n ,n))) :in)

> 		 (lda (c-ptr int))
> 		 (ipiv (c-ptr (c-array double-float ,n)) :out)
> 		 ; b is strictly :in-out, but we throw away the result here
> 		 (b (c-ptr (c-array double-float ,n)))
> 		 (ldb (c-ptr int))
> 		 (info (c-ptr int) :out))
>      (:language :stdc))))
>
> (defun solve (matrice y)
>   ;; This code is bad, because the result is dereferenced even in case of
>   ;; error. In general, doing so could result in run-time conversion erro=
rs.
>   (multiple-value-bind (ipiv info)
>       (funcall
>        (foreign-function #'dgesv (memofn (length y)))
>        n 1 matrice n y n)
>     (if (=3D info 0) ipiv)))

So here also ipiv is automatically allocated ?

> In order to not dereference :out variables upon failures, we have to do w=
hat clisp does internally by hand: stack-allocated variables are manipulate=
d explicitly then, with-foreign resurfaces:
>
> (defun memofn (n) ; actual memoization not shown
>   (parse-c-type
>    `(c-function
>      (:return-type nil)
>      (:arguments (n (c-ptr int))
> 		 (nrhs (c-ptr int))
> 		 ; a is strictly :in-out, but we throw away the result here
> 		 (a (c-ptr (c-array double-float ,n ,n)) :in)
> 		 (lda (c-ptr int))
> 		 (ipiv c-pointer) ; conceptually :out (c-ptr double-float n)
> 		 ; b is strictly :in-out, but we throw away the result here
> 		 (b (c-ptr (c-array double-float ,n)))
> 		 (ldb (c-ptr int))
> 		 (info (c-ptr int) :out))
>      (:language :stdc))))
> ; this differs from the above for ipiv.
>
> (defun solve (matrice y)
>   (let ((n (length y)))
>     (with-foreign-object (ipiv `(c-array double-float ,n))
>       (let ((info (funcall (foreign-function
> 			    #'dgesv (memofn (length y)))
> 			   n 1 matrice n ipiv y n)))
> 	(if (=3D info 0) (foreign-value ipiv))))))

Nice, but see above my remark about info return values.
I guess you still need a def-call-out definition, as #'dgesv is
not known with the above code alone.

> I expect this memoized function declaration stuff to beat any custom nest=
ing of with-foreign-object any time -- if memoized. Please report.

If n varies from call to call, as is typically the case, it will not
help. I am wondering why in clisp foreign arrays must have fixed
dimensions ?

> One could also memoize the result of foreign-function, but that's less li=
kely to work across saved images.
>
> BTW, UFFI is very far from such conciseness. But then, it's just a matter=
 of piling up macros.
>
>
> Remembers me that some month ago, I talked about a macro to provide exact=
ly such variable-length function definition ability, with :guards (for :out=
put checks) and bells and whistles. Any volunteer?
>
> Regards,
> 	J=F6rg H=F6hle
>

Thank you a lot for your comments, my code has already improved much.

--=20

Bernard Urban

Hi,
[sorry if you already got this, there's a hopefully intermittent =
problem with the dumb MS-Exchange->SMTP gateway]

Bernard Urban wrote:
>I guess you still need a def-call-out definition, as #'dgesv is
>not known with the above code alone.

Sure, you need a "prototype" definition for #'dgesv, and that will be =
redefined (cast) dynamically via the FOREIGN-FUNCTION constructor.


>>     (with-foreign-object (mat-alien `(c-array double-float=20
>,(* n n)) matrice)

>It does not work, due to rank differences. The conversion works if we
>write:
>(with-foreign-object (mat-alien `(c-array double-float (,n=20
>,n)) matrice)
>Here mat-alien is a 2-dimensional C array, that means a pointer to a
>pointer to a double. But 2-D arrays in Fortran are always 1-D behind
>the scene. So this won't work either.

I don't understand what you mean. 2D arrays in CLISP (and cmucl) are =
not pointer rows of pointers to doubles, (like char*argv[]), they are =
n*m doubles in sequence.

Otherwise, the following would yield garbage:
(with-c-var (m '(c-array uint8(2 3)) #2a((1 2 3)(4 5 6)))
  (cast m '(c-array uint8 6)))
-> #(1 2 3 4 5 6)

BTW, I forgot to recommend you investigate ROW-MAJOR-AREF instead of
using nested loops across the matrix (as in your cmucl code).

Maybe you wanted to express that the Fortran matrix does not use =
row-major order like Lisp (and C)? Indeed, your nested dotimes code is =
not expressing row-major order.
>	    (dotimes (i n)
>	      (dotimes (j n)
>		(setf (deref mat-alien (+ i (* n j))) (aref matrice i j))))
Alternatively, you could transpose the matrix?
(There was some mention of array element ordering varying between C and =
FORTRAN in the URLs I visited, maybe I now understand).

An alternative could be a (:language :fortran) extension to the FFI, =
that would treat matrices that way. Don't count on this.

>Nice, you do not need to pass the :out argument !
Same with cmucl, BTW, IIRC (that was ten years ago).


>If there is an error, Lapack quits (hence my question about how to
>intercept the Fortran STOP statement).

I have no idea. The two URLs I mentioned do not document the API you =
have (they pass integers directly, not pointers to :in integers). You =
have to consult the API. Maybe there's a callback or some such that is =
called in exceptional situations.

Throwing out of foreign code is often dangerous (unreliable in general =
-- think about unwind-protect across language barriers), you have to =
know the details of the concerned languages. Possibly there's some =
callback that you could install.


>You mean
> (a (c-ptr (c-array double-float (,n ,n))) :in)
Thank you very much! Indeed, that's the syntax,
not (c-array type n1 n2 ... nm). Silly me!

>So here also ipiv is automatically allocated ?
Yes, allocated on the stack prior to call, and convert to a Lisp array =
upon exit, thanks to:
		 (ipiv (c-ptr (c-array double-float ,n)) :out)

BTW, in the example where you allocate it yourself, I wrote:
		 (ipiv c-pointer) ; conceptually :out (c-ptr double-float n)
You could also write (ipiv (c-pointer (c-array double-float n))) =
;conceptually :out
This will give you better documentation and type-checks (note it's =
(c-pointer xyz), not (c-ptr xyz)). But it's a minor point (and a little =
slower).

Regards,
	J=F6rg H=F6hle.

Hi,

Bernard Urban wrote:
>> I don't understand what you mean. 2D arrays in CLISP (and 
>cmucl) are not pointer rows of pointers to doubles, (like 
>char*argv[]), they are n*m doubles in sequence.
>Yes, but in C, non-static 2D arrays are like char*argv[],
What's the definition of a non-static 2D array in C?

>and you are converting to such a beast.
No, I'm not! [cf. Monty Python]

>In Fortran, arrays are organized as in Lisp.
Row-major first, are your sure? Your code does it the other way round.

>> Otherwise, the following would yield garbage:
>> (with-c-var (m '(c-array uint8(2 3)) #2a((1 2 3)(4 5 6)))
>>   (cast m '(c-array uint8 6)))
>> -> #(1 2 3 4 5 6)
>Are you sure cast does not revert from the char*argv[] representtaion
>to the original ? My experiences with such things did not work in my
>lapack tests. 
FFI:CAST does nothing but check that SIZEOF and ALIGNOF match.

I'm sure I'm talking about 6 consecutive uint8 in memory. No *argv[].
There's a huge difference between
 (c-array (c-ptr (c-array x n) m) and
 (c-array (c-array x n) m).

I'm talking solely about the latter (2nd):
(with-c-var(m '(c-array uint8 (2 3)) #2a((1 2 3)(4 5 6)))
  (cast m '(c-array (c-array uint8 3) 2)))
-> #(#(1 2 3) #(4 5 6))
And it's an array of this kind that your initial cmucl and clisp code is filling.

(with-c-var(m '(c-array uint8 (2 3)) #2a((1 2 3)(4 5 6)))
  (cast m '(c-array (c-array uint8 2) 3)))
-> #(#(1 2) #(3 4) #(5 6))

>By the way, this is a method to coerce the rank of matrices !
It's exactly like the CLHS example on ARRAY-ROW-MAJOR-INDEX, showing off a displaced array using a different rank (giving a linear view on a matrix).

Regards,
	Jorg Hohle

[Reposted to the list]
Bernard Urban wrote:
>Some comments about your recents emails:
>
>1) Stack overflow in clisp: to avoid the problem, I switched to a
>allocate-{deep,shallow} formulation. The similarities with cmucl are
>now striking, see code below.=20
>I can now solve linear systems of size 1000x1000 in a matter
>of seconds.=20
>I remember that I had also (rare) stack overflow problems in the
>past in clisp on other programs, but not in cmucl. I think clisp
>must have a lower default stack limit, but I never investigated that.

>2) Array ordering on Lisp, Fortran, and C
>In Fortran and Lisp, all multidimensional arrays are behind the scene
>one-dimensional, with some syntaxic sugar to access them
>"naturally". But in your example of #2a((1 2 3)(4 5 6), if in Lisp the
>1D version is  #(1 2 3 4 5 6), in Fortran it will be  #(1 4 2 5 3 6).
>I have used ARRAY-ROW-MAJOR-INDEX in the code below=20
>as you suggested, and the difference is very visible now:=20
>you must transpose the input matrix to feed Fortran.
>
>I said I was not able to have correct result when initializing foreign
>variable mat-alien with a 2D Lisp matrix: I missed the
>transpose step, so final result was obviously false.
>
>Which leads us the case of C array semantic. Consider the following
>C code:
>
>/* in memory, 6 consecutive elements */
>double test[2][3] =3D { 1, 2, 3, 4, 5, 6 };
>
>and the following code:
>
>/* not sure that elements are consecutives */=20
>double * test_dyn[2];
>int i, j;
>  for (i =3D 0; i < 2; ++i) {
>    test_dyn[i] =3D calloc(3, sizeof(double));
>    for (j =3D 0; j < 3; ++j) {
>      test_dyn[i][j] =3D 1 + j + 3*i;
>    }
>  }
>
>test_dyn is what I called "dynamic" C arrays. The first "constant"
>ones are of rare use in my experience. So it did not make sense to me
>that FFI supports them, I was sure before your remarks we were
>speaking about the "dynamic" variant. This explains also the
>insistance on having constant matrix dimensions in the FFI.

Ah, I see. Indeed, test_dyn[][] and similar-looking test[][] are =
confusing.


>3) I got your UFFI patches, I will make some tests, and inform you
>   about the results.
>
>4) I do not intend to use your memofn function in that FFI case,=20
>   as it probably allocates on stack.

Yes, :out uses the stack, like with-foreign-object.

>5) In the following code, I added also the generic Lapack eigenvalue
>   solver (in case you wonder: in f2c Fortran, length of strings are
>   passed as hidden arguments). Stack limitation is even more critical
>   for that.=20

[cmucl variant]
are you sure you need
(multiple-value-bind (ret info)
  (declare ignore ret)
even for a function with no/void return type declaration?
I thought (I would have bet) cmucl and clisp behave the same and return =
no value, thus the primary value would be that of the first :out or =
:in-out parameter (as in CLISP).

>		(setf (deref mat-alien (array-row-major-index matrice j i))
>		      (aref matrice i j))))
I'd suggest you add a commment to your code about this transposition =
issue, otherwise it looks more like a typo.

[clisp variant]
>	  (x-alien (allocate-deep 'double-float y :count n))

>	       (let ((resul (make-array `(,n)=20
Just use (make-array n ...) for 1 dimension.
>					:element-type 'double-float)))
>		 (dotimes (i n (values resul (=3D info 0)))
>		   (setf (aref resul i)=20
>			 (element (foreign-value x-alien) i))))))

Just use
(values (foreign-value x-alien) (=3D info 0))
CLISP knows about the foreign array and can convert for you. For =
1-dimension, there's no problem with transposition.


>;;; Recherche de valeurs et vecteurs propres
>(defun valpro (matrice)
>  "Recherche des valeurs propres d'une matrice r=E9elle avec=20
[not reviewed]

>Again, thanks a lot for your feedback, I begin to see the light !

Regards,
	J=F6rg H=F6hle.

Hi,

Bernard Urban wrote:
[please always reply to the list as well]
>I have done some tests with your UFFI for clisp. I have choosen
>something different than lapack, namely the PROJ library  (a 
>geographic projection
>package; there are 'proj' and 'proj-ps-doc' Debian packages).

Thanks a lot for your report. Actually, it's the first one I ever received (I hope I forgot nobody, sorry if I do)!

>There is a one-line patch to apply to uffi.lisp (file patch-uffi).
Ah I see, I took advantage of the fact that the :library argument is actually evaluated, you maybe have an older clisp where this was not the case.

The advantage of the original code is that the library paths can be different between compilation and load&play times. No fixed path gets compiled in.

>The most annoying thing here is the need of the C-wrapper, it lowers a
>lot the usefulness of UFFI. Notice that clisp native ffi can avoid it
>easily, but cmucl seems to need it anyway, as returning struct seems
>here problematic.

Browsing uffi-ized code, you'll see that a lot of packages use some custom stubs, e.g. cl-sdl, clsql etc. IIRC...

Regards,
	Jorg Hohle.