Smalll equation benchmark 
Author Message
 Smalll equation benchmark

 >
 > Libraries like BLAS and LAPACK are optimzed for solving large systems of
 > linear equations, where cache management is a major issue.
 >
 > Would anyone have recommendations for code for efficient solving of
SMALL,
 > i.e. with a size of typically 10-20, sets of linear equations, with
symmetric,
 > as well as with unsymmetric matrices of coefficients. The cost for each
set is
 > modest, but they may have to be called millions of times.

Below source may be adaptable to your problem, note it incorporates a very
fast "built-in" inverse to speed things up.
Compiling below benchmark with CVF-6.6b and running it using my
Celeron/833Mhz computer produces below output as (approx) expected..

  0.50
  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
CRUNCH SEC =  0.35

PROGRAM CRUNCH   ! benchmark utilizes least squares solution for
IMPLICIT NONE    ! random data sets,   y = b0 + b1x1 + .... + bnxn
REAL(8),ALLOCATABLE :: a(:,:), ai(:,:), b(:), temp(:), coef(:)
REAL(8),ALLOCATABLE :: xxx(:,:), yyy(:)
INTEGER,ALLOCATABLE :: ipvt(:)
REAL(8) :: c, d, x, y
INTEGER :: set, i, j, k, n, m, imax(1), clock1, clock2, rate

INTEGER :: nc = 20         ! #variables/set of independent obs
INTEGER :: nobs = 100000   ! #obs sets

n = nc+1
ALLOCATE (a(n,n), ai(n,n), b(n), temp(n), ipvt(n), &
          coef(0:nc), xxx(nc,nobs), yyy(nobs) )

CALL RANDOM_SEED()
CALL RANDOM_NUMBER (yyy)  ! get dependent obs 0->1.
CALL RANDOM_NUMBER (xxx)  ! get independent obs sets 0->1.

CALL SYSTEM_CLOCK (clock1,rate)  ! get start time  (clock1)

! accum. least squares matrices a,b
a = 0. ; b = 0.
DO set = 1,nobs
   a(1,1) = a(1,1) + 1
   y = yyy(set)
   b(1) = b(1) + y
   DO k = 1,nc
      x = xxx(k,set)
      j = k+1
      a(j,1) = a(j,1) + x
      a(1,j) = a(j,1)
      DO i = 2,nc+1
         a(i,j) = a(i,j) + x * xxx(i-1,set)
      END DO
      b(k+1) = b(k+1) + xxx(k,set) * y
   END DO
END DO

! compute inverse of "a" matrix
n = nc+1
ipvt = [ (i,i=1,n) ]
DO k = 1,n
   imax = MAXLOC(ABS(a(k:n,k)))
   m = k-1+imax(1)
   IF (m /= k) THEN
      ipvt( [m,k] ) = ipvt( [k,m] )
      a( [m,k],:) = a( [k,m],:)
   END IF
   d = 1/a(k,k)
   temp = a(:,k)
   DO j = 1,n
      c = a(k,j)*d
      a(:,j) = a(:,j)-temp*c
      a(k,j) = c
   END DO
   a(:,k) = temp*(-d)
   a(k,k) = d
END DO
ai(:,ipvt) = a    ! "a" inverse

! compute least-squares coefficients
DO k = 0,nc
   coef(k) = ai(k+1,1) * b(1)
   DO i = 1,nc
      coef(k) = coef(k) +ai(k+1,i+1) * b(i+1)
   END DO
END DO
CALL SYSTEM_CLOCK (clock2,rate)  ! get stop time

WRITE (*,'(F6.2)') coef(0)       ! bias coef. approx = 0.50 ?
WRITE (*,'(10F6.2)') coef(1:)    ! remaining coef. approx = 0.00 ?
WRITE (*,'(A,F5.2)') 'CRUNCH SEC = ', (clock2-clock1)/DBLE(rate)
WRITE (*,*) '<enter> to exit program'
READ  (*,*)
STOP
END PROGRAM



Tue, 14 Jun 2005 22:23:46 GMT  
 Smalll equation benchmark
A couple thoughts after posting my crunch source:
I recently became aware that right-side array value list can be more neatly
expressed using following syntax,  ipvt = [ 1:n ]    instead of using
implied do syntax I posted in crunch source yesterday,      ipvt = [
(i,i=1,n) ]
How about reader(s) with hotter computers posting faster results than my
0.35 sec ?

Heres the link to get source,
http://home.cfl.rr.com/davegemini/crunch.f90
 and  a link to get my CVF-compiled exe
http://home.cfl.rr.com/davegemini/crunch.exe
 to compare runtimes..

Let us know if your compiler produces a  faster runtime than my simple
CVF-compiled exe using below command line statement.

Quote:
>df crunch



Wed, 15 Jun 2005 22:32:31 GMT  
 Smalll equation benchmark

Quote:

> I recently became aware that right-side array value list can be more neatly
> expressed using following syntax,  ipvt = [ 1:n ]

I don't recognize that syntax....well, ok, I recall seeing something like
it on one of the early f8x drafts.  I also recall being disappointed
when it changed in later drafts.

Quote:
> instead of using implied do syntax I posted in crunch source yesterday,
>      ipvt = [ (i,i=1,n) ]

That one I recognize, though one part of it has been the subject of a
recent controversy (which in turn is a repeat of an old one).  I've
been a pretty strong advocate of allowing square brackets to be used
for array constructors.  So if I have my way, the second form above
will be legal f2k.  (I also personally have some sympathy for the
first form, but it isn't currently in the draft and I'd rate the
chances of it getting in as negligable).  As of the last vote on the
subject in J3, I'm with a majority on allowing the square brackets,
but there is clearly some disagreement from some people.  See the
recent thread about the German vote on the draft f2k CD.

--
Richard Maine                       |  Good judgment comes from experience;
email: my last name at host.domain  |  experience comes from bad judgment.
host: altair, domain: dfrc.nasa.gov |        -- Mark Twain



Fri, 17 Jun 2005 05:28:57 GMT  
 Smalll equation benchmark


Quote:

> Let us know if your compiler produces a  faster runtime than my simple
> CVF-compiled exe using below command line statement.
> >df crunch

Since there is no response here or by email, I assume anyone who tried this
benchmark with their non-CVF compiler came in with slower runtime vs.
running cvf-created exe on their computer..


Sat, 18 Jun 2005 22:18:29 GMT  
 Smalll equation benchmark

Quote:



> > Let us know if your compiler produces a  faster runtime than my simple
> > CVF-compiled exe using below command line statement.
> > >df crunch

> Since there is no response here or by email, I assume anyone who tried this
> benchmark with their non-CVF compiler came in with slower runtime vs.
> running cvf-created exe on their computer..

Or maybe they're all still off doing holidays instead of running
benchmarks...

J
--

****************************************************
** Facilior veniam posterius quam prius capere!   **
****************************************************
** James F. Cornwall, sole owner of all opinions  **
**  expressed in this message...                  **
****************************************************



Sun, 19 Jun 2005 01:40:42 GMT  
 
 [ 6 post ] 

 Relevant Pages 

1. Equation evaluator

2. Pell Equation (was: Continued Fractions.)

3. APL or J Solutions of Complex Partial Differential Equations

4. cubic equation in J

5. Deriving Equation out of a set of numbers

6. Equation Compiler Plugin

7. numerical equation solver

8. Equations solution

9. OLE in VisualWorks (HotDraw, Equation Editor)

10. Q: Differential Equations ?

11. Help Needed with equations in Reports

12. linear equations-we need help

 

 
Powered by phpBB® Forum Software