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

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

pde2 floating point benchmark code

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 1 of 2  [1] 2  Next page →

#20272 — pde2 floating point benchmark code

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

--

\ pde2.4th
\
\ Numerically solve the 1-D diffusion equation using the method
\ of finite differences:
\
\    u_t = D*u_xx
\
\ K. Myneni, 2013-02-26
\
\ Notes:
\
\  The function at t=0 is given by,
\
\     u(x,0) = 1
\
\  The boundary conditions are:
\
\     u(10, t) = 20  for t > 0
\     u_x(0, t) = 0
\
\ References:
\
\ 1. S.J. Farlow, Partial Differential Equations for Scientists
\    and Engineers, Dover Publications (1982); see Lesson 38.

fvariable D        0.5e D f!      \ Diffusion coefficient
fvariable u_ext    20e  u_ext f!

0.01e fconstant dx
4e-5  fconstant dt

variable nx
10e dx f/ f>d d>s 1+ nx !

create x      nx @ FLOATS allot
create u_x_0  nx @ FLOATS allot
create u_j    nx @ FLOATS allot
create u_jp1  nx @ FLOATS allot

\ syntactic sugar for simple fp arrays
: []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
    FLOATS + f@ ;
: []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
    FLOATS + f! ;

: 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

fvariable fk
D f@ dt f* dx fdup f* f/  fk f! \ choose dt so that fk < 0.5

\ 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
;


\ Compute the solution at t=1
ms@ 
25000 evolve
ms@ swap - .

[toc] | [next] | [standalone]

#20277

On Mar 5, 2:02 am, krishna.myn...@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
>
> --
>
> \ pde2.4th
> \
> \ Numerically solve the 1-D diffusion equation using the method
> \ of finite differences:
> \
> \    u_t = D*u_xx
> \
> \ K. Myneni, 2013-02-26
> \
> \ Notes:
> \
> \  The function at t=0 is given by,
> \
> \     u(x,0) = 1
> \
> \  The boundary conditions are:
> \
> \     u(10, t) = 20  for t > 0
> \     u_x(0, t) = 0
> \
> \ References:
> \
> \ 1. S.J. Farlow, Partial Differential Equations for Scientists
> \    and Engineers, Dover Publications (1982); see Lesson 38.
>
> fvariable D        0.5e D f!      \ Diffusion coefficient
> fvariable u_ext    20e  u_ext f!
>
> 0.01e fconstant dx
> 4e-5  fconstant dt
>
> variable nx
> 10e dx f/ f>d d>s 1+ nx !
>
> create x      nx @ FLOATS allot
> create u_x_0  nx @ FLOATS allot
> create u_j    nx @ FLOATS allot
> create u_jp1  nx @ FLOATS allot
>
> \ syntactic sugar for simple fp arrays
> : []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
>     FLOATS + f@ ;
> : []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
>     FLOATS + f! ;
>
> : 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
>
> fvariable fk
> D f@ dt f* dx fdup f* f/  fk f! \ choose dt so that fk < 0.5
>
> \ 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
> ;
>
> \ Compute the solution at t=1
> ms@
> 25000 evolve
> ms@ swap - .

Krishna,

Welcome back to C.L.F. I'd noted that we hadn't seen you for quite a
while, and was hoping you were okay. Good to see you back on C.L.F!

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

#20290

On Tuesday, March 5, 2013 1:23:16 AM UTC-6, Mark Wills wrote:

> 
> 
> Krishna,
> 
> 
> 
> Welcome back to C.L.F. I'd noted that we hadn't seen you for quite a
> 
> while, and was hoping you were okay. Good to see you back on C.L.F!


Hi Mark,

Thanks. I've been focused on a few other interests and family activities recently, but I've been paying attention to c.l.f. also. Haven't given up on Forth, but development had temporarily slowed. Good to see the activity maintained here!

Cheers,
Krishna

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

#20299

Using a modified version of Julian Noble's FORmula TRANslator and 
FSL-style arrays (with optional index checking):

fvalues{ D u_ext dx dt k }

f" D = 0.5    "
f" u_ext = 20 "

f" dx = 0.01  "
f" dt = 4e-5  "
f" (10/dx)+1  " f>s value nx

nx 1 floats 1array x{
nx 1 floats 1array u_i{
nx 1 floats 1array u_ip1{

f" k = D*dt/dx^2"

: init ( -- )
     nx 1- 0 DO f"   x{ i 1+ } = x{ i } + dx " LOOP
     nx 0 DO    f"    u_i{ i } = 1           " LOOP
;
init

: evolve ( n -- )
   0 DO
     nx 1- 1 DO
     f" u_ip1{ i } =
     k*( u_i{ i 1+ } - 2*u_i{ i } + u_i{ i 1- } ) + u_i{ i }"
     LOOP
   f"     u_ip1{ 0 } = u_ip1{ 1 }  "
   f" u_ip1{ nx 1- } = u_ext       "

   u_ip1{ 0 } u_i{ 0 } nx FLOATS move
   LOOP
;

timer-reset
25000 evolve
.elapsed

Runs about 33% faster (Vfx). YMMV

In my work, the ability to quickly/accurately translate conventional 
engineering equations is a productivity enhancer.

-Doug

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

#20326

On Tuesday, March 5, 2013 11:11:26 AM UTC-6, Doug Hoffman wrote:
> Using a modified version of Julian Noble's FORmula TRANslator and 
> 
> FSL-style arrays (with optional index checking):
> 
> 
> 
> fvalues{ D u_ext dx dt k }
> 
...
> Runs about 33% faster (Vfx). YMMV
> 
> 
> 
> In my work, the ability to quickly/accurately translate conventional 
> 
> engineering equations is a productivity enhancer.
> 
> 
> 
> -Doug

Hi Doug,

A metacompiler such as FTRAN or variants certainly have their place. Mostly, I think they're useful when they provide a high level of abstraction. For short, low-level number crunching, I have to say I prefer coding directly in Forth rather than using FTRAN.

I believe that your speedup is likely a result of array operations being inlined, and not as a result of the algebraic expression translator performing other optimizations (but I may be wrong!). The same effect should be achieved by making the words, []F@  and  []F! , be macros, i.e.

\ syntactic sugar for simple fp arrays
: []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
    s" FLOATS + f@" evaluate ; immediate 
: []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
    s" FLOATS + f!" evaluate ; immediate

(yes, I know it's preferable to use POSTPONE instead).

Cheers,
Krishna

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

#20339

On 3/5/13 7:11 PM, krishna.myneni@ccreweb.org wrote:
> On Tuesday, March 5, 2013 11:11:26 AM UTC-6, Doug Hoffman wrote:
>> Using a modified version of Julian Noble's FORmula TRANslator and
>>
>> FSL-style arrays (with optional index checking):
>>
>>
>>
>> fvalues{ D u_ext dx dt k }
>>
> ...
>> Runs about 33% faster (Vfx). YMMV
>>
>>
>>
>> In my work, the ability to quickly/accurately translate conventional
>>
>> engineering equations is a productivity enhancer.
>>
>>
>>
>> -Doug
>
> Hi Doug,
>
> A metacompiler such as FTRAN or variants certainly have their place. Mostly, I think they're useful when they provide a high level of abstraction. For short, low-level number crunching, I have to say I prefer coding directly in Forth rather than using FTRAN.
>
> I believe that your speedup is likely a result of array operations being inlined, and not as a result of the algebraic expression translator performing other optimizations (but I may be wrong!). The same effect should be achieved by making the words, []F@  and  []F! , be macros, i.e.
>
> \ syntactic sugar for simple fp arrays
> : []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
>      s" FLOATS + f@" evaluate ; immediate
> : []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
>      s" FLOATS + f!" evaluate ; immediate
>
> (yes, I know it's preferable to use POSTPONE instead).
>
> Cheers,
> Krishna
>

Hi Krishna,

My point was not the speed difference as it was plenty fast either way, 
IMO.  I shouldn't have mentioned it.

Rather, conversion from the original published code was much easier:

As published:
k = D*dt/dx^2
u_ip1[j]  = k*( u_i[j+1] - 2*u_i[j] + u_i[j-1] ) + u_i[j]

Julian Noble's FTRAN:
k = D*dt/dx^2
u_ip1{ i } = k*( u_i{ i 1+ } - 2*u_i{ i } + u_i{ i 1- } ) + u_i{ i }

As a bonus, the FSL-style arrays I use do index range checking (catch 
indexes too large or too small).  This could be important because I 
noticed that you used a slightly different array size and index scheme, 
which is a common thing to do.  Though the FSL indexes are expressed in 
RPN (and must use spaces), not infix, which I suppose is an inconsistency.

After using Julian's code for many years, a bug has yet to appear.

Regards,
-Doug

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

#20343

On Wednesday, March 6, 2013 3:26:55 AM UTC-6, Doug Hoffman wrote:
> On 3/5/13 7:11 PM, krishna.myneni@ccreweb.org wrote:
> 
> > On Tuesday, March 5, 2013 11:11:26 AM UTC-6, Doug Hoffman wrote:
> 
> >> Using a modified version of Julian Noble's FORmula TRANslator and
> 
> >>
> 
> >> FSL-style arrays (with optional index checking):
> 
> >>
> 
> >>
> 
> >>
> 
> >> fvalues{ D u_ext dx dt k }
> 
> >>
> 
> > ...
> 
> >> Runs about 33% faster (Vfx). YMMV
> 
> >>
> 
> >>
> 
> >>
> 
> >> In my work, the ability to quickly/accurately translate conventional
> 
> >>
> 
> >> engineering equations is a productivity enhancer.
> 
> >>
> 
> >>
> 
> >>
> 
> >> -Doug
> 
> >
> 
> > Hi Doug,
> 
> >
> 
> > A metacompiler such as FTRAN or variants certainly have their place. Mostly, I think they're useful when they provide a high level of abstraction. For short, low-level number crunching, I have to say I prefer coding directly in Forth rather than using FTRAN.
> 
> >
> 
> > I believe that your speedup is likely a result of array operations being inlined, and not as a result of the algebraic expression translator performing other optimizations (but I may be wrong!). The same effect should be achieved by making the words, []F@  and  []F! , be macros, i.e.
> 
> >
> 
> > \ syntactic sugar for simple fp arrays
> 
> > : []F@ ( a n -- ) ( F: -- r) \ ( a n -- r)
> 
> >      s" FLOATS + f@" evaluate ; immediate
> 
> > : []F! ( a n -- ) ( F: r -- ) \ ( r a n -- )
> 
> >      s" FLOATS + f!" evaluate ; immediate
> 
> >
> 
> > (yes, I know it's preferable to use POSTPONE instead).
> 
> >
> 
> > Cheers,
> 
> > Krishna
> 
> >
> 
> 
> 
> Hi Krishna,
> 
> 
> 
> My point was not the speed difference as it was plenty fast either way, 
> 
> IMO.  I shouldn't have mentioned it.
> 
> 

I'm glad you brought up the issue, though. A 30% efficiency increase is substantial! For an application like this, making the array access words be macros is something I should have done. I didn't use FSL arrays, as defined in the FSL utilities module, here since their access words are rather inefficient, but the FSL arrays are just fine for general purpose use.

> 
> Rather, conversion from the original published code was much easier:
> 
> 
> 
> As published:
> 
> k = D*dt/dx^2
> 
> u_ip1[j]  = k*( u_i[j+1] - 2*u_i[j] + u_i[j-1] ) + u_i[j]
> 
> 
> 
> Julian Noble's FTRAN:
> 
> k = D*dt/dx^2
> 
> u_ip1{ i } = k*( u_i{ i 1+ } - 2*u_i{ i } + u_i{ i 1- } ) + u_i{ i }
> 
>

For such a simple formula, the preference of one notation over another, algebraic over RPN, is mostly a matter of taste and training. The meaning is easy enough to unravel in either notation.
 
> 
> As a bonus, the FSL-style arrays I use do index range checking (catch 
> 
> indexes too large or too small).  This could be important because I 
> 
> noticed that you used a slightly different array size and index scheme, 
> 
> which is a common thing to do.  Though the FSL indexes are expressed in 
> 
> RPN (and must use spaces), not infix, which I suppose is an inconsistency.
> 
>

Can you remove the index bounds checking? I bet it will give you a dramatic increase in efficiency.
 
> 
> After using Julian's code for many years, a bug has yet to appear.
>

I'm certainly not disparaging Julian's FTRAN. I'd like to see a version of it which compiles to native code!

> 
> 
> Regards,
> 
> -Doug

Krishna

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

#20344

On 3/6/13 8:22 AM, krishna.myneni@ccreweb.org wrote:
> On Wednesday, March 6, 2013 3:26:55 AM UTC-6, Doug Hoffman wrote:

>> As published:
>> u_ip1[j]  = k*( u_i[j+1] - 2*u_i[j] + u_i[j-1] ) + u_i[j]
>>
>> Julian Noble's FTRAN:
>> u_ip1{ i } = k*( u_i{ i 1+ } - 2*u_i{ i } + u_i{ i 1- } ) + u_i{ i }
>
> For such a simple formula, the preference of one notation over another,
 > algebraic over RPN, is mostly a matter of taste and training. The meaning
 > is easy enough to unravel in either notation.

Yes.  Horses for courses.  I prefer RPN, this is Forth after all, but 
there are some classes of problems that lend well to algebraic.

>> As a bonus, the FSL-style arrays I use do index range checking (catch
>> indexes too large or too small).  This could be important because I
>> noticed that you used a slightly different array size and index scheme,
>> which is a common thing to do.

> Can you remove the index bounds checking?

Yes.

> I bet it will give you a dramatic increase in efficiency.

It is pretty slow with checking on.  But oh the headaches it can 
prevent.  IMO, all non-trivial array access should/can be done with 
checking during development.

-Doug

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

#20328

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

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

#20360

krishna.myneni@ccreweb.org writes Re: pde2 floating point benchmark code

> 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
[..]

> 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

Compare what to what?

Here is a version using Wil Baden's OPG translator. It converts to native
code (in iForth), and it can switch (without ANY change) between sfloat, 
dfloat, xfloat, and arbitrary precision. I suppose XFLOAT (80-bit) is most 
appropriate in this case.

NEEDS -xopg

.xfloat.p
\ .arbitrary.p

LET 0.5e:      GVAL D	-- Diffusion coefficient
LET 20.0:      GVAL u_ext
LET 0.01:      GVAL dx
LET 4e-5:      GVAL dt
LET D*dt/dx^2: GVAL kk	-- choose dt so that kk < 0.5

LET INTEGER(1+10/dx): =: nx

nx GARRAY x
nx GARRAY u_x_0
nx GARRAY u_i
nx GARRAY u_ip1

: INIT ( -- )
	nx 0 DO  LET x[i]=double[i]*dx:  
		 LET u_x_0[i]=1:  
		 LET u_i[i]=1:  
	   LOOP ; INIT

-- Evolve the solution by n time steps, i.e. by an elapsed time of n*dt
: EVOLVE ( n -- )
	u_i u_ip1 GCOPY
	0 ?DO	nx 1- 1 DO  LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:  LOOP
		LET  u_ip1[0]=u_ip1[1]:	\ b.c.: u_x(0, t) = 0
		LET  u_ip1[nx-1]=u_ext:	\ b.c.: u(10, t) = u_ext
		u_ip1 u_i GCOPY
	 LOOP ;

-- Compute the solution at t=1
TIMER-RESET  #25000 evolve  .ELAPSED  LET. u_i[nx-2]:

( 50% of the time is needed to copy the arrays around )
\ xfloat:     0.192 seconds elapsed. 1.9848403070438166106e+0001
\ arbitrary: 43.132 seconds elapsed. 1.9848403070438166105436153141988337094009e+0001

$01480940  : EVOLVE
$0148094A  push          $01480428 d#
$0148094F  push          $01480468 d#
$01480954  lea           rbp, [rbp -8 +] qword
$01480958  mov           [rbp 0 +] qword, $01480965 d#
$01480960  jmp           GCOPY+10 ( $0147B28A ) offset NEAR

$01480965  xor           rbx, rbx
$01480968  pop           rcx
$01480969  call          (?DO) offset NEAR
$01480973  lea           rax, [rax 0 +] qword
$01480978  push          rbx

$01480979  mov           rcx, #1000 d#
$01480980  mov           rbx, 1 d#
$01480987  call          (DO) offset NEAR
$01480991  lea           rax, [rax 0 +] qword

\ LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:

$01480998  fld           $0147FAE0 tbyte-offset
$0148099E  mov           rdi, [rbp 0 +] qword
$014809A2  lea           rdi, [rdi 1 +] qword
$014809A6  shl           rdi, 4 b#
$014809AA  fld           [rdi $064DBE90 +] tbyte
$014809B0  fld           $011BAED0 tbyte-offset
$014809B6  mov           rdi, [rbp 0 +] qword
$014809BA  shl           rdi, 4 b#
$014809BE  fld           [rdi $064DBE90 +] tbyte
$014809C4  fmulp         ST(1), ST
$014809C6  fsubp         ST(1), ST
$014809C8  mov           rdi, [rbp 0 +] qword
$014809CC  lea           rdi, [rdi -1 +] qword
$014809D0  shl           rdi, 4 b#
$014809D4  fld           [rdi $064DBE90 +] tbyte
$014809DA  faddp         ST(1), ST
$014809DC  fmulp         ST(1), ST
$014809DE  mov           rdi, [rbp 0 +] qword
$014809E2  shl           rdi, 4 b#
$014809E6  fld           [rdi $064DBE90 +] tbyte
$014809EC  faddp         ST(1), ST
$014809EE  mov           rdi, [rbp 0 +] qword
$014809F2  shl           rdi, 4 b#
$014809F6  fstp          [rdi $064DFD40 +] tbyte

$014809FC  add           [rbp 0 +] qword, 1 b#
$01480A01  add           [rbp 8 +] qword, 1 b#
$01480A06  jno           $01480998 offset NEAR
$01480A0C  add           rbp, #24 b#

\ LET  u_ip1[0]=u_ip1[1]:	\ b.c.: u_x(0, t) = 0

$01480A10  fld           $064DFD50 tbyte-offset
$01480A16  fstp          $064DFD40 tbyte-offset

\ LET  u_ip1[nx-1]=u_ext:	\ b.c.: u(10, t) = u_ext

$01480A1C  fld           $0147FA80 tbyte-offset
$01480A22  fstp          $064E3BC0 tbyte-offset

$01480A28  push          rbx
$01480A29  push          $01480468 d#
$01480A2E  push          $01480428 d#
$01480A33  lea           rbp, [rbp -8 +] qword
$01480A37  mov           [rbp 0 +] qword, $01480A44 d#
$01480A3F  jmp           GCOPY+10 ( $0147B28A ) offset NEAR
...

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

#20392

On Wednesday, March 6, 2013 3:06:09 PM UTC-6, Marcel Hendrix wrote:
> krishna.myneni@ccreweb.org writes Re: pde2 floating point benchmark code
> 
> 
> 
> > 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
> 
> [..]
> 
> 
> 
> > 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
> 
> 
> 
> Compare what to what?
> 
> 

I meant, compare the execution time of the R program's calculation (for diffusion-fast.R) with the execution time of the compiled Forth program ...

> 
> Here is a version using Wil Baden's OPG translator. It converts to native
> 
> code (in iForth), and it can switch (without ANY change) between sfloat, 
> 
> dfloat, xfloat, and arbitrary precision. I suppose XFLOAT (80-bit) is most 
> 
> appropriate in this case.
> 
>

No doubt more bits would help the truncation and roundoff errors in the calculation, with increasing t (number of time steps). There's an analytic solution to this pde, so it should be possible to compute the error vs number of time steps for the finite difference method. Also, R provides a package (ReacTran) which implements the "method of lines" for solving this, and more general variants of this type of PDE. The method of lines allows taking much larger time steps with a higher accuracy. I've done a quick visual comparison with the solution from the ReacTran calculation and the solution from the simple finite-difference method, and they agree, but I haven't determined the number of significant digits of agreement.

The method of lines calculation program may be found at

ftp://ccreweb.org/software/R/pde/pde-example.R

It outputs a rather nice color image of the spatial-temporal solution, and worth looking at if you have installed or can install R on your system.


> 
> ( 50% of the time is needed to copy the arrays around )
>

Not sure how to get around this problem. Even though the calculation might be done in-place with some extra variables, I don't know if the memory shuffling can be minimized significantly.
 
> \ xfloat:     0.192 seconds elapsed. 1.9848403070438166106e+0001
>
 
> \ arbitrary: 43.132 seconds elapsed. 1.9848403070438166105436153141988337094009e+0001
> 


Krishna

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

#20393

On Thursday, March 7, 2013 1:57:36 PM UTC+1, krishna...@ccreweb.org wrote:
> On Wednesday, March 6, 2013 3:06:09 PM UTC-6, Marcel Hendrix wrote: 
> krishna.myneni@ccreweb.org writes Re: pde2 floating point benchmark code 
>>>>> On Monday, March 4, 2013 8:02:48 PM UTC-6, krishna...@ccreweb.org wrote: >>>> Here's a useful floating point calculation, 
[..]
>> ( 50% of the time is needed to copy the arrays around ) 
> Not sure how to get around this problem. 

(not tested):
1) Indirect access: swap pointers to the data, not the data itself
2) Duplicate the innerloop code: a) buf1->buf2 b) buf2->buf1

-marcel

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

#20412

krishna.myneni@ccreweb.org writes Re: pde2 floating point benchmark code
[..]
>>>                       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
>
>> Compare what to what?

> I meant, compare the execution time of the R program's calculation 
> (for diffusion-fast.R) with the execution time of the compiled Forth program ...

Can you quote numbers? I don't have R installed.
[..]
> The method of lines calculation program may be found at

> ftp://ccreweb.org/software/R/pde/pde-example.R

> It outputs a rather nice color image of the spatial-temporal solution, and 
> worth looking at if you have installed or can install R on your system.

Could you post a picture? I can't read the R syntax to do it with my native 
plotters :-)

>> ( 50% of the time is needed to copy the arrays around )

More like 25% for the fast floating-point < 80 bits.

> Not sure how to get around this problem. Even though the calculation might 
> be done in-place with some extra variables, I don't know if the memory 
> shuffling can be minimized significantly.

Swapping pointers works.

-marcel

-- --------------------
NEEDS -xopg

.dfloat.p
\ .xfloat.p
\ .arbitrary.p

LET 0.5e:      GVAL D	-- Diffusion coefficient
LET 20.0:      GVAL u_ext
LET 0.01:      GVAL dx
LET 4e-5:      GVAL dt
LET D*dt/dx^2: GVAL kk	-- choose dt so that fk < 0.5

LET INTEGER(1+10/dx): =: nx

nx GARRAY x
nx GARRAY u_x_0
nx GARRAY u_i
nx GARRAY u_ip1

: INIT ( -- )
	nx 0 DO  LET x[i]=double[i]*dx:  
		 LET u_x_0[i]=1:  
		 LET u_i[i]=1:  
	   LOOP ; INIT

: EVOLVE ( n -- )
	u_ip1 u_i LOCALS| u_ip1 u_i |
	0 ?DO	nx 1- 1 DO  LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:  LOOP
		LET  u_ip1[0]=u_ip1[1]:	\ b.c.: u_x(0, t) = 0
		LET  u_ip1[nx-1]=u_ext:	\ b.c.: u(10, t) = u_ext
		u_ip1  u_i TO u_ip1  TO u_i
	 LOOP ;

-- Compute the solution at t=1
TIMER-RESET  #25000 evolve  CR .ELAPSED  LET. u_i[nx-2]:

\ xfloat:     0.192 seconds elapsed. 1.9848403070438166106e+0001
\  (new)      0.156 seconds elapsed. 1.9848403070438166106e+0001
\ arbitrary: 43.132 seconds elapsed. 1.9848403070438166105436153141988337094009e+0001
\  (new)     41.382 seconds elapsed. 1.9848403070438166105436153141988337094009e+0001

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

#20420

On Thursday, March 7, 2013 1:12:56 PM UTC-6, Marcel Hendrix wrote:
> krishna.myneni@ccreweb.org writes Re: pde2 floating point benchmark code
> 
> [..]
> 
> >>>                       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
> 
> >
> 
> >> Compare what to what?
> 
> 
> 
> > I meant, compare the execution time of the R program's calculation 
> 
> > (for diffusion-fast.R) with the execution time of the compiled Forth program ...
> 
> 
> 
> Can you quote numbers? I don't have R installed.
> 
> [..]
> 

I'll have to figure out the timing words used to benchmark code in R and get back with you on that. Your run times look pretty efficient (what's your hardware?).

> > The method of lines calculation program may be found at
> 
> 
> 
> > ftp://ccreweb.org/software/R/pde/pde-example.R
> 
> 
> 
> > It outputs a rather nice color image of the spatial-temporal solution, and 
> 
> > worth looking at if you have installed or can install R on your system.
> 
> 
> 
> Could you post a picture? I can't read the R syntax to do it with my native 
> 
> plotters :-)
> 
> 

The image may be found at,

http://ccreweb.org/documents/programming/pde-example.png

It's a 2-D grid, with discrete times in the interval (1, 100) along one axis, and discrete positions in the interval (0, 10) along the other axis. The value of u(x,t) is mapped to a color. One physical problem that maps to the PDE is that of a metal rod, starting out at a low uniform temperature, and then held at a fixed higher temperature at one end, and insulated at the other end. Then, u(x, t) is the temperature at a given position and time. The image shows the heat diffusing from one end to the other. The diffusion coefficient, D, represents the thermal conductivity of the material. There are, of course, other types of physical problems represented by the same 1-D diffusion equation.


> 
> >> ( 50% of the time is needed to copy the arrays around )
> 
> 
> 
> More like 25% for the fast floating-point < 80 bits.
> 
> 
> 
> > Not sure how to get around this problem. Even though the calculation might 
> 
> > be done in-place with some extra variables, I don't know if the memory 
> 
> > shuffling can be minimized significantly.
> 
> 
> 
> Swapping pointers works.
> 
> 

Ok. Thanks. My implementation of pointers actually slows down the code when I use pointer swapping instead of copying memory. However, I get a significant speedup if I don't use indexed array access, but just manipulate a pointer. The code becomes less transparent though.

> 
> -marcel
> 
>

Krishna

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

#20510

On 3/6/13 4:06 PM, Marcel Hendrix wrote:

> : EVOLVE ( n -- )
> 	u_i u_ip1 GCOPY
> 	0 ?DO	nx 1- 1 DO  LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:  LOOP

I like the infix notation for array indexes.

How are 2-dim arrays handled, like this?:

a[[i+1,i*i+2]]   ??


Do you have valid index checking?

Is your version of Baden's OPG available?

-Doug

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

#20514

Doug Hoffman <glidedog@gmail.com> writes Re: pde2 floating point benchmark code

> On 3/6/13 4:06 PM, Marcel Hendrix wrote:

>> : EVOLVE ( n -- )
>> 	u_i u_ip1 GCOPY
>> 	0 ?DO	nx 1- 1 DO  LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:  LOOP

> I like the infix notation for array indexes.

> How are 2-dim arrays handled, like this?:

> a[[i+1,i*i+2]]   ??

Like   a[i+1,i*i+2] . (Thanks for asking, I have not had a need to use this yet!
There a many efficient ways to rewrite matrix ops to vectors.)

XOPG works for integers, floats, dfloats, sfloats, complex, and arbitrary float/complex.
Therefore, for index operations (between [ and ]) the interpreter switches to
the integer XOPG tables. If the current type *is* integer, one can write 
a(i+1,i*i+2), too.

Wil Baden let the "comma" handle lists, so a function/index calculation can return
multiple values, which fits the syntax nicely.

> Do you have valid index checking?

No, the idea is to make XOPG code as fast as possible. But of course, checks can
be added. In the ideal case the compiler will warn during compilation - XOPG is 
a JIT. It would be good idea to warn that the address *IS* calculated at run-time.

> Is your version of Baden's OPG available?

I am continuously working on it, but a snapshot will be in the upcoming iForth
release. Currently, some matrix operations are being added. At some point 
type conversion will need to be implemented. For the time being there are 
functions like double() and integer() that handle the most frequent cases.

-marcel

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

#20533

On 3/10/13 12:33 PM, Marcel Hendrix wrote:
> Doug Hoffman <glidedog@gmail.com> writes Re: pde2 floating point benchmark code
>
>> On 3/6/13 4:06 PM, Marcel Hendrix wrote:
>
>>> : EVOLVE ( n -- )
>>> 	u_i u_ip1 GCOPY
>>> 	0 ?DO	nx 1- 1 DO  LET u_ip1[i]=kk*(u_i[i+1]-2*u_i[i]+u_i[i-1])+u_i[i]:  LOOP
>
>> I like the infix notation for array indexes.
>
>> How are 2-dim arrays handled, like this?:
>
>> a[[i+1,i*i+2]]   ??
>
> Like   a[i+1,i*i+2] .

OK.  So 3-dim would be a[x,y,z] and so on.


> XOPG works for integers, floats, dfloats, sfloats, complex, and arbitrary float/complex.

Julian set up FTRAN with vectored operations so it can do the same. 
It's a nice flexibility.

> Therefore, for index operations (between [ and ]) the interpreter switches to
> the integer XOPG tables. If the current type *is* integer, one can write
> a(i+1,i*i+2), too.
>
> Wil Baden let the "comma" handle lists, so a function/index calculation can return
> multiple values, which fits the syntax nicely.

Sounds the same as FTRAN, though it isn't obvious to me how to invoke it 
for array indices.


>> Do you have valid index checking?
>
> No, the idea is to make XOPG code as fast as possible.

Only use it for development.  Conditionally compile.


>> Is your version of Baden's OPG available?
>
> I am continuously working on it, but a snapshot will be in the upcoming iForth
> release.

Good.  I'll look for it when ready.


> Currently, some matrix operations are being added. At some point
> type conversion will need to be implemented. For the time being there are
> functions like double() and integer() that handle the most frequent cases.

I have found it convenient to automatically convert I,J,K to floats but 
only when not used as array indexes.

-Doug

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

#20425

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

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

#20441

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

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

#20463

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

> 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.

FORTH> see fiberONE
Flags: ANSI
$0140DE00  : fiberONE
$0140DE0A  fld           $01405130 tbyte-offset
$0140DE10  fld           $014098D0 qword-offset
$0140DE16  fld           $014098C0 qword-offset
$0140DE1C  faddp         ST(1), ST
$0140DE1E  fmulp         ST(1), ST
$0140DE20  fld           $014098C8 qword-offset
$0140DE26  fld           $01405150 tbyte-offset
$0140DE2C  fmulp         ST(1), ST
$0140DE2E  faddp         ST(1), ST
$0140DE30  fstp          $0140B828 qword-offset
$0140DE36  fld           $01405130 tbyte-offset
$0140DE3C  fld           $014098D8 qword-offset
... etcetera.

> 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 :-)

> 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?

FORTH> in
nx: 1001
fiber size:44008

0.135 seconds elapsed.
0.059 seconds elapsed. ok

| 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

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

Page 1 of 2  [1] 2  Next page →

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

csiph-web