[GRASS-dev] gmath BLAS/LAPACK wrapper

["Sören Gebbert"]

Hi Brad, Dear devs,
i am currently working on a ATLAS - BLAS/LAPACK implementation
within the gpde library.
The idea is to provide full access to ATLAS BLAS level 1-3
functions and have multi-threaded BLAS 1-3 functions implemented in GRASS
as backup. I have almost finished the BLAS level 1 vector routines and
working on
level 2 the vector - matrix and level 3 matrix - matrix implementation.

The grass BLAS level 1-3 implementation is of course not as fast as
the ATLAS library.
ATLAS is cache optimized and uses recursive functions.
Additionally, ATLAS uses pthreads to provide multi threaded level 3 functions.

But for each specialised ATLAS function, a more general grass functions
will exists. Eg: there are level 3 functions for quadratic, symmetric and hermitian
matrices in BLAS, grass will provide only one matrix-matrix function.
If the user compiles grass without ATLAS support, all the modules
which are using
BLAS functions will still work ... but not that fast. ;)

I try to keep the interface as simple as i can. The vectors are 1d
float or double arrays
and matrices are 1d float or double arrays with additional row
pointers (using the G_alloc_vector and
G_alloc_matrix functions). I use the same low level implementation
as the ATLAS interface. The user must take care of correct allocation
and row, column counting.

Additionally i will try to create a wrapper to the LAPACK solver 
provided in ATLAS.
IMHO the gpde lu decomposition with row pivoting
and the gpde bandwidth optimized cholesky solver are the backup routines
for most of the LAPACK sover.

The krylov space solver (cg and bicgstab) will use the ATLAS BLASS 1-2
functions.

So the gpde interface is much more low level than the gmath LAPACK/BLAS
wrapper. And i think we can use the gpde routines instead of the native
lapack routines in the gmath warpper. But the matrix functions will work 
in row major order.
The blas functions are named like the (IMHO horrible)
netlib-blas naming convention, accept that i have put a N_ before the name.
Eg:
N_cblas_ddot, N_cblas_sdot ... .

If you use the gpde interface you will automatically use ATLAS. And if a user
don't compile grass with ATLAS support, 
the modules which are using the gmath structures
will still work. :)

What do you think? Unfortunately Glynn did not respond yet about this topic.
But i think renaming the gmath interface is a good idea. I had also in mind
to rename it in Gmath_ or G_math_ ... .

Here an ATLAS-gpde Interface example:

/*grass blas level 1 implementation of the dot product*/
double N_blas1_d_x_dot_y(double *x, double *y, int rows)
{
    int i;
    double s = 0.0;
#pragma omp parallel for schedule (static) reduction(+:s)
    for (i = rows - 1; i >= 0; i--) {
	s += x[i] * y[i];
    }
    return s;
}

/*The atlas wrapper with grass backup routine*/
double N_cblas_ddot(double *x, double *y, int rows)
{
#if defined(HAVE_ATLAS)
     return cblas_ddot(rows, x, 1, y, 1);
#else 
     return N_blas1_d_x_dot_y(x, y, rows)
#endif
}


Best regards
Soeren


Btw.:
I have updated the allocation routines in the gmath library but did not
committed it to cvs. I think i will have a look on your changes
before i commit.

-- 
Der GMX SmartSurfer hilft bis zu 70% Ihrer Onlinekosten zu sparen! 
Ideal für Modem und ISDN: http://www.gmx.net/de/go/smartsurfer