newlispfanclub.alh.net

Posted: **Mon Dec 17, 2012 6:44 am**

I want to compute a matrix with lapacke.
The example is here http://www.netlib.org/lapack/lapacke.ht ... t_dgels_tt

here is my code:

Code: Select all

(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 all

(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.

Posted: **Mon Dec 17, 2012 3:54 pm**

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 all

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

and perhaps use:

Code: Select all

(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.

Posted: **Mon Dec 17, 2012 6:12 pm**

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 all

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 all

(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 all

> (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 all

>(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 all

> (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))
"\000\000\000\000\000\000�C\000\000\000\000\000\000�C\000\000\000\000\000\000�C\000\000\000\000\000\000�C\000\000\000\000\000\000�C\000\000\000\000\000\000(@\000\000\000\000\000\0000@\000\000\000\000\000\0000@\000\000\000\000\000\0002@\000\000\000\000\000\0000@"
> (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 all

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 all

#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 all

*> 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 ||.
*>

Posted: **Mon Dec 17, 2012 8:58 pm**

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 all

> 
(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
"\000\000\000\000\000\000??\000\000\000\000\000\000??\000\000\000\000\000\000??\000\000\000\000\000\000\000@\000\000\000\000\000\000\008@\000\000\000\000\000\000\016@\000\000\000\000\000\000\008@\000\000\000\000\000\000\020@\000\000\000\000\000\000\000@\000\000\000\000\000\000\016@\000\000\000\000\000\000\000@\000\000\000\000\000\000\020@\000\000\000\000\000\000\020@\000\000\000\000\000\000\016@\000\000\000\000\000\000\008@"
> (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 all

> (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 all

> (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))
"\000\000\000\000\000\000?\000\000\000\000\000\000??\000\000\000\000\000\000?\000\000\000\000\000\000\000@\000\000\000\000\000\000\008?\000\000\000\000\000\000\016@\000\000\000\000\000\000\008?\000\000\000\000\000\000\020@\000\000\000\000\000\000\000?\000\000\000\000\000\000\016@\000\000\000\000\000\000\000?\000\000\000\000\000\000\020@\000\000\000\000\000\000\020?\000\000\000\000\000\000\016@\000\000\000\000\000\000\008?"
> (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 all

(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))

Posted: **Tue Dec 18, 2012 4:46 am**

Thank you, lutz.
I transform "a" to float explicitly:

Code: Select all

(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 "\000".

am I right?

PS:

Code: Select all

> (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.

Posted: **Tue Dec 18, 2012 8:12 pm**

there is a way to check for that in a future version:
http://www.newlisp.org/downloads/develo ... nprogress/

newlispfanclub.alh.net

newLISP compute matrix with lapacke and some problems

newLISP compute matrix with lapacke and some problems

Re: newLISP compute matrix with lapacke and some problems

Re: newLISP compute matrix with lapacke and some problems

Re: newLISP compute matrix with lapacke and some problems

Re: newLISP compute matrix with lapacke and some problems

Re: newLISP compute matrix with lapacke and some problems