>

> 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