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

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

pde2 floating point benchmark code

Started bykrishna.myneni@ccreweb.org
First post2013-03-04 18:02 -0800
Last post2013-03-08 16:43 -0800
Articles 9 on this page of 29 — 8 participants

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

Contents

  pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-04 18:02 -0800
    Re: pde2 floating point benchmark code Mark Wills <markrobertwills@yahoo.co.uk> - 2013-03-04 23:23 -0800
      Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 05:22 -0800
    Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-05 12:11 -0500
      Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 16:11 -0800
        Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-06 04:26 -0500
          Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-06 05:22 -0800
            Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-06 09:44 -0500
    Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 16:17 -0800
      Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-06 23:06 +0200
        Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-07 04:57 -0800
          Re: pde2 floating point benchmark code m.a.m.hendrix@tue.nl - 2013-03-07 05:49 -0800
          Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-07 21:12 +0200
            Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-07 16:07 -0800
        Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-10 09:17 -0400
          Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-10 18:33 +0200
            Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-11 04:54 -0400
      Re: pde2 floating point benchmark code humptydumpty <ouatubi@gmail.com> - 2013-03-07 22:32 -0800
        Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 05:19 -0800
          Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-08 21:04 +0200
            Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 16:31 -0800
              Re: pde2 floating point benchmark code Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2013-03-09 00:56 +0000
                Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 19:52 -0800
                  Re: pde2 floating point benchmark code Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2013-03-09 09:09 +0000
              Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-09 13:37 +0200
                Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-09 06:28 -0800
                Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-09 17:58 +0200
          Re: pde2 floating point benchmark code humptydumpty <ouatubi@gmail.com> - 2013-03-08 12:23 -0800
            Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 16:43 -0800

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

#20474

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2013-03-08 16:31 -0800
Message-ID<5a243597-c8bf-4e7c-8611-a5a88245429c@l16g2000yqe.googlegroups.com>
In reply to#20463
On Mar 8, 1:04 pm, m...@iae.nl (Marcel Hendrix) wrote:
> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
>

Ok. Here are some timing results for execution time of 25,000 time
steps of the diffusion calculation, on my system:

diffusion.R (optimized R byte code):
-----------
user time: 60.189 s

diffusion-fast.R (inline compiled Fortran code):
-----------------
user time: 0.308 s

kforth:
-------
elapsed: 9.412 s

kforth-fast:
-----------
elapsed: 6.779 s

I'm having difficulty with building gforth from source, so unable to
provide those numbers at present. My system has a Athlon 64 X2 Dual
Core Processor 5200+ at 1.000 GHz.


> > On Mar 8, 12:32am, humptydumpty <ouat...@gmail.com> wrote:
>
> [..]
>
> > That's certainly an interesting approach for this problem, to unroll
> > the inner loop. A factor of two speedup is significant -- it may not
> > hold up for a native code compiler, where the loop overhead is likely
> > to be substantially lower.
>
> Not because of the loop overhead, but because the indexed
> address is calculated at compile time.
>

Thanks! I did not notice that earlier.

> ...
> > For raw speed improvement, I would consider
> > writing the inner loop as an assembler code word. That should provide
> > a huge improvement on indirect threaded code systems. Also, I found
> > that replacing the array indexing with simpler pointer arithmetic in
> > the inner loop also led to close to a factor of two improvement.
>
> Or use a different Forth, when speed is really a concern :-)
>

No reason why native code Forth systems should have all the fun :) .
The inner loop of evolve should be easy to code in assembler.

> > However, speed is only one concern with this method of calculation.
> > The grid size, set by dt and dx, must be small to obtain reasonable
> > accuracy of the solution, even for short time intervals such as t = 1.
>
> Yes, that is really a problem. What is the error of this method?
>

Don't know yet. I hope to investigate this presently. It's a good
method to illustrate the use of computing for solving partial
differential equations, since the method is very easy to understand
(just replace the partial derivatives with finite difference
approximations).

> FORTH> in
> nx: 1001
> fiber size:44008
>
> 0.135 seconds elapsed.
> 0.059 seconds elapsed. ok

Speedup of a factor of ~2.

>
> | Output:
> | ~/wrk <> gforth-fast pde.fs
> | redefined DT with dt
> | nx: 1001
> | fiber size:87940
>
> | 1451235
> | 1116868
>
> Microseconds, so 1.451 and 1.117 seconds.
> Amazingly, the Gforth code is almost twice *larger* than iForth native code.
>
> -marcel

Nice. By the way, there's an extraneous line of code in the original
Forth program, which may be dispensed with. It was left over from
translation of the R program, where it was needed to initialize a
local array.

: evolve ( n -- )
    u_j u_jp1 nx @ FLOATS move   <--- this line is not needed.
    0 ?DO
      nx @ 1- 1 DO
       :


Cheers,
Krishna

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

#20476

FromZbiggy <zbigniew2011REMOVE@gmail.REMOVE.com>
Date2013-03-09 00:56 +0000
Message-ID<slrnkjl5cg.bf4.zbigniew2011REMOVE@Tichy.myhome.org>
In reply to#20474
In comp.lang.forth, Krishna Myneni wrote:

> I'm having difficulty with building gforth from source, so unable to
> provide those numbers at present. My system has a Athlon 64 X2 Dual
> Core Processor 5200+ at 1.000 GHz.

It builds just fine - but use newest 7.0.2
-- 
The consensus was, as usual in this community, that there is no consensus. (RA)

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

#20478

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2013-03-08 19:52 -0800
Message-ID<ff2f61ea-d046-459e-834f-5ecbbf921115@h17g2000yqe.googlegroups.com>
In reply to#20476
On Mar 8, 6:56 pm, Zbiggy <zbigniew2011REM...@gmail.REMOVE.com> wrote:
> In comp.lang.forth, Krishna Myneni wrote:
>
> > I'm having difficulty with building gforth from source, so unable to
> > provide those numbers at present. My system has a Athlon 64 X2 Dual
> > Core Processor 5200+ at 1.000 GHz.
>
> It builds just fine - but use newest 7.0.2
> --
> The consensus was, as usual in this community, that there is no consensus. (RA)

On my system, the build process fails at some point. However, with
7.0.2, it fails after making the gforth and gforth-fast executables.
The execution times for "25000 evolve", using my originally posted
source on my system are,

gforth:      5.332 s
gforth-fast: 1.910 s


With H.D.'s loop unrolling method, those times could be cut by more
than half.

Krishna

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

#20487

FromZbiggy <zbigniew2011REMOVE@gmail.REMOVE.com>
Date2013-03-09 09:09 +0000
Message-ID<slrnkjm295.1st.zbigniew2011REMOVE@Tichy.myhome.org>
In reply to#20478
In comp.lang.forth, Krishna Myneni wrote:

> On my system, the build process fails at some point. However, with
> 7.0.2, it fails after making the gforth and gforth-fast executables.

You must be facing "install-info" problem. It's enough to change the contents
of Makefile.in a little, before compilation. There is a line:

INSTALL_INFO = @INSTALL_INFO@

Just delete the right side (@INSTALL_INFO@), and it won't try to execute
install-info anymore.
-- 
The consensus was, as usual in this community, that there is no consensus. (RA)

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

#20489

Frommhx@iae.nl (Marcel Hendrix)
Date2013-03-09 13:37 +0200
Message-ID<19669094008434@frunobulax.edu>
In reply to#20474
Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: pde2 floating point benchmark code

> On Mar 8, 1:04 pm, m...@iae.nl (Marcel Hendrix) wrote:
>> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
[..]
>>> However, speed is only one concern with this method of calculation.
>>> The grid size, set by dt and dx, must be small to obtain reasonable
>>> accuracy of the solution, even for short time intervals such as t = 1.
>>
>> Yes, that is really a problem. What is the error of this method?
>>

> Don't know yet. I hope to investigate this presently. It's a good
> method to illustrate the use of computing for solving partial
> differential equations, since the method is very easy to understand
> (just replace the partial derivatives with finite difference
> approximations).

Here's an experimental approach.

The outcome on my system is that the *difference between iterations* 
goes down by a factor of four when dx is halved and dt adjusted so that 
fx = 0.5. Each of these iteration takes eight times longer than the 
previous one.

This does not say anything about the absolute / relative accuracy, but
allows one to stop when the solver makes too little progress. E.g. 
when u_ext is 20 degrees C, a good stop would be when the difference
between iterations is less than 0.01 degree.

FORTH> 1e-5 GO
dx = 1.000000e-1 dt = 1.000000e-2 error 4.153411e0  0.001 seconds elapsed.
dx = 5.000000e-2 dt = 2.500000e-3 error 6.851562e-2  0.002 seconds elapsed.
dx = 2.500000e-2 dt = 6.250000e-4 error 1.693734e-2  0.004 seconds elapsed.
dx = 1.250000e-2 dt = 1.562500e-4 error 4.222571e-3  0.023 seconds elapsed.
dx = 6.250000e-3 dt = 3.906250e-5 error 1.054911e-3  0.168 seconds elapsed.
dx = 3.125000e-3 dt = 9.765625e-6 error 2.636820e-4  1.349 seconds elapsed.
dx = 1.562500e-3 dt = 2.441406e-6 error 6.591764e-5  10.810 seconds elapsed.
dx = 7.812500e-4 dt = 6.103516e-7 error 1.647923e-5  86.704 seconds elapsed.
dx = 3.906250e-4 dt = 1.525879e-7 error 4.119797e-6  696.044 seconds elapsed. ok

-marcel

-- ---------

ANEW -diffuse2

0.01e		 FVALUE dx 	\ 0.01 m
0.5e 		 FVALUE D	\ Diffusion coefficient
20e 		 FVALUE u_ext    
0.5e 		 FVALUE fk	\ choose dt so that fk < 0.5
dx FSQR f2/ D F/ FVALUE dt
1e fk F2* F-     FVALUE 1-2fk 

#10 		  VALUE nx

1e-5 FCONSTANT dxmin
10e dxmin F/ F>S 1+ =: nxmax  

CREATE x      nxmax DFLOATS ALLOT
CREATE u_x_0  nxmax DFLOATS ALLOT
CREATE u_j    nxmax DFLOATS ALLOT
CREATE u_jp1  nxmax DFLOATS ALLOT
CREATE golden 0e DF, 0e DF, 0e DF,

\ Compute the solution at t=1
\ Compute the absolute difference with the previous iteration
: res? ( F: -- err ) 
	CR ." dx = "  dx E. ." dt = "  dt E. ." error "  
	golden 0 DFLOAT[] DF@   u_j nx 7 #10 */ DFLOAT[] DF@ FDUP golden 0 DFLOAT[] DF!  F- FSQR
	golden 1 DFLOAT[] DF@   u_j nx 8 #10 */ DFLOAT[] DF@ FDUP golden 1 DFLOAT[] DF!  F- FSQR F+
	golden 2 DFLOAT[] DF@   u_j nx 9 #10 */ DFLOAT[] DF@ FDUP golden 2 DFLOAT[] DF!  F- FSQR F+ 
	3e F/ FSQRT FDUP E. SPACE ;

: evolve ( n -- ) ( F: -- res )
	2/ 
	0 ?DO
		nx 1- 1 DO
			  u_j   I 1+ DFLOAT[] DF@  
			  u_j   I 1- DFLOAT[] DF@ F+ fk F*
			  u_j   I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_jp1 I    DFLOAT[] DF!
		      LOOP
		u_jp1 1 DFLOAT[] DF@ u_jp1 0     DFLOAT[] DF!   \ b.c.: u_x(0, t) = 0
		u_ext                u_jp1 nx 1- DFLOAT[] DF!   \ b.c.: u(10, t) = u_ext

		nx 1- 1 DO
			  u_jp1 I 1+ DFLOAT[] DF@  
			  u_jp1 I 1- DFLOAT[] DF@ F+ fk F*
			  u_jp1 I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_j   I    DFLOAT[] DF!
		      LOOP
		u_j   1 DFLOAT[] DF@ u_j   0     DFLOAT[] DF!   \ b.c.: u_x(0, t) = 0
		u_ext                u_j   nx 1- DFLOAT[] DF!   \ b.c.: u(10, t) = u_ext
	  LOOP  res? ;

: INIT ( F: dx -- )
	TO dx
	dx FSQR F2/ D F/ TO dt
	10e dx F/ F>S 1+ TO nx
    	nx nxmax > ABORT" init :: dx too small for preallocated array size" 
	0e nx 0 DO FDUP x     I DFLOAT[] DF! dx F+  LOOP  FDROP
	   nx 0 DO  1e  u_x_0 I DFLOAT[] DF!  	    LOOP
	u_x_0 u_j nx DFLOATS MOVE ;

-- Stop when difference between two last iterations smaller than ERROR
: GO ( F: error -- )
	     FLOCAL error
	0.1e FLOCAL xstart 
	golden 3 DFLOATS ERASE
	BEGIN  
	   xstart INIT  
	   TIMER-RESET 1e ( second) dt F/ F>S EVOLVE .ELAPSED
	   error F>  EKEY? 0= AND
	WHILE
	   xstart F2/ TO xstart
	REPEAT ;

0.01e init CR TIMER-RESET 1e dt F/ F>S evolve FDROP .ELAPSED

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

#20490

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2013-03-09 06:28 -0800
Message-ID<fa7388dd-f8e0-488f-9ac5-845e1907e513@u7g2000yqg.googlegroups.com>
In reply to#20489
On Mar 9, 5:37 am, m...@iae.nl (Marcel Hendrix) wrote:
> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
>
>
>
> > On Mar 8, 1:04 pm, m...@iae.nl (Marcel Hendrix) wrote:
> >> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
> [..]
> >>> However, speed is only one concern with this method of calculation.
> >>> The grid size, set by dt and dx, must be small to obtain reasonable
> >>> accuracy of the solution, even for short time intervals such as t = 1.
>
> >> Yes, that is really a problem. What is the error of this method?
>
> > Don't know yet. I hope to investigate this presently. It's a good
> > method to illustrate the use of computing for solving partial
> > differential equations, since the method is very easy to understand
> > (just replace the partial derivatives with finite difference
> > approximations).
>
> Here's an experimental approach.
>
> The outcome on my system is that the *difference between iterations*
> goes down by a factor of four when dx is halved and dt adjusted so that
> fx = 0.5. Each of these iteration takes eight times longer than the
> previous one.
>
> This does not say anything about the absolute / relative accuracy, but
> allows one to stop when the solver makes too little progress. E.g.
> when u_ext is 20 degrees C, a good stop would be when the difference
> between iterations is less than 0.01 degree.
>
> FORTH> 1e-5 GO
> dx = 1.000000e-1 dt = 1.000000e-2 error 4.153411e0  0.001 seconds elapsed.
> dx = 5.000000e-2 dt = 2.500000e-3 error 6.851562e-2  0.002 seconds elapsed.
> dx = 2.500000e-2 dt = 6.250000e-4 error 1.693734e-2  0.004 seconds elapsed.
> dx = 1.250000e-2 dt = 1.562500e-4 error 4.222571e-3  0.023 seconds elapsed.
> dx = 6.250000e-3 dt = 3.906250e-5 error 1.054911e-3  0.168 seconds elapsed.
> dx = 3.125000e-3 dt = 9.765625e-6 error 2.636820e-4  1.349 seconds elapsed.
> dx = 1.562500e-3 dt = 2.441406e-6 error 6.591764e-5  10.810 seconds elapsed.
> dx = 7.812500e-4 dt = 6.103516e-7 error 1.647923e-5  86.704 seconds elapsed.
> dx = 3.906250e-4 dt = 1.525879e-7 error 4.119797e-6  696.044 seconds elapsed. ok
>
> -marcel
>

In the output above, "error" is a misnomer, as you indicated. It's
really an indicator of convergence. It looks like the finite
difference method will be really sluggish if one wants to achieve more
than 3 significant digits in the solution. Now, you've really got me
curious about how many significant digits are in the ReacTran
calculation, which is very fast, even out to t=100.

How did you choose the particular positions x={7,8,9} to measure the
convergence?

Krishna

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

#20494

Frommhx@iae.nl (Marcel Hendrix)
Date2013-03-09 17:58 +0200
Message-ID<93771894008434@frunobulax.edu>
In reply to#20489
mhx@iae.nl (Marcel Hendrix) wrote Re: pde2 floating point benchmark code

> Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: pde2 floating point benchmark code

>> On Mar 8, 1:04 pm, m...@iae.nl (Marcel Hendrix) wrote:
>>> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
[..]
> Here's an experimental approach.

[..]FORTH> 1e-5 GO
> dx = 1.000000e-1 dt = 1.000000e-2 error 4.153411e0  0.001 seconds elapsed.
> dx = 5.000000e-2 dt = 2.500000e-3 error 6.851562e-2  0.002 seconds elapsed.
> dx = 2.500000e-2 dt = 6.250000e-4 error 1.693734e-2  0.004 seconds elapsed.
> dx = 1.250000e-2 dt = 1.562500e-4 error 4.222571e-3  0.023 seconds elapsed.
> dx = 6.250000e-3 dt = 3.906250e-5 error 1.054911e-3  0.168 seconds elapsed.
> dx = 3.125000e-3 dt = 9.765625e-6 error 2.636820e-4  1.349 seconds elapsed.
> dx = 1.562500e-3 dt = 2.441406e-6 error 6.591764e-5  10.810 seconds elapsed.
> dx = 7.812500e-4 dt = 6.103516e-7 error 1.647923e-5  86.704 seconds elapsed.
> dx = 3.906250e-4 dt = 1.525879e-7 error 4.119797e-6  696.044 seconds elapsed. ok
[..]

This one works for 2 threads (but not for 4) and is indeed two times faster still.
I tried 4 threads, but apparently the interactions become too complex then and the
error doesn't decrease as fast as expected.

FORTH> 1e-5 GO
dx = 1.000000e-1 dt = 1.000000e-2 error 1.050422 1.874837 7.074819 4.268937e0  0.003 seconds elapsed.
dx = 5.000000e-2 dt = 2.500000e-3 error 1.051084 1.867074 7.040390 2.038034e-2  0.003 seconds elapsed.
dx = 2.500000e-2 dt = 6.250000e-4 error 1.051243 1.865146 7.031773 5.098856e-3  0.006 seconds elapsed.
dx = 1.250000e-2 dt = 1.562500e-4 error 1.051283 1.864665 7.029618 1.274944e-3  0.022 seconds elapsed.
dx = 6.250000e-3 dt = 3.906250e-5 error 1.051293 1.864545 7.029079 3.187503e-4  0.108 seconds elapsed.
dx = 3.125000e-3 dt = 9.765625e-6 error 1.051295 1.864515 7.028945 7.968847e-5  0.746 seconds elapsed.
dx = 1.562500e-3 dt = 2.441406e-6 error 1.051296 1.864508 7.028911 1.992217e-5  5.732 seconds elapsed.
dx = 7.812500e-4 dt = 6.103516e-7 error 1.051296 1.864506 7.028902 4.980547e-6  45.438 seconds elapsed.
dx = 3.906250e-4 dt = 1.525879e-7 error 1.051296 1.864505 7.028900 1.245137e-6  460.515 seconds elapsed. ok

-marcel

-- ------------
NEEDS -threads
ANEW -diffuse3

0.01e		 FVALUE dx 	\ 0.01 m
0.5e 		 FVALUE D	\ Diffusion coefficient
20e 		 FVALUE u_ext    
0.5e 		 FVALUE fk	\ choose dt so that fk < 0.5
dx FSQR f2/ D F/ FVALUE dt
1e fk F2* F-     FVALUE 1-2fk 

#10 		  VALUE nx

1e-5 FCONSTANT dxmin
10e dxmin F/ F>S 1+ =: nxmax  

CREATE x      nxmax DFLOATS ALLOT
CREATE u_x_0  nxmax DFLOATS ALLOT
CREATE u_j    nxmax DFLOATS ALLOT
CREATE u_jp1  nxmax DFLOATS ALLOT
CREATE golden 0e DF, 0e DF, 0e DF,

\ Compute the solution at t=1
\ Compute the absolute difference with the previous iteration
: res? ( F: -- err ) 
	CR ." dx = "  dx E. ." dt = "  dt E. ." error "  
	golden 0 DFLOAT[] DF@   u_j nx 7 #10 */ DFLOAT[] DF@ FDUP golden 0 DFLOAT[] DF!  FDUP F.  F- FSQR
	golden 1 DFLOAT[] DF@   u_j nx 8 #10 */ DFLOAT[] DF@ FDUP golden 1 DFLOAT[] DF!  FDUP F.  F- FSQR F+
	golden 2 DFLOAT[] DF@   u_j nx 9 #10 */ DFLOAT[] DF@ FDUP golden 2 DFLOAT[] DF!  FDUP F.  F- FSQR F+ 
	3e F/ FSQRT FDUP E. SPACE ;

: evolve1 ( -- )
	1e ( second) dt F/ F>S	2/  
	0 ?DO
		nx 2/ 1- 1 DO
			  u_j   I 1+ DFLOAT[] DF@  
			  u_j   I 1- DFLOAT[] DF@ F+ fk F*
			  u_j   I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_jp1 I    DFLOAT[] DF!
		      LOOP
		u_jp1 1 DFLOAT[] DF@  u_jp1 0     DFLOAT[] DF!	\ b.c.: u_x(0, t) = 0
		nx 2/ 1- 1 DO
			  u_jp1 I 1+ DFLOAT[] DF@  
			  u_jp1 I 1- DFLOAT[] DF@ F+ fk F*
			  u_jp1 I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_j   I    DFLOAT[] DF!
		      LOOP
		u_j   1 DFLOAT[] DF@  u_j   0     DFLOAT[] DF!  \ b.c.: u_x(0, t) = 0
	  LOOP  ;

: evolve2 ( -- )
	1e ( second) dt F/ F>S	2/  
	0 ?DO
		nx 1-  nx 2/ 1- 
			DO
			  u_j   I 1+ DFLOAT[] DF@  
			  u_j   I 1- DFLOAT[] DF@ F+ fk F*
			  u_j   I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_jp1 I    DFLOAT[] DF!
		      LOOP
		u_jp1 1 DFLOAT[] DF@  u_jp1 0     DFLOAT[] DF!	\ b.c.: u_x(0, t) = 0
		nx 1-  nx 2/ 1-
			DO
			  u_jp1 I 1+ DFLOAT[] DF@  
			  u_jp1 I 1- DFLOAT[] DF@ F+ fk F*
			  u_jp1 I    DFLOAT[] DF@ 1-2fk F* F+  
			  u_j   I    DFLOAT[] DF!
		      LOOP
		u_j   1 DFLOAT[] DF@  u_j   0     DFLOAT[] DF!	\ b.c.: u_x(0, t) = 0
	  LOOP  ;

: evolve ( -- ) ( F: -- error )
	PAR
	  STARTP  evolve1  ENDP
	  STARTP  evolve2  ENDP
	ENDPAR	res? ;

: INIT ( F: dx -- )
	TO dx
	dx FSQR F2/ D F/ TO dt
	10e dx F/ F>S 1+ TO nx
    	nx nxmax > ABORT" init :: dx too small for preallocated array size" 
	0e nx 0 DO FDUP x     I DFLOAT[] DF! dx F+  LOOP  FDROP
	   nx 0 DO  1e  u_x_0 I DFLOAT[] DF!  	    LOOP
	u_ext  u_x_0 nx 1- DFLOAT[] DF!	\ b.c.: u(10, t) = u_ext
	u_x_0  u_j   nx  DFLOATS MOVE 
	u_x_0  u_jp1 nx  DFLOATS MOVE ;


-- Stop when difference between two last iterations not smaller than ERROR
: GO ( F: error -- )
	     FLOCAL error
	0.1e FLOCAL xstart 
	golden 3 DFLOATS ERASE
	BEGIN  
	   xstart INIT  
	   TIMER-RESET  evolve .ELAPSED
	   error F>  EKEY? 0= AND
	WHILE
	   xstart F2/ TO xstart
	REPEAT ;

0.01e init CR TIMER-RESET evolve FDROP .ELAPSED

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

#20464

Fromhumptydumpty <ouatubi@gmail.com>
Date2013-03-08 12:23 -0800
Message-ID<1f2fefc3-3d56-45b8-bd69-ca5278ee802f@googlegroups.com>
In reply to#20441
On Friday, March 8, 2013 3:19:09 PM UTC+2, Krishna Myneni wrote:
> On Mar 8, 12:32 am, humptydumpty <ouat...@gmail.com> wrote:
> 
> > On Wednesday, March 6, 2013 2:17:29 AM UTC+2, krishna...@ccreweb.org wrote:
> 
> > > On Monday, March 4, 2013 8:02:48 PM UTC-6, krishna...@ccreweb.org wrote:
> 
> >
> 
> > > > Here's a useful floating point calculation, solution of the 1-D diffusion equation, which may be useful for benchmarking non-native code Forth systems. A corresponding version in the higher level language, R, may be found at
> 
> >
> 
> > > >ftp://ccreweb.org/software/R/pde/diffusion.R
> 
> >
> 
> > > > KM
> 
> >
> 
> > > ...
> 
> >
> 
> > > The above R program uses a byte-code optimizer available within R. I find that the Forth version, using kforth-fast, is considerably faster than this partially optimized version.
> 
> >
> 
> > > Also, for reference, there is an R version of this program which uses the R package, "inline", to permit the "evolve()" function to be written in Fortran within the R program. The Fortran code is compiled and used transparently, and the speedup is enormous. Native code Forth compilers could use this hybrid R program for comparison. The link is
> 
> >
> 
> > >ftp://ccreweb.org/software/R/pde/diffusion-fast.R
> 
> >
> 
> > > KM
> 
> >
> 
> > Hi Krishna!
> 
> > Is possible to speed-up unrolling inner-loop
> 
> > into one word (a "calculus-fiber"):
> 
> > --
> 
> > \ syntactic sugar for simple fp arrays
> 
> > : []F   ( a n -- a' )
> 
> >         FLOATS + ;
> 
> > : []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
> 
> >     ]] FLOATS + f@ [[ ; immediate
> 
> > : []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
> 
> >     ]] FLOATS + f! [[ ; immediate
> 
> >
> 
> > fvariable D        0.5e D f!      \ Diffusion coefficient
> 
> > fvariable u_ext    20e  u_ext f!
> 
> >
> 
> > 0.01e fconstant dx
> 
> > 4e-5  fconstant dt
> 
> >
> 
> > D f@ dt f* dx fdup f* f/        FVARIABLE fk fk F!      \ choose dt so that fk < 0.5
> 
> > 1e0 2e0 fk f@ f* f-             FVARIABLE 1-2fk 1-2fk F!
> 
> >
> 
> > 10e dx f/ f>d d>s 1+ constant nx
> 
> > cr ." nx: " nx . cr
> 
> >
> 
> > create x      nx FLOATS allot
> 
> > create u_x_0  nx FLOATS allot
> 
> > create u_j    nx FLOATS allot
> 
> > create u_jp1  nx FLOATS allot
> 
> >
> 
> > : calculus-fiber:
> 
> >         :
> 
> >         nx 1- 1 DO
> 
> >                 fk           postpone LITERAL postpone F@
> 
> >                 u_j I 1+ []F postpone LITERAL postpone F@
> 
> >                 u_j I 1- []F postpone LITERAL postpone F@
> 
> >                 postpone F+ postpone F*
> 
> >                 u_j I    []F postpone LITERAL postpone F@
> 
> >                 1-2fk        postpone LITERAL postpone F@
> 
> >                 postpone F* postpone F+
> 
> >                 u_jp1 I  []F postpone LITERAL postpone F!
> 
> >         LOOP
> 
> >         postpone ;
> 
> > ;
> 
> >
> 
> > here calculus-fiber: fiberONE here swap - ." fiber size:" . cr
> 
> > \ see fiberONE
> 
> > : evolve2
> 
> >     u_j u_jp1 nx FLOATS move
> 
> >     0 ?DO
> 
> >         fiberONE
> 
> >         u_jp1 1 []F@ u_jp1 F!             \ b.c.: u_x(0, t) = 0
> 
> >         u_ext f@     u_jp1 nx 1- []F!   \ b.c.: u(10, t) = u_ext
> 
> >
> 
> >         u_jp1 u_j nx FLOATS move
> 
> >     LOOP
> 
> > ;
> 
> >
> 
> > \ Evolve the solution by n time steps, i.e. by an elapsed time
> 
> > \   of n*dt
> 
> > : evolve ( n -- )
> 
> >     u_j u_jp1 nx FLOATS move
> 
> >     0 ?DO
> 
> >       nx 1- 1 DO
> 
> >         u_j I 1+ []F@     u_j I []F@ 2e f* f-  u_j I 1- []F@ f+  fk f@ f*
> 
> >         u_j I    []F@ f+  u_jp1 I []F!
> 
> >       LOOP
> 
> >       u_jp1 1 []F@ u_jp1 F!             \ b.c.: u_x(0, t) = 0
> 
> >       u_ext f@     u_jp1 nx 1- []F!   \ b.c.: u(10, t) = u_ext
> 
> >
> 
> >       u_jp1 u_j nx FLOATS move
> 
> >     LOOP
> 
> > ;
> 
> >
> 
> > : init ( -- )
> 
> >     0e nx 0 DO fdup x I []F! dx f+  LOOP  fdrop
> 
> >     nx 0 DO  1e u_x_0 I []F!  LOOP
> 
> >     u_x_0 u_j nx FLOATS move
> 
> > ;
> 
> > init
> 
> >
> 
> > \ Compute the solution at t=1
> 
> > utime
> 
> > 25000 evolve
> 
> > utime 2swap d- cr d.
> 
> > utime
> 
> > 25000 evolve2
> 
> > utime 2swap d- cr d.
> 
> > cr .s cr bye
> 
> > --
> 
> > Output:
> 
> > ~/wrk <> gforth-fast pde.fs
> 
> > redefined DT with dt
> 
> > nx: 1001
> 
> > fiber size:87940
> 
> >
> 
> > 1451235
> 
> > 1116868
> 
> > <0>
> 
> > ~/wrk <> gforth pde.fs
> 
> > redefined DT with dt
> 
> > nx: 1001
> 
> > fiber size:87940
> 
> >
> 
> > 2036871
> 
> > 1082571
> 
> > <0>
> 
> >
> 
> > That "calculus-fiber" could be used to implement a TDMA.
> 
> > Then domenium of problem could be splitted into "fibers",
> 
> > then iterate domenium through "fibers".
> 
> > ( Convergence is another story. )
> 
> >
> 
> > Have a nice day,
> 
> > humptydumpty
> 
> 
> 
> 
> 
> Hi H.D.,
> 
> 
> 
> That's certainly an interesting approach for this problem, to unroll
> 
> the inner loop. A factor of two speedup is significant -- it may not
> 
> hold up for a native code compiler, where the loop overhead is likely
> 
> to be substantially lower. For raw speed improvement, I would consider
> 
> writing the inner loop as an assembler code word. That should provide
> 
> a huge improvement on indirect threaded code systems. Also, I found
> 
> that replacing the array indexing with simpler pointer arithmetic in
> 
> the inner loop also led to close to a factor of two improvement.
> 
> However, speed is only one concern with this method of calculation.
> 
> The grid size, set by dt and dx, must be small to obtain reasonable
> 
> accuracy of the solution, even for short time intervals such as t = 1.
> 
> I haven't tested the method for calculating the solution out to t=100,
> 
> as is done in the R example using the ReacTran package, which, by the
> 
> way, computes the solution across the entire grid very quickly (~5 sec
> 
> on my system). I suspect the finite difference method is going to run
> 
> into serious numerical difficulties in trying to evolve the solution
> 
> over 2,500,000 steps.
> 
> 
> 
> By the way, why do use the word "fiber" to describe your approach?  I
> 
> don't understand the analogy with optical fiber, if that's what you
> 
> intended.
> 
> 
> 
> Cheers,
> 
> Krishna


That is how I feel Forth: a straightforward way to express ideas/concepts
into language of addresses&arithmetics.

A fine lecture I had some decade ago - 'Numerical heat transfer with
fluid flow' of mr. Patankar describe SIMPLE/SIMPLER algorithms to solve
that kinds of problems.
In that book is presented that a straightforward way to propagate
'boundary constraints' into domain through unidimensional slices into
domain(lines,columns) to which is applied TDMA(TriDiagonalMatrixAlgorithm).
Of course is an iterative method on domain, then on time.

A few years ago I wrote an unrolled loop TDMA in Forth with ideea to
split domain into array of pointers = "fibers" to compile them into
"TDMA calculus fibers" to obtain some speedup.  That "fibers" could
be of any form, not only rectilinear. AFAIR it run faster than
looped TDMA even on some native forth compilers (apropos, some of them 
have a limit for a definition size).

Remained to write a domain editor, then modules that describe
kind of problem to solve -> calculus of coefficients for TDMA, 
etc..
Now project is dormant yet for lack of time&motivation.

Have a nice day,
humptydumpty

P.S. If you could/have time to read that book
I think you'll be delighted.

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

#20475

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2013-03-08 16:43 -0800
Message-ID<c33a3257-71ae-46e6-a20c-3d5b0d0cc785@p5g2000yqj.googlegroups.com>
In reply to#20464
On Mar 8, 2:23 pm, humptydumpty <ouat...@gmail.com> wrote:
> On Friday, March 8, 2013 3:19:09 PM UTC+2, Krishna Myneni wrote:
> > On Mar 8, 12:32 am, humptydumpty <ouat...@gmail.com> wrote:
>
> > > On Wednesday, March 6, 2013 2:17:29 AM UTC+2, krishna...@ccreweb.org wrote:
>
> > > > On Monday, March 4, 2013 8:02:48 PM UTC-6, krishna...@ccreweb.org wrote:
>
> > > > > Here's a useful floating point calculation, solution of the 1-D diffusion equation, which may be useful for benchmarking non-native code Forth systems. A corresponding version in the higher level language, R, may be found at
>
> > > > >ftp://ccreweb.org/software/R/pde/diffusion.R
>
> > > > > KM
>
> > > > ...
>
> > > > The above R program uses a byte-code optimizer available within R. I find that the Forth version, using kforth-fast, is considerably faster than this partially optimized version.
>
> > > > Also, for reference, there is an R version of this program which uses the R package, "inline", to permit the "evolve()" function to be written in Fortran within the R program. The Fortran code is compiled and used transparently, and the speedup is enormous. Native code Forth compilers could use this hybrid R program for comparison. The link is
>
> > > >ftp://ccreweb.org/software/R/pde/diffusion-fast.R
>
> > > > KM
>
> > > Hi Krishna!
>
> > > Is possible to speed-up unrolling inner-loop
>
> > > into one word (a "calculus-fiber"):
>
> > > --
>
> > > \ syntactic sugar for simple fp arrays
>
> > > : []F   ( a n -- a' )
>
> > >         FLOATS + ;
>
> > > : []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
>
> > >     ]] FLOATS + f@ [[ ; immediate
>
> > > : []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
>
> > >     ]] FLOATS + f! [[ ; immediate
>
> > > fvariable D        0.5e D f!      \ Diffusion coefficient
>
> > > fvariable u_ext    20e  u_ext f!
>
> > > 0.01e fconstant dx
>
> > > 4e-5  fconstant dt
>
> > > D f@ dt f* dx fdup f* f/        FVARIABLE fk fk F!      \ choose dt so that fk < 0.5
>
> > > 1e0 2e0 fk f@ f* f-             FVARIABLE 1-2fk 1-2fk F!
>
> > > 10e dx f/ f>d d>s 1+ constant nx
>
> > > cr ." nx: " nx . cr
>
> > > create x      nx FLOATS allot
>
> > > create u_x_0  nx FLOATS allot
>
> > > create u_j    nx FLOATS allot
>
> > > create u_jp1  nx FLOATS allot
>
> > > : calculus-fiber:
>
> > >         :
>
> > >         nx 1- 1 DO
>
> > >                 fk           postpone LITERAL postpone F@
>
> > >                 u_j I 1+ []F postpone LITERAL postpone F@
>
> > >                 u_j I 1- []F postpone LITERAL postpone F@
>
> > >                 postpone F+ postpone F*
>
> > >                 u_j I    []F postpone LITERAL postpone F@
>
> > >                 1-2fk        postpone LITERAL postpone F@
>
> > >                 postpone F* postpone F+
>
> > >                 u_jp1 I  []F postpone LITERAL postpone F!
>
> > >         LOOP
>
> > >         postpone ;
>
> > > ;
>
> > > here calculus-fiber: fiberONE here swap - ." fiber size:" . cr
>
> > > \ see fiberONE
>
> > > : evolve2
>
> > >     u_j u_jp1 nx FLOATS move
>
> > >     0 ?DO
>
> > >         fiberONE
>
> > >         u_jp1 1 []F@ u_jp1 F!             \ b.c.: u_x(0, t) = 0
>
> > >         u_ext f@     u_jp1 nx 1- []F!   \ b.c.: u(10, t) = u_ext
>
> > >         u_jp1 u_j nx FLOATS move
>
> > >     LOOP
>
> > > ;
>
> > > \ Evolve the solution by n time steps, i.e. by an elapsed time
>
> > > \   of n*dt
>
> > > : evolve ( n -- )
>
> > >     u_j u_jp1 nx FLOATS move
>
> > >     0 ?DO
>
> > >       nx 1- 1 DO
>
> > >         u_j I 1+ []F@     u_j I []F@ 2e f* f-  u_j I 1- []F@ f+  fk f@ f*
>
> > >         u_j I    []F@ f+  u_jp1 I []F!
>
> > >       LOOP
>
> > >       u_jp1 1 []F@ u_jp1 F!             \ b.c.: u_x(0, t) = 0
>
> > >       u_ext f@     u_jp1 nx 1- []F!   \ b.c.: u(10, t) = u_ext
>
> > >       u_jp1 u_j nx FLOATS move
>
> > >     LOOP
>
> > > ;
>
> > > : init ( -- )
>
> > >     0e nx 0 DO fdup x I []F! dx f+  LOOP  fdrop
>
> > >     nx 0 DO  1e u_x_0 I []F!  LOOP
>
> > >     u_x_0 u_j nx FLOATS move
>
> > > ;
>
> > > init
>
> > > \ Compute the solution at t=1
>
> > > utime
>
> > > 25000 evolve
>
> > > utime 2swap d- cr d.
>
> > > utime
>
> > > 25000 evolve2
>
> > > utime 2swap d- cr d.
>
> > > cr .s cr bye
>
> > > --
>
> > > Output:
>
> > > ~/wrk <> gforth-fast pde.fs
>
> > > redefined DT with dt
>
> > > nx: 1001
>
> > > fiber size:87940
>
> > > 1451235
>
> > > 1116868
>
> > > <0>
>
> > > ~/wrk <> gforth pde.fs
>
> > > redefined DT with dt
>
> > > nx: 1001
>
> > > fiber size:87940
>
> > > 2036871
>
> > > 1082571
>
> > > <0>
>
> > > That "calculus-fiber" could be used to implement a TDMA.
>
> > > Then domenium of problem could be splitted into "fibers",
>
> > > then iterate domenium through "fibers".
>
> > > ( Convergence is another story. )
>
> > > Have a nice day,
>
> > > humptydumpty
>
> > Hi H.D.,
>
> > That's certainly an interesting approach for this problem, to unroll
>
> > the inner loop. A factor of two speedup is significant -- it may not
>
> > hold up for a native code compiler, where the loop overhead is likely
>
> > to be substantially lower. For raw speed improvement, I would consider
>
> > writing the inner loop as an assembler code word. That should provide
>
> > a huge improvement on indirect threaded code systems. Also, I found
>
> > that replacing the array indexing with simpler pointer arithmetic in
>
> > the inner loop also led to close to a factor of two improvement.
>
> > However, speed is only one concern with this method of calculation.
>
> > The grid size, set by dt and dx, must be small to obtain reasonable
>
> > accuracy of the solution, even for short time intervals such as t = 1.
>
> > I haven't tested the method for calculating the solution out to t=100,
>
> > as is done in the R example using the ReacTran package, which, by the
>
> > way, computes the solution across the entire grid very quickly (~5 sec
>
> > on my system). I suspect the finite difference method is going to run
>
> > into serious numerical difficulties in trying to evolve the solution
>
> > over 2,500,000 steps.
>
> > By the way, why do use the word "fiber" to describe your approach?  I
>
> > don't understand the analogy with optical fiber, if that's what you
>
> > intended.
>
> > Cheers,
>
> > Krishna
>
> That is how I feel Forth: a straightforward way to express ideas/concepts
> into language of addresses&arithmetics.
>
> A fine lecture I had some decade ago - 'Numerical heat transfer with
> fluid flow' of mr. Patankar describe SIMPLE/SIMPLER algorithms to solve
> that kinds of problems.
> In that book is presented that a straightforward way to propagate
> 'boundary constraints' into domain through unidimensional slices into
> domain(lines,columns) to which is applied TDMA(TriDiagonalMatrixAlgorithm).
> Of course is an iterative method on domain, then on time.
>
> A few years ago I wrote an unrolled loop TDMA in Forth with ideea to
> split domain into array of pointers = "fibers" to compile them into
> "TDMA calculus fibers" to obtain some speedup.  That "fibers" could
> be of any form, not only rectilinear. AFAIR it run faster than
> looped TDMA even on some native forth compilers (apropos, some of them
> have a limit for a definition size).
>
> Remained to write a domain editor, then modules that describe
> kind of problem to solve -> calculus of coefficients for TDMA,
> etc..
> Now project is dormant yet for lack of time&motivation.
>
> Have a nice day,
> humptydumpty
>
> P.S. If you could/have time to read that book
> I think you'll be delighted.

Ah, thanks. I'm afraid I don't see quite what the solution approach
is, but your reference sounds like a good one and will likely clarify
this terminology.

Krishna

[toc] | [prev] | [standalone]

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

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

csiph-web