win32forth-cvs

[George H. <geo...@us...>]

Update of /cvsroot/win32forth/win32forth-stc/demos
In directory sc8-pr-cvs9.sourceforge.net:/tmp/cvs-serv13333/win32forth-stc/demos

Added Files:
	Testpde.f matrix.f pde1.f 
Log Message:
gah:Added pde1 benchmark

--- NEW FILE: pde1.f ---
\ pde1.f  ( win32forth version)
\
\ Numerical Solution of Electrostatics Boundary-Value Problems.
\ Solve Laplace's Equation in 2 Dimensions:
\
\	D_xx u(x,y) + D_yy u(x,y) = 0
\
\ Copyright (c) 2003 Krishna Myneni, Creative Consulting for Research 
\ and Education
\
\ Provided under the terms of the GNU General Public License.
\
\ This program demonstrates a method of solving one kind of a partial 
\ differential equation (PDE) for a function u(x,y), a function
\ of the two variables x and y. In Laplace's Equation above, 
\ D_xx represents taking the second partial derivative with respect to 
\ x of u(x,y), and D_yy the second partial derivative w.r.t. y. This 
\ equation holds for the electrostatic potential u(x,y) inside
\ a charge-free two dimensional region. If we know the values of
\ u(x,y) along a boundary enclosing the region, Laplace's equation
\ may be solved to obtain the values of u(x,y) at all interior points
\ of the region. 
\
\ In this demonstration, we can setup two different bounding regions:
\
\ 1) a hollow rectangular box with voltages defined on the edges,
\
\ 2) a hollow circular region with the top half boundary at one voltage,
\      and the bottom half boundary at a second voltage.
\
\ Very thin insulators are assumed to be separating the regions which 
\ are at different potentials on the bounding region.
\ 
\ Laplace's equation is solved by an iterative application of the 
\ "mean value theorem for the electrostatic potential" (see 
\ "Classical Electrodynamics", 2nd ed, by J.D. Jackson) to each grid 
\ point inside the boundary until the solution converges. For more 
\ information on solving PDEs and boundary value problems, 
\ see "Partial differential equations for engineers and scientists", 
\ by Stanley J. Farlow, 1982, Dover. The method of solving Laplace's 
\ equation used in this example is known as Liebmann's method.
\
\
\ K. Myneni, 1998-10-23
\
\ Adapted for gforth on 2003-12-15; graphics output removed. KM
\

include matrix.f
: FROUND>S FROUND F>D D>S ;

: fmat_copy ( a1 a2 -- | copy fmatrix a1 into a2)
    over mat_size@ * dfloats cell+ cell+ cmove ;

\ Create a floating pt matrix to hold the grid values

64 constant GRIDSIZE
GRIDSIZE dup fmatrix grid
GRIDSIZE dup fmatrix last_grid	\ copy of last grid values for convergence test

\ Rectangular Region Boundary Values

100e  FCONSTANT  TOP_EDGE	\ Top edge at   100.0 V
0e    FCONSTANT  RIGHT_EDGE	\ Right edge at   0.0 V
0e    FCONSTANT  BOTTOM_EDGE	\ Bottom edge at  0.0 V
50e   FCONSTANT  LEFT_EDGE	\ Left edge at   50.0 V

: inside_rectangle? ( row col -- flag | inside rectangular boundary?)
    dup 1 > swap GRIDSIZE < AND swap
    dup 1 > swap GRIDSIZE < AND AND
;

: set_rectangular_bvs ( -- | setup the rectangular boundary values)
    GRIDSIZE 1+ 1 do  TOP_EDGE    1 i grid fmat! loop
    GRIDSIZE 1+ 1 do  RIGHT_EDGE  i GRIDSIZE grid fmat! loop
    GRIDSIZE 1+ 1 do  BOTTOM_EDGE GRIDSIZE i grid fmat! loop
    GRIDSIZE 1+ 1 do  LEFT_EDGE   i 1 grid fmat! loop
;

: init_rectangular_grid ( -- | set up the starting grid values )
    set_rectangular_bvs
    TOP_EDGE BOTTOM_EDGE RIGHT_EDGE LEFT_EDGE f+ f+ f+ 4e f/
    GRIDSIZE 1+ 1 do
      GRIDSIZE 1+ 1 do
        j i inside_rectangle? IF fdup j i grid fmat! THEN
      loop
    loop fdrop ;

\ Circular Region Boundary Values

100e  FCONSTANT  TOP_HALF	\ Top half of boundary region at 100. V
0e    FCONSTANT  BOTTOM_HALF    \ Bottom half at 0.0 V
GRIDSIZE 2 - 2/ CONSTANT RADIUS  \ Radius of boundary region

: inside_circle? ( row col -- flag | inside circular boundary? )
     GRIDSIZE 2/ - dup * swap 
     GRIDSIZE 2/ - dup * + s>f fsqrt fround>s
     RADIUS < ;
 
: set_circular_bvs ( -- | setup the circular boundary region )
    GRIDSIZE 1+ 1 do
      GRIDSIZE 1+ 1 do
        j i inside_circle? 0= IF
	  j GRIDSIZE 2/ < IF TOP_HALF ELSE BOTTOM_HALF THEN
	  j i grid fmat!
	THEN 
      LOOP
    LOOP ;

: init_circular_grid ( -- | set starting values of the grid)
    set_circular_bvs
    TOP_HALF BOTTOM_HALF f+ 2e f/
    GRIDSIZE 1+ 1 do
      GRIDSIZE 1+ 1 do
        j i inside_circle? IF fdup j i grid fmat! THEN
      loop
    loop fdrop ;
	    
defer inside?

: circ ( -- | use the two semi-circle boundary values )
    grid fmat_zero
    ['] inside_circle? is inside? 
    init_circular_grid ;

: rect ( -- | use rectangular boundary values )
    grid fmat_zero
    ['] inside_rectangle? is inside?
    init_rectangular_grid ;


: nearest@ ( i j -- f1 f2 f3 f4 | fetch the nearest neighbor grid values )
    2>R
    2R@ 1- 1 MAX grid fmat@         \ fetch left nearest neighbor
    2R@ 1+ GRIDSIZE MIN grid fmat@  \ fetch right nearest neighbor
    2R@ SWAP 1- 1 MAX SWAP grid fmat@  \ fetch up nearest neighbor
    2R> SWAP 1+ GRIDSIZE MIN SWAP grid fmat@ \ fetch down nearest neighbor
;	    	  

\ Apply the mean value theorem once to each of the interior grid values:
\   Replace each grid value with the average of the four nearest
\   neighbor values.

: iterate ( -- ) 
	GRIDSIZE 1+ 1 ?do
	  GRIDSIZE 1+ 1 ?do
	    j i inside? IF
	      j i nearest@	\ fetch four nearest neighbors
	      f+ f+ f+ 4e f/	\ take average of the four values
	      j i grid fmat!	\ store at this position
	    THEN
	  loop
	loop
;



fvariable tol	\ tolerance for solution
1e-3 tol f!

: converged? ( -- flag | test for convergence between current and last grid)
    GRIDSIZE 1+ 1 do
      GRIDSIZE 1+ 1 do
        j i inside? IF j i grid fmat@ j i last_grid fmat@ f-
	               fabs tol f@ f> IF FALSE unloop unloop EXIT THEN
		    THEN
      loop
    loop TRUE ;


\ Iterate until the solution converges to the specified tolerance 
\ at all interior points.

: solve ( -- )
	begin
          grid last_grid fmat_copy
	  iterate
	  converged?
	until
;


fvariable temp 

: grid_minmax ( -- fmin fmax | find min and max of grid values )
	1 1 grid fmat@ fdup
	GRIDSIZE 1+ 1 do
	  GRIDSIZE 1+ 1 do
	    j i grid fmat@ fswap fover fmax ( 2>r) temp f! fmin ( 2r>) temp f@
	  loop
	loop
;

: display_grid ( -- | display the grid values as a character map )
	grid_minmax
	fover f- 
	15e fswap f/	\ scale factor to scale grid value from 0 to 15
	fswap

	GRIDSIZE 1+ 1 ?do
	  GRIDSIZE 1+ 1 ?do
	    fover fover
	    j i grid fmat@ fswap f- f*
	    fround>s dup 9 >
	    if 55 + else 48 + then emit
	  loop
	  cr
	loop

	fdrop fdrop
;




rect
CR CR
.( Numerical Solution of Electrostatics Boundary-Value Problems ) CR 
GRIDSIZE dup 3 .r char x emit .    
.( grid has been setup. Type: ) CR CR
.(    rect           to use the rectangular boundary values) CR
.(    circ           to use the circular boundary values) CR
.(    solve          to find the solution) CR
.(    display_grid   to view grid as a character map) CR
CR
    

--- NEW FILE: matrix.f ---
\
\ matrix.f  (Win32Forth version)
\
\ Integer and floating point matrix manipulation routines for Forth
\   systems which use a separate floating point stack and which do NOT
\   have kForth style data typing.
\
\ Copyright (c) 1998--2002 Krishna Myneni
\
\ Revisions:
\
\	12-29-1998
\	3-29-1999   added rc>frc KM
\	12-25-1999  updated  KM
\	05-10-2000  fixed determ for singular matrix  KM
\	05-17-2000  added defining words for matrices  KM
\	08-10-2001  improved efficiency of several matrix words;
\	              about 10% faster execution in real apps  KM
\       12-09-2002  begin port to Forths with separate fp stack;
\		      changed all references to dfloats to floats  KM
\       12-14-2002  finished port to pfe and gforth. cleaned up
\                     determ and matinv by using values and the
\		      loop indexing word K	 KM
\       12-16-2002  fixed fmat_addr when size of float is not 2 cells KM
\
\ Notes:
\
\   Usage:
\	n m matrix name
\	n m fmatrix name
\
\   Examples:
\
\	3 5 matrix alpha 	( create a 3 by 5 integer matrix called alpha )
\	3 3 fmatrix beta	( create a 3 by 3 floating pt matrix, beta )

\ Memory storage format for matrices:
\   The first four bytes contains ncols and the next four bytes contains
\   nrows. The matrix data is stored next in row order.
\
\ Indexing Convention:
\	Top left element is 1, 1
\	Bottom right element is nrows, ncols
\
\ matinv and determ are based on routines from P.R. Bevington,
\   "Data Reduction and Error Analysis for the Physical Sciences".
\
\ The word K is assumed to be defined. It provides the second outer
\   loop index, as an extension of I, J, ...
\
\ Some Forths may require the following defs:
\
\ : s>f s>d d>f ;
: ?allot here swap allot ;
: float- [ 1 floats ] literal - ;

: rc_index ( n -- rc | generate an rc with a running index 1 to n )
	dup 1+ 1 do i swap loop ;

: rc_neg ( rc1 -- rc2 | negate the values in the rc )
	sp@ cell+ over 0 do dup @ negate over ! cell+ loop drop ;

: rc_dup ( rc -- rc rc | duplicate an rc on the stack )
	dup 1+ dup 0 do dup pick swap loop drop ;

: rc_max ( rc -- n | find max value in rc )
	1- dup 0> if 0 do max loop else drop then ;

: rc_min ( rc -- n | find min value in rc )
	1- dup 0> if 0 do min loop else drop then ;

: frc_index ( n -- | generate fp running index )
	    ( F: -- 1e 2e ... ne )
	dup 1+ 1 do i s>f loop ;

: frc_neg ( n -- n | negate the values in the frc )
	  ( F: f1 f2 ... fn -- -f1 -f2 ... -fn )
	dup dup 1- floats floatsp + swap
	0 do dup f@ fnegate dup f! float- loop drop ;

: frc_dup ( n -- n n | duplicate an frc on the stacks )
	  ( F: f1 f2 ... fn -- f1 f2 ... fn f1 f2 ... fn )
	dup dup 1- floats floatsp + swap 0 do dup f@ float- loop drop dup ;

: frc_max ( n --  | find max value in frc )
	  ( F: f1 f2 ... fn -- fmax )
	1- dup 0> if 0 do fmax loop else drop then ;

: frc_min ( n --  | find min value in frc )
	  ( F: f1 f2 ... fn -- fmin )
	1- dup 0> if 0 do fmin loop else drop then ;

: rc>frc ( m1 m2 ... mn n -- n | convert integer rc to frc )
	 ( F: -- f1 f2 ... fn )
	dup 0 do dup i - roll s>f loop ;

: mat_size@ ( a -- nrows ncols | gets the matrix size )
	dup cell+ @ swap @ ;

: mat_size! ( nrows ncols a -- | set up the matrix size )
	dup >r ! r> cell+ ! ;

: mat_addr ( i j a -- a2 | returns address of the i j element of a )
	>r cells swap 1- r@ @ * cells + cell+ r> + ;

: mat@ ( i j a -- n | returns the i j element of a )
	mat_addr @ ;

: mat! ( n i j a -- | store n as the i j element of a )
	mat_addr ! ;

: mat_zero ( a -- | zero all entries in matrix )
	dup mat_size@ * >r 1 1 rot mat_addr r>
	0 do 0 over ! cell+ loop drop ;

: row@ ( i a -- rc | fetch row i onto the stack as an rc )
	dup @ >r 1 swap mat_addr r> dup
	0 do over @ -rot swap cell+ swap loop nip ;

: row! ( rc i a -- | store rc as row i of matrix a )
	dup @ dup >r swap mat_addr r>
	0 do rot over ! 4 - loop 2drop ;

: col@ ( j a -- rc | fetch column j onto the stack as an rc )
	dup mat_size@ cells 2>r 1 -rot mat_addr 2r>
	swap dup >r
	0 do over @ -rot swap over + swap loop 2drop r> ;

: col! ( rc j a -- | store rc as column j of matrix a )
	dup mat_size@ cells >r dup >r -rot mat_addr r> r>
	swap
	0 do >r rot over ! r@ - r> loop 2drop drop ;

: row_swap ( i j a -- | swap rows i and j of matrix a )
	tuck 2dup 2>r 2over 2>r
	2>r row@ 2r> row@
	2r> row! 2r> row! ;

: col_swap ( i j a -- | swap columns i and j of matrix a )
	tuck 2dup 2>r 2over 2>r
	2>r col@ 2r> col@
	2r> col! 2r> col! ;

: mat. ( a -- | print out the matrix )
	dup mat_size@ 1+
	swap 1+ 1 do dup 1 do over j i rot mat@ . 9 emit loop cr loop
	2drop ;

: fmat_addr ( i j a -- a2 | returns address of the i j element of a )
	    ( F: -- )
	>r 1- floats swap 1- r@ @ * floats + r> +
	[ 2 cells ] literal + ;

: fmat@ ( i j a -- | returns the i j element of a ) ( F: -- f )
	fmat_addr f@ ;

: fmat! ( i j a -- | store f as the i j element of a ) ( F: f -- )
	fmat_addr f! ;

: fmat_zero ( a -- | zero all entries in fp matrix ) ( F: -- )
	dup mat_size@ * >r 1 1 rot fmat_addr r>
	0 do dup 0e f! float+ loop drop ;

: frow@ ( i a -- n | fetch row i of fp matrix a as an frc )
	( F: -- f1 f2 ... fn )
	dup @ >r 1 swap fmat_addr r@
	0 do dup f@ float+ loop drop r> ;

: fcol@ ( j a -- n | fetch column j of fp matrix a )
	( F: -- f1 f2 ... fn )
	dup mat_size@ floats 2>r 1 -rot fmat_addr 2r>
	swap dup >r
	0 do over f@ swap over + swap loop 2drop r> ;

: frow! ( n i a -- | store frc as row i of fp matrix a )
	( F: f1 f2 ... fn -- )
	dup @ dup >r swap fmat_addr r>
	0 do dup f! float- loop 2drop ;

: fcol! ( n j a -- | store frc as column j of fp matrix a )
	( F: f1 f2 ... fn -- )
	dup mat_size@ floats >r dup >r -rot fmat_addr r> r>
	swap
	0 do over f! dup >r - r> loop 2drop drop ;

: frow_swap ( i j a -- | interchange rows i and j for fp matrix a )
	tuck 2dup 2>r 2over 2>r
	2>r frow@ 2r> frow@
	2r> frow! 2r> frow! ;

: fcol_swap ( i j a -- | interchange columns i and j for a )
	tuck 2dup 2>r 2over 2>r
	2>r fcol@ 2r> fcol@
	2r> fcol! 2r> fcol! ;

: fmat. ( a -- | print out the fp matrix )
	cr
	dup mat_size@ 1+
	swap 1+
	1 do
	  dup
	  1 do
	    over j i rot fmat@ f. 9 emit
	  loop
	  cr
	loop
	2drop
;

\ Defining words for matrices

: matrix ( nrows ncols -- | allocate space and initialize size )
	create 2dup * cells 2 cells + ?allot mat_size! ;

: fmatrix ( nrows ncols -- | allocate space for fp matrix and initialize size )
	create 2dup * floats 2 cells + ?allot mat_size! ;



0 value norder		\ order of matrix for determ and matinv
0 value arr		\ address of array for determ and matinv
fvariable det

\ Calculate the determinant of a square floating pt matrix
\ Destroys the input matrix

: determ ( a -- | a is the fmatrix address )
	 ( F: -- fdet )
	to arr
	arr mat_size@ to norder drop
	1e det f!

	norder 1+ 1 do
	  i dup arr fmat@
	  f0= if
	    \ Find next element in row which is non-zero
	    i norder < if
	      i dup 1+
	      begin
	        2dup arr fmat@ f0=
	        over norder < and
	      while
	        1+
	      repeat
	      2dup arr fmat@
	      f0= if
	        2drop 0e fdup det f!
	        unloop exit
	      then
	      nip
	      norder 1+ i do
	        i over arr fmat@
	        i j arr fmat@
	        i   over arr fmat!
	        i j arr fmat!
	      loop
	      drop
	      det f@ fnegate det f!
	    else
	      0e fdup det f! unloop exit
	    then
	  then

\ Subtract row k from lower rows to get diagonal matrix

	  i dup dup arr fmat@ det f@ f* det f!

	  norder < if
	    norder 1+ i 1+ do
	      norder 1+ j 1+ do
	        j i arr fmat@
	        j k arr fmat@
	        k i arr fmat@
	        k dup arr fmat@
	        f/ f* f-
	        j i arr fmat!
	      loop
	    loop
	  then
	loop

	det f@
;



fvariable amax
64 1 matrix ik		\ norder can be up to 64
64 1 matrix jk
0 value kk

\ matinv computes the inverse of a symmetric matrix and returns its determinant
\ The input matrix is replaced by its inverse

: matinv ( a  -- )  ( F: -- fdet )
    dup to arr
    @ dup to norder
    dup 1 ik mat_size! 1 jk mat_size!
    1e det f!

    norder 1+ 1 do

\ Find largest element in rest of matrix

      0e amax f!

      begin
        begin
          norder 1+ i do
            norder 1+ j do
	      j i arr fmat@
	      fdup fabs amax f@ fabs
	      f>= if
	        amax f!
	        j k 1 ik mat!
	        i k 1 jk mat!
	      else
                fdrop
	      then
            loop
         loop

\ Interchange rows and columns to put amax on diagonal

         amax f@
	 f0= if 0e fdup det f! exit then

         i 1 ik mat@ i >=
       until

       i 1 ik mat@ i > if
	 i arr frow@ frc_neg
	 i 1 ik mat@ arr frow@
	 i arr frow!
	 i 1 ik mat@ arr frow!
       then

       i 1 jk mat@ i >=
     until

     i 1 jk mat@ i > if
       i arr fcol@ frc_neg
       i 1 jk mat@ arr fcol@
       i arr fcol!
       i 1 jk mat@ arr fcol!
     then

\ Accumulate elements of inverse matrix

     norder 1+ 1 do
       i j <>
       if
         i j arr fmat@ fnegate amax f@ f/
	 i j arr fmat!
       then
     loop

     norder 1+ 1 do
       norder 1+ 1 do
         j k <> if
           i k <> if
             k i arr fmat@
             j k arr fmat@
             f*
             j i arr fmat@ f+
	     j i arr fmat!
           then
         then
       loop
     loop

     norder 1+ 1 do
       i j <> if
         j i arr fmat@ amax f@ f/
         j i arr fmat!
       then
     loop

     1e amax f@ f/ i dup arr fmat!
     det f@ amax f@ f* det f!

   loop

\ Restore ordering of matrix

   norder 0 do
     norder i - to kk
     kk 1 ik mat@
     kk > if
       kk arr fcol@
       kk 1 ik mat@ arr fcol@ frc_neg
       kk arr fcol!
       kk 1 ik mat@ arr fcol!
     then
     kk 1 jk mat@
     kk > if
       kk arr frow@
       kk 1 jk mat@ arr frow@ frc_neg
       kk arr frow!
       kk 1 ik mat@ arr frow!
     then
   loop

   det f@
;














--- NEW FILE: Testpde.f ---

only forth also definitions decimal

\ defined b/float nip 0= [if] 8 constant b/float [then]

\ needs stc/float

\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\ Timing Routines                     \
\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

create TIME-BUF
here
        0 w,    \ +0  year
        0 w,    \ +2  month
        0 w,    \ +4  day of week
        0 w,    \ +6  day of month
        0 w,    \ +8  hour
        0 w,    \ +10 minute
        0 w,    \ +12 second
        0 w,    \ +14 milliseconds
here swap - constant TIME-LEN

create date$ 32 allot
create time$ 32 allot

: get-local-time ( -- )                 \ get the local computer date and time
        time-buf call GetLocalTime drop ;

: time&date     ( -- sec min hour day month year )
                get-local-time
                time-buf 12 + w@        \ seconds
                time-buf 10 + w@        \ minutes
                time-buf  8 + w@        \ hours
                time-buf  6 + w@        \ day of month
                time-buf  2 + w@        \ month of year
                time-buf      w@ ;      \ year

: .#"           ( n1 n2 -- a1 n3 )
                >r 0 <# r> 0 ?do # loop #> ;

: >date"        ( time_structure -- )
                >r 31 date$ null \ z" ddddd',' MMMM dd yyyy"
                r> null LOCALE_USER_DEFAULT
                call GetDateFormat date$ swap 1- ;

: .date         ( -- )
                get-local-time time-buf >date" type ;

: >month,day,year" ( time_structure -- )
                >r 31 date$  z" ddddd',' MMMM dd yyyy"
                r> null LOCALE_USER_DEFAULT
                call GetDateFormat date$ swap 1- ;


: .month,day,year ( -- )
                get-local-time time-buf >month,day,year" type ;

: >time"        ( time_structure -- )
                >r 31 time$ null
                r> null LOCALE_USER_DEFAULT
                call GetTimeFormat time$ swap 1- ;

: .time         ( -- )
                get-local-time time-buf >time" type ;

: >am/pm"       ( time_structure -- )
                >r 31 time$ z" h':'mmtt"
                r> null LOCALE_USER_DEFAULT
                call GetTimeFormat time$ swap 1- ;

: .am/pm        ( -- )
                get-local-time time-buf >am/pm" type ;

: ms@           ( -- ms )
                get-local-time
                time-buf
                dup   8 + w@     60 *           \ hours
                over 10 + w@ +   60 *           \ minutes
                over 12 + w@ + 1000 *           \ seconds
                swap 14 + w@ + ;                \ milli-seconds

0 value start-time

: time-reset    ( -- )
                ms@ to start-time ;

' time-reset alias timer-reset

: .elapsed      ( -- )
                ." Elapsed time: "
                ms@ start-time -
                1000 /mod
                  60 /mod
                  60 /mod 2 .#" type ." :"
                          2 .#" type ." :"
                          2 .#" type ." ."
                          3 .#" type ;

: elapse        ( -<commandline>- )
                time-reset interpret cr .elapsed ;

\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\ Other utilities                                \
\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

: ROLL          ( n1 n2 .. nk k -- n2 n3 .. nk n1 )
\  Rotate k values on the stack, bringing the deepest to the top.
                DUP>R PICK SP@ DUP CELL+ R> CELLS CELL+ MOVE DROP ;

code k          ( -- n )
                mov     -4 [ebp], eax
                mov     eax, 20 [esp]
                add     eax, 24 [esp]
                lea     ebp, -4 [ebp]
                next ;c

\ Need to set directory until fload etc work with paths.
s" demos" Prepend<home>\ "chdir

fload pde1