[Axiom-developer] Re: BLAS and GCL

[Camm Maguire]

Greetings!

Tim Daly  <address@hidden> writes:

> In response to the numericIntegration question I looked at what
> would be required to rebuild the missing NAG functionality. The
> base layer seems to be BLAS (Basic Linear Algebra Subroutines).
> I downloaded BLAS over the weekend and rebuilt it as pamphlets.
> Now I'm going to try to link it into the GCL image. Have you 
> ever seen a fortran linkage or should I use the C cover functions
> (CBLAS)?
> 

Well, fortunately or not, I also happen to maintain blas/lapack/atlas
for Debian, and we use the blas heavily here at work.  I've always
intended to give GCL blas functionality -- we even discussed on the
mailing list how to instruct the compiler to optimize written out lisp
loops into blas calls.  On Debian, such calls are automatically
redirected at runtime to the most efficient blas implementation
installed for the running cpu, taking advantage of various subarch
extensions such as sse2, etc.  There is matlisp, which I have not
looked at, but my guess is that it would be far more overhead than GCL
needs, as if there is anything at which GCL excels, it is in its C
interface. 

My only uncertainty in how to proceed here lies in how to modularize
this capability.  One surely does not want all gcl images to depend on
a blas shared external lib.  On the other hand, one wants to be able
to load a module with the lisp interface, dlopen the lib, save the
system, and have the resulting image still depend on blas and to work
properly.  One can do this at present via (compiler::link
"lisp_blas.o" "new_image" "" "-lblas -lg2c"), but the beauty of
save-system over this is that it maintains the state of the heap
created in the session up to that point, wheras the link function
would have to be supplied an optional argument containing code to
recreate this state in the newly linked image.  However, the dynamic
linker ld.so will not redirect symbols in loaded compiled lisp objects
after a save-system and re-execution -- GCL's internal linker is
'static'.  I suppose the solution ultimately is to rewrite sfaslbfd.c
and the like to leave such symbols undefined but marked for processing
by ld.so -- should be possible, I just don't know how to do it yet.

In any case, this is a lot more than you wanted to know no doubt.  In
short, you can use either the fortran or C interface functions.  On
Debian, cblas wrappers are rolled into the same single libblas with
the fortran.

An example fortran call would look like:

extern void
dgemv_    (char *,int *,int *,double  *,
           const double  *,int *,const double  *,
           int *,double  *,double  *,int *);

dgemv_    (&t,&n,&m,&alpha,a,&lda,x,&ix,&beta,y,&iy);

i.e. all functions have an underscore postpended, all arguments are
passed by reference, and one can link with the normal C linker adding
-lg2c.  A home grown C wrapper for this would look like:

int
dgemv    (char t,int m,int n,double  alpha,const double  *a,
          int lda,const double  *x,int ix,
          double  beta,double  *y,int iy);

One way to proceed would therefore be to place the following in
dgemv.lisp (obviously more error checking is needed, but you get the
idea):

=============================================================================
(clines "
extern void
dgemv_    (char *,int *,int *,double  *,
           const double  *,int *,const double  *,
           int *,double  *,double  *,int *);

int
dgemv    (char t,int m,int n,double  alpha,object  a,
          int lda,object  x,int ix,
          double  beta,object  y,int iy) {

    t=t=='T' ? 'N' : 'T'; /* Fortran arrays are column-major*/
    dgemv_    
(&t,&n,&m,&alpha,a->lfa.lfa_self,&lda,x->lfa.lfa_self,&ix,&beta,y->lfa.lfa_self,&iy);
  return 0;

}
init_dgemv() {}

")


(defentry %dgemv (char int int double object int object int double object int) 
(int dgemv))
=============================================================================

and then in gcl:

>(compile-file "/tmp/dgemv.lisp")
>(compiler::link '("/tmp/dgemv.o") "/tmp/new_image" "" "-lblas -lg2c")


One could similarly use the cblas interface if desired.

Come to think of it, we should really have compiler::link append the
new libs to compiler::*ld-libs* for subsequent invocations in the new
image.

Take care,

> Tim
> 
> 
> 

-- 
Camm Maguire                                            address@hidden
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah