Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #59451 > unrolled thread

matrix multiplication

Back to article view | Back to comp.lang.forth

Contents

  matrix multiplication krishna.myneni@ccreweb.org - 2017-05-17 15:13 -0700
    Re: matrix multiplication stephen@mpeforth.com (Stephen Pelc) - 2017-05-17 22:31 +0000
      Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-17 15:39 -0700
        Re: matrix multiplication rickman <gnuarm@gmail.com> - 2017-05-18 03:22 -0400
          Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-18 00:29 -0700
            Re: matrix multiplication rickman <gnuarm@gmail.com> - 2017-05-18 03:48 -0400
            Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-18 08:42 +0000
    Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-18 00:25 -0700
      Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-19 15:13 +0000
        Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-19 17:26 +0000
          Re: matrix multiplication Alex <alex@rivadpm.com> - 2017-05-19 18:56 +0100
          Re: matrix multiplication Paul Rubin <no.email@nospam.invalid> - 2017-05-19 21:48 -0700
          Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-20 10:21 +0000
            Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-20 06:12 -0700
            Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-20 09:47 -0700
              Re: matrix multiplication Paul Rubin <no.email@nospam.invalid> - 2017-05-20 15:31 -0700
                Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-21 04:46 -0700
                  Re: matrix multiplication mhx@iae.nl - 2017-05-21 13:41 -0700
                    Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-21 15:32 -0700
              Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-21 07:19 +0000
                Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-21 05:10 -0700
                  Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-21 14:29 +0000
                    Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-21 09:17 -0700
        Re: matrix multiplication stephen@mpeforth.com (Stephen Pelc) - 2017-05-19 20:54 +0000
        Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-19 14:33 -0700
          Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-19 16:12 -0700
          Re: matrix multiplication anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-05-20 14:48 +0000
      Re: matrix multiplication krishna.myneni@ccreweb.org - 2017-05-19 18:38 -0700

Page 2 of 2 — ← Prev page 1 [2]

#59644

On Sunday, May 21, 2017 at 2:51:44 AM UTC-5, Anton Ertl wrote:
> krishna...org writes:
> >A next step could be to do the computations in parallel on a multi-core system, following the example I posted earlier:
> >
> >ftp://ccreweb.org/software/kforth/examples/parallel-mm.4th
> >
> >Do gforth, vfxlin, and swift forth provide fork and mmap syscalls, or equivalents?
> 
> There is multi-threading support in these systems (for Gforth in the
> development version) that you could use; then you also don't need to
> share the memory with mmap.
> 

But I think the idea of forking the process, rather than using threads has some potentially significant benefits to the programmer. The child process(es) inherit a copy of all the code and data memory of the parent, and this eliminates a lot of worries about code re-entrancy and mutexes. The big data would be held in shared memory, and, for some cases, mutexes for the shared memory may be required, but not for matrix multiplication. Each process writes into a separate area of the shared memory even when the output matrix is contiguous.

One thing I'm unsure about are the conditions under which the OS decides to assign a child process to another core. The same applies to threads. From my experiments with fork under Linux thus far, it appears that when cores are idle, the forked process is assigned to a different core. However, it would be useful to determine in advance that this will be the case.
 

> However, my next step would be to implement cache blocking.  While the
> slowdown from going beyond the cache is (to me) surprisingly low for
> the FAXPY-based version, that is on a system that is otherwise idle
> (interestingly both on my PC and on the Odroid C2).  If the system
> makes more use of its multiple cores, and especially if you want to
> have several instances of matrix multiply running on them,using less
> main memory bandwidth by using cache blocking should help.
> 

Dealing with the processor cache is one area in which I'm not well versed. Are there instructions for blocking and unblocking the cache? I assumed the cache was largely independent of programmer control.

> Another thing that I would (in general) do before using multiple cores
> is vectorization.  Unfortunately, we don't have vector support in
> Forth, not even in a system-specific way, yet (AFAIK).
> 

Isn't explicit use of the SSE2 instructions considered vectorization? 

> And if I am really interested in multiplying a lot of matrices
> quickly, I would use a library specialized for that, i.e., ATLAS, or
> one provided by the CPU manufacturer.
> 

There are applications problems which range across a vast scale, and Forth implementations can be convenient for some scales. As Marcel pointed out in another thread, matrix multiplication is only one of different types of computations which, together, are useful in an application. Another one is finding eigenvalues of large matrices. Indeed, use of library routines specially written for such purpose is warranted when the focus is on solving the application problem in a short time. But, it's also useful to have Forth code available for important algorithms, for the purpose of studying their behavior.

> > We should be able to realize a speedup of ~n times by distributing the computing over n cores. Matrix multiplication is an easy computation to parallelize.
> 
> But you realize this speedup only if your computation does not become
> memory-bandwidth bound.
> 

Yes, I realize it after thinking about it. That was a naive statement on my part, based on experience with n=2. I don't have a good idea of how well the scaling degrades at higher n because of cache and memory-bandwidth issues. It would be interesting to experiment with an example such as parallel-mm.4th on a machine with 8 cores.

Krishna

[toc] | [prev] | [next] | [standalone]

#59653

krishna.myneni@ccreweb.org writes:
>On Sunday, May 21, 2017 at 2:51:44 AM UTC-5, Anton Ertl wrote:
>> krishna...org writes:
>> >A next step could be to do the computations in parallel on a multi-core =
>system, following the example I posted earlier:
>> >
>> >ftp://ccreweb.org/software/kforth/examples/parallel-mm.4th
>> >
>> >Do gforth, vfxlin, and swift forth provide fork and mmap syscalls, or eq=
>uivalents?
>>=20
>> There is multi-threading support in these systems (for Gforth in the
>> development version) that you could use; then you also don't need to
>> share the memory with mmap.
>>=20
>
>But I think the idea of forking the process, rather than using threads has =
>some potentially significant benefits to the programmer. The child process(=
>es) inherit a copy of all the code and data memory of the parent, and this =
>eliminates a lot of worries about code re-entrancy and mutexes. The big dat=
>a would be held in shared memory, and, for some cases, mutexes for the shar=
>ed memory may be required, but not for matrix multiplication. Each process =
>writes into a separate area of the shared memory even when the output matri=
>x is contiguous.

Use user variables for the stuff you don't want to be shared, and you
are fine.

>> However, my next step would be to implement cache blocking.  While the
>> slowdown from going beyond the cache is (to me) surprisingly low for
>> the FAXPY-based version, that is on a system that is otherwise idle
>> (interestingly both on my PC and on the Odroid C2).  If the system
>> makes more use of its multiple cores, and especially if you want to
>> have several instances of matrix multiply running on them,using less
>> main memory bandwidth by using cache blocking should help.
>>=20
>
>Dealing with the processor cache is one area in which I'm not well versed. =
>Are there instructions for blocking and unblocking the cache? I assumed the=
> cache was largely independent of programmer control.

Cache blocking means that you process the data in blocks that fit in
the cache.  There are no special instructions for it.  The cache is
handled mostly transparently (apart from performance) by the CPU, but
programmers who know how they work can make use of that knowledge to
get better performance.

In the present case, if we look at your inner-product-based matrix
multiplication, you compute

C[i,j] = A[i,.] * B[.,j]

B[.,j] makes good use of spatial locality (typically, 8 elements of
B[.,j] fit in one cache line), A[i,.] doesn't.

But if we assume that all elements of A[i,.] fit in the cache level of
interest here, you can reuse the elements of A[i,.] by first computing
several elements of C[i,.] before changing the i.

However, if you compute all elements of C[i,.] before changing the i,
the whole B matrix will be read, so you have no temporal locality on
B.  So let's compute several elements of C[i,0..15] first, then
C[i+1,0..15] etc. until C[i+15,0..15]; in this way you need only
A[i..i+15,.] and B[.,0..15] in the cache, for computing 256 elements
of C.  Then compute the next 256 elements of C, e.g.,
C[i..i+15,16..31]; this will reuse the data from A, but load the next
block of B data into the cache.

The block sizes can be adapted to the cache sizes, or you can use
recursion to have blocking at a lot of sizes (some of which will be
good fits for the cache sizes).

>Isn't explicit use of the SSE2 instructions considered vectorization?=20

If you are using SIMD/vector instructions (e.g., MULPD), yes.  If you
are using scalar instructions (e.g., MULSD), no.

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: http://www.forth200x.org/forth200x.html
   EuroForth 2017: http://www.euroforth.org/ef17/

[toc] | [prev] | [next] | [standalone]

#59668

On Sunday, May 21, 2017 at 9:54:02 AM UTC-5, Anton Ertl wrote:
...
> Cache blocking means that you process the data in blocks that fit in
> the cache.  There are no special instructions for it.  The cache is
> handled mostly transparently (apart from performance) by the CPU, but
> programmers who know how they work can make use of that knowledge to
> get better performance.
> 
> In the present case, if we look at your inner-product-based matrix
> multiplication, you compute
> 
> C[i,j] = A[i,.] * B[.,j]
> 
> B[.,j] makes good use of spatial locality (typically, 8 elements of
> B[.,j] fit in one cache line), A[i,.] doesn't.
> 
> But if we assume that all elements of A[i,.] fit in the cache level of
> interest here, you can reuse the elements of A[i,.] by first computing
> several elements of C[i,.] before changing the i.
> 
> However, if you compute all elements of C[i,.] before changing the i,
> the whole B matrix will be read, so you have no temporal locality on
> B.  So let's compute several elements of C[i,0..15] first, then
> C[i+1,0..15] etc. until C[i+15,0..15]; in this way you need only
> A[i..i+15,.] and B[.,0..15] in the cache, for computing 256 elements
> of C.  Then compute the next 256 elements of C, e.g.,
> C[i..i+15,16..31]; this will reuse the data from A, but load the next
> block of B data into the cache.
> 
> The block sizes can be adapted to the cache sizes, or you can use
> recursion to have blocking at a lot of sizes (some of which will be
> good fits for the cache sizes).
> ...

Thanks for the detailed explanation. This is helpful to me.

Krishna

[toc] | [prev] | [next] | [standalone]

#59566

>Forth-94                            Forth-2012
>nc1 @ 0 DO                          nc1 @ 0 DO        
>  2dup 2>r >r f@ r> f@ f* f+          2dup f@ f@ f* f+
>  2r> roffs2 @ + >r                   roffs2 @ + >r   
>  dfloat+ r>                          dfloat+ r>      
>LOOP                                LOOP              

>Here's the code that VFX generates for the inner loop
>
>( 080C6E60    DB2B )                  FLD       TBYTE 0 [EBX]
>( 080C6E62    8B5500 )                MOV       EDX, [EBP]
>( 080C6E65    DB2A )                  FLD       TBYTE 0 [EDX]
>( 080C6E67    DEC9 )                  FMULP     ST(1), ST
>( 080C6E69    DEC1 )                  FADDP     ST(1), ST
>( 080C6E6B    031D70340A08 )          ADD       EBX, [080A3470]
>( 080C6E71    53 )                    PUSH      EBX
>( 080C6E72    8B5D00 )                MOV       EBX, [EBP]
>( 080C6E75    83C308 )                ADD       EBX, 08
>( 080C6E78    5A )                    POP       EDX
>( 080C6E79    83042401 )              ADD       [ESP], 01
>( 080C6E7D    8344240401 )            ADD       [ESP+04], 01
>( 080C6E82    895D00 )                MOV       [EBP], EBX
>( 080C6E85    8BDA )                  MOV       EBX, EDX
>( 080C6E87    71D7 )                  JNO       080C6E60

For VFX, try replacing
  >r dfloat r>
with
  swap dfloat+ swap

Stephen


-- 
Stephen Pelc, stephenXXX@mpeforth.com
MicroProcessor Engineering Ltd - More Real, Less Time
133 Hill Lane, Southampton SO15 5AF, England
tel: +44 (0)23 8063 1441, fax: +44 (0)23 8033 9691
web: http://www.mpeforth.com - free VFX Forth downloads

[toc] | [prev] | [next] | [standalone]

#59568

On Friday, May 19, 2017 at 12:07:34 PM UTC-5, Anton Ertl wrote:
> krishna...org writes:
> >On Wednesday, May 17, 2017 at 5:13:50 PM UTC-5, krishna...@ccreweb.org wrot=
> >e:
> >> ... I spent an hour writing a full matrix multiplier for an integrated st=
> >ack system, source below. ...
> >
> >The row offsets may be precomputed for the two matrices to be multiplied. T=
> >he revised Forth source is given below.
> 
> I took this code, replaced occurences of A@ with @ and replaced
> 
> >        J I df_mul_r1c2 2 pick f!
> >        dfloat+
> 
> with
> 
>         J I rot >r df_mul_r1c2 r@ f!
>         r> dfloat+
> 
> and the result seems to be working; I did not find this replacement
> trivial, BTW. ...

A temporary variable (or a local) for storing the address would make this simpler, cleaner in appearance and portable to all standard Forth systems which have fp support. I don't know if it would be slower to do this, though.

> ...  I also produced a version of the inner loop that works
> only on separate-stack systems (i.e., Forth-2012):
> 
> Forth-94                            Forth-2012
> nc1 @ 0 DO                          nc1 @ 0 DO        
>   2dup 2>r >r f@ r> f@ f* f+          2dup f@ f@ f* f+
>   2r> roffs2 @ + >r                   roffs2 @ + >r   
>   dfloat+ r>                          dfloat+ r>      
> LOOP                                LOOP              
> 

This works under gforth. See the following link for additional test code:

ftp://ccreweb.org/software/gforth/fsl/extras/mmul.fs


> Results for a 500x500 matrix multiply on a Core i5-6600K:
> 
> Forth-94 Forth-2012 Forth-2012 Forth-94 Forth-94 
> Myneni  Myneni/Ertl     Ertl   Rees/Ertl   pahihu 
> 109966M     107831M    1046M       1202M    1517M vfxlin Ndp387.fth
>   2986M       2260M    2220M       2945M    2764M gforth-fast
>   4028M       4076M    4844M       5067M    5017M SwiftForth HW stack, no WAIT
> ...

Can you make available the programs you ran for these tests? They may have been discussed in the other thread, but it would be helpful to collect them in a single public folder, for comparison.
 
Krishna

[toc] | [prev] | [next] | [standalone]

#59569

On Friday, May 19, 2017 at 4:33:06 PM UTC-5, krishna...@ccreweb.org wrote:
...
> Can you make available the programs you ran for these tests? They may have been discussed in the other thread, but it would be helpful to collect them in a single public folder, for comparison.
>  

I found them. I overlooked them in my first glance at the matmul package.

Krishna

[toc] | [prev] | [next] | [standalone]

#59599

krishna.myneni@ccreweb.org writes:
>Can you make available the programs you ran for these tests?

The version of your program that I used is

<http://www.complang.tuwien.ac.at/forth/programs/mmul.4th>

Note the 

: matmulr df_mmul ;

at the end:-).

For benchmarking I used
<http://theforth.net/package/matmul/current-view/benchmark.4th>

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: http://www.forth200x.org/forth200x.html
   EuroForth 2017: http://www.euroforth.org/ef17/

[toc] | [prev] | [next] | [standalone]

#59570

On Thursday, May 18, 2017 at 2:25:45 AM UTC-5, krishna...org wrote:

I have to fix a problem with the examples for the df_mmul tests. I inadvertently allotted double integer matrices rather than double precision floats. On a 32-bit system, it works because they are the same size. DOUBLE is defined as a constant equal to 2 CELLS in the fsl-util.x file. The data size should be equal to 1 DFLOATS. Strangely enough, we don't have a constant named DFLOAT to represent that size in the fsl-util.x file. To fix the examples, we should first define

1 DFLOATS constant DFLOAT

> 
> 100 100 double matrix a{{
>  ok
> 100 100 double matrix b{{
>  ok
> 100 100 double matrix c{{
>  ...

should be

  100 100 dfloat matrix a{{
  100 100 dfloat matrix b{{
  100 100 dfloat matrix c{{

> 
> 2 2 double matrix a{{
> 2 2 double matrix b{{
> 2 2 double matrix c{{
> ...

should be

  2 2 dfloat matrix a{{
  2 2 dfloat matrix b{{
  2 2 dfloat matrix c{{

I'm going to revise the fsl-util.x file to include the constant DFLOAT.

Krishna

[toc] | [prev] | [standalone]

Page 2 of 2 — ← Prev page 1 [2]

Back to top | Article view | comp.lang.forth

csiph-web