Hello,
After some tests, for the intended use (multiply a matrix by a
scalar), dgemm is not faster
that dscal, but in the C code of "iMultiRealScalarByRealMatrix",
the part which takes the most
of the CPU time is the call to "dcopy". For example, on my
machine, for a 10000x10000 matrix,
the call to dcopy takes 540 milliseconds and the call to dscal 193
milliseconds. Continuing my
explorations today, I tried to see how Scilab expressions such as
Y+2*X
are parsed and executed. To this purpose I have written an
interface (daxpy.sci and daxpy.c
attached) to the BLAS function "daxpy" which does "y<-y+a*x" and a
script comparing the above
expression to
daxpy(2,X,Y)
for two 10000x10000 matrices. Here are the results (MacBook air
core i7@1,7GHz):
daxpy(2,X,Y)
(dcopy: 582 ms)
(daxpy: 211 ms)
elapsed time: 793 ms
Y+2*X
elapsed time: 1574 ms
Considered the above value, the explanation is that in "Y+2*X"
there are *two* copies of a
10000x10000 matrix instead of only one in "daxpy(2,X,Y)". In
"Y+2*X+3*Z" there will be three
copies, although there could be only one if daxpy was used twice.
I am not blaming Scilab here, I am just blaming "vectorization",
which can be inefficient when
large objects are used. That's why explicits loops can
sometimes be faster than
vectorized operations in Matlab or Julia (which both use JIT compilation).
S.
Le 15/02/2018 à 17:11, Antoine ELIAS a écrit :
> Hello Stéphane,
>
> Interesting ...
>
> In release, we don't ship the header of BLAS/LAPACK functions.
> But you can define them in your C file as extern. ( and let the
linker do his job )
>
> extern int C2F(dgemm) (char *_pstTransA, char *_pstTransB, int
*_piN, int *_piM, int *_piK,
> double *_pdblAlpha, double *_pdblA, int *_piLdA,
> double *_pdblB, int *_piLdB, double
*_pdblBeta, double *_pdblC, int
> *_piLdC);
> and
>
> extern int C2F(dscal) (int *_iSize, double *_pdblVal, double
*_pdblDest, int *_iInc);
>
> Others BLAS/LAPACK prototypes can be found at
http://cgit.scilab.org/scilab/tree/scilab/modu
> les/elementary_functions/includes/elem_common.h?h=6.0
>
> Regards,
> Antoine
> Le 15/02/2018 à 16:50, Stéphane Mottelet a écrit :
> > Hello all,
> >
> > Following the recent discussion with fujimoto, I discovered
that Scilab does not (seem to)
> > fully use SIMD operation in BLAS as it should. Besides the
bottlenecks of its code, there
> > are also many operations of the kind
> >
> > scalar*matrix
> >
> > Althoug this operation is correctly delegated to the DSCAL
BLAS function (can be seen in C
> > function iMultiRealMatrixByRealMatrix in
> > modules/ast/src/c/operations/matrix_multiplication.c) :
> >
> > > int iMultiRealScalarByRealMatrix(
> > > double _dblReal1,
> > > double *_pdblReal2, int _iRows2, int _iCols2,
> > > double *_pdblRealOut)
> > > {
> > > int iOne = 1;
> > > int iSize2 = _iRows2 * _iCols2;
> > >
> > > C2F(dcopy)(&iSize2, _pdblReal2, &iOne, _pdblRealOut, &iOne);
> > > C2F(dscal)(&iSize2, &_dblReal1, _pdblRealOut, &iOne);
> > > return 0;
> > > }
> > in the code below the product "A*1" is likely using only one
processor core, as seen on
> > the cpu usage graph and on the elapsed time,
> >
> > A=rand(20000,20000);
> > tic; for i=1:10; A*1; end; toc
> >
> > ans =
> >
> > 25.596843
> >
> > but this second piece of code is more than 8 times faster and
uses 100% of the cpu,
> >
> > ONE=ones(20000,1);
> > tic; for i=1:10; A*ONE; end; toc
> >
> > ans =
> >
> > 2.938314
> >
> > with roughly the same number of multiplications. This second
computation is delegated to
> > DGEMM (C<-alpha*A*B + beta*C, here with alpha=1 and beta=0)
> >
> > > int iMultiRealMatrixByRealMatrix(
> > > double *_pdblReal1, int _iRows1, int _iCols1,
> > > double *_pdblReal2, int _iRows2, int _iCols2,
> > > double *_pdblRealOut)
> > > {
> > > double dblOne = 1;
> > > double dblZero = 0;
> > >
> > > C2F(dgemm)("n", "n", &_iRows1, &_iCols2, &_iCols1, &dblOne,
> > > _pdblReal1 , &_iRows1 ,
> > > _pdblReal2, &_iRows2, &dblZero,
> > > _pdblRealOut , &_iRows1);
> > > return 0;
> > > }
> > Maybe my intuition is wrong, but I have the feeling that
using dgemm with alpha=0 will be
> > faster than dscal. I plan to test this by making a quick and
dirty code linked to Scilab
> > so my question to devs is : which are the #includes to add on
top of the source (C) to be
> > able to call dgemm and dscal ?
> >
> > Thanks for your help
> >
> > S.
> >
> >
>
> _______________________________________________
> dev mailing list
> [email protected]
> http://lists.scilab.org/mailman/listinfo/dev
_______________________________________________
dev mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/dev
