newLISP compute matrix with lapacke and some problems

[kuan]

I want to compute a matrix with lapacke.

The example is here http://www.netlib.org/lapack/lapacke.html#_calling_tt_dgels_tt">http://www.netlib.org/lapack/lapacke.ht ... t_dgels_tt">http://www.netlib.org/lapack/lapacke.html#_calling_tt_dgels_tt



here is my code:

Code Select Expand

(set 'LAPACKE "/usr/local/lib/liblapacke.so")
(import LAPACKE "LAPACKE_dgels" "int" "int" "char" "int" "int" "int" "void*" "int" "void*" "int")
(set 'a '(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3))
(set 'b '(-10 -3 12 14 14 12 16 16 18 16))

(setq m 5 n 3 nrhs 2 lda 3 ldb 2)

(set 'aptr (pack (dup "lf" (* m n)) a))

(set 'bptr (pack (dup "lf" (* m nrhs)) b))

(set 'info (LAPACKE_dgels 101 78 m n nrhs (address aptr) lda (address bptr) ldb))


but the rsult is wrong.

Code Select Expand

(unpack (dup "lf" (* m nrhs)) bptr)
-> (1380.761318 1380.761318 1.844674407e+17 1.844674407e+17 1.844674407e+17 1.844674407e+17
 -1.415684734e+19 -1.415684734e+19 -1.153518066e+19 -1.153518066e+19)


I don't know where is wrong, can anyone give me some advice?

Thank you.

kuan.

[Lutz]

Two things:



In the LAPACKE C API documentation it says in the "Array Arguments" chapter: "Arrays are passed as pointers, not as a pointer to pointers.". The C-array used in the vars a and b is a "pointer to pointer", but they do not get passed as such, but get passed as simple pointers relying on the fact that all numbers will be contiguous in memory. I would try simply:

Code Select Expand


(LAPACKE_dgels 101 78 m n nrhs aptr lda bptr ldb)


and perhaps use:

Code Select Expand


(import LAPACKE "LAPACKE_dgels" "int" "int" "int" "int" "int" "int" "void*" "int" "void*" "int")

On the "int" versus "char" issue, I am not sure, and your guess is as good as mine. In their documentation they treat it as a "char" passing 'N', but I could not verify that the C call pattern really is "char". Unfortunately the docs do not contain the call patterns and the info they give in the code is incomplete. So, I would try "int" too. Perhaps you can find the exact call pattern for LAPACKE_dgels in the file lapacke.h.

[kuan]

I think "char" or "int" doesn't matter, when I change  "int" to/from "char", the result is the same as before.

And, when I do as you said, replace (address aptr) with aptr, strangely, the result doesn't change.

But when I test C array API with my own function below, the result is good. Here is code:

Code Select Expand

int sm(int m, int n, double * p)
{
    double *p_end;
    p_end = p + (m * n) - 1;
    for (;p <= p_end; p++)
        (*p) = (* p) + 1;
    return 0;
}

here is my test.lsp:

Code Select Expand

(set 'LIB "scalam.so")
(import LIB "sm" "int" "int" "int" "void*")
(set 'a '(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3))
(setq m 5 n 3)
(set 'ap (pack (dup "lf" (* m n)) a))
(set 'info (sm m n ap))
output is right

Code Select Expand

> (unpack (dup "lf" 15) ap)
(2 2 2 3 4 5 4 6 3 5 3 6 6 5 4)

My all questions are:

1.in newLISP, list is contiguous in memory? just as array? or anything else?

2.To compute with address or not doesn't matter, I am confused.

3.When I test my own sm function with b which has negative number, there are error with the negative elements:

Code Select Expand

>(set 'bp (pack (dup "lf" 10) b)
>(set 'info (sm 5 2 (address bp)))
> (unpack (dup "lf" 10) bp)
(1.844674407e+19 1.844674407e+19 13 15 15 13 17 17 19 17)
so, the reason is the negative number? and the more negative numbers, the more error

Code Select Expand

> (set 'b '(-10 -3 -12 -14 -14 12 16 16 18 16))
(-10 -3 -12 -14 -14 12 16 16 18 16)
> (set 'bp (pack (dup "lf" 10) b))
"000000000000�C000000000000�C000000000000�C000000000000�C000000000000�C000000000000(@0000000000000@0000000000000@0000000000002@0000000000000@"
> (set 'info (sm 5 2 (address bp)))
0
> (unpack (dup "lf" 10) bp)
(1.844674407e+19 1.844674407e+19 1.844674407e+19 1.844674407e+19 1.844674407e+19
13 17 17 19 17)
am I clear? Sorry for my terrible English.



PS: the dgels functioin in lapacke.h

Code Select Expand

lapack_int LAPACKE_dgels( int matrix_order, char trans, lapack_int m,
                          lapack_int n, lapack_int nrhs, double* a,
                          lapack_int lda, double* b, lapack_int ldb );

Code Select Expand

#ifndef lapack_int
#if defined(LAPACK_ILP64)
#define lapack_int              long
#else
#define lapack_int              int
#endif
#endif

#define LAPACK_ROW_MAJOR               101
#define LAPACK_COL_MAJOR               102

in dgels fortran source code, I found:

Code Select Expand

*> DGELS solves overdetermined or underdetermined real linear systems
*> involving an M-by-N matrix A, or its transpose, using a QR or LQ
*> factorization of A.  It is assumed that A has full rank.
*>
*> The following options are provided:
*>
*> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
*>    an overdetermined system, i.e., solve the least squares problem
*>                 minimize || B - A*X ||.
*>
*> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
*>    an underdetermined system A * X = B.
*>
*> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
*>    an undetermined system A**T * X = B.
*>
*> 4. If TRANS = 'T' and m < n:  find the least squares solution of
*>    an overdetermined system, i.e., solve the least squares problem
*>                 minimize || B - A**T * X ||.
*>

[Lutz]

I cannot see yet why it doesn't work using the LAPACK function, but can clarify the memory situation and make your C function work for negative values.



What you pass to the function is not a newLISP list but a contiguous piece of memory, put into a byte buffer of length 120 in ap. When passing that byte buffer of newLISP string-type to a C function the memory address of that buffer is passed. Applying the 'address' function to  a string, the same number is passed, so taking 'address' of string-type is redundant:

Code Select Expand


>
(set 'a '(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3))
(setq m 5 n 3)
(set 'ap (pack (dup "lf" (* m n)) a))

(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3)
3
"000000000000??000000000000??000000000000??00000000000000@00000000000008@00000000000016@00000000000008@00000000000020@00000000000000@00000000000016@00000000000000@00000000000020@00000000000020@00000000000016@00000000000008@"
> (string? ap)
true
> (address ap)
2028063120
> (unpack (dup "lf" (* m n)) (address ap))
(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3)
> (unpack (dup "lf" (* m n)) ap)
(1 1 1 2 3 4 3 5 2 4 2 5 5 4 3)
>


but in the following code, we have to use 'address':

Code Select Expand


> (get-float ap)
1
> (get-float (+ (address ap) 8))
1
> (get-float (+ (address ap) 16))
1
> (get-float (+ (address ap) 24))
2
> (get-float (+ (address ap) 32))
3
>

in this case we have to use 'address' to be able to do arithmetic on it.



As numbers are not in floating point format - no decimal point - convert them to floats, then negative numbers will work:

Code Select Expand


> (set 'a (map float '(-1 1 -1 2 -3 4 -3 5 -2 4 -2 5 -5 4 -3)))
(-1 1 -1 2 -3 4 -3 5 -2 4 -2 5 -5 4 -3)
> (setq m 5 n 3)
3
> (set 'ap (pack (dup "lf" (* m n)) a))
"000000000000?000000000000??000000000000?00000000000000@00000000000008?00000000000016@00000000000008?00000000000020@00000000000000?00000000000016@00000000000000?00000000000020@00000000000020?00000000000016@00000000000008?"
> (unpack (dup "lf" (* m n)) (address ap))
(-1 1 -1 2 -3 4 -3 5 -2 4 -2 5 -5 4 -3)
>


this would work too:

Code Select Expand

(set 'a '(-1.0 1.0 -1.0 2.0 -3.0 4.0 -3.0 5.0 -2.0 4.0 -2.0 5.0 -5.0 4.0 -3.0))

[kuan]

Thank you, lutz.

I transform "a" to float explicitly:

Code Select Expand


(set 'af (map float a))

then, the LAPACK_dgels works well.



"pack" looks like malloc or other allocate memory functions in C. It allocates continuous memory and assigns value of some type (float, int etc) to every position. This is so-called "packs expressions into a binary format". The return value of pack is string filled up with "00".



am I right?



PS:

Code Select Expand


> (file? "/usr/")
true
> (directory? "/usr/")
true

"file?" function "will also return true for directories". We've already had "directory?". This is why "directory?" function exits. When I want to test a path whether is a file, it can't give me an answer. Certainly, I can use "(not (directory? path))" , but it is not direct.

[Lutz]

there is a way to check for that in a future version:

http://www.newlisp.org/downloads/development/inprogress/">http://www.newlisp.org/downloads/develo ... nprogress/">http://www.newlisp.org/downloads/development/inprogress/