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