my loop is too slow

[Michael Haggerty]

SC Zhong <sai at dft.chem.cuhk.edu.hk> writes:

> 	for k in range(dimension):
> 		for i in range(dimension):
> 			for j in range(dimension):
> 				for h in range(dimension):
> 					a[k,i]=a[k,i]+c[k,h]*c[k,j]*\
> 						f[j,i]*f[i,h]

Always look at algorithms first, before worrying about optimization!

The problem with this code is that you have nested your matrix
multiplies, giving you an O(N**4) algorithm.  In matrix notation you
have

    A[k,i] = (C*F)[k,i] * (C*transpose(F))[k,i]

You should always compute this by computing C*F and C*transpose(F)
first (each O(N**3)) then compute the term-by-term product of those
(O(N**2)).  This algorithm will be about N/3 times faster, which in your
case is a factor of almost 50.

In this case, you should use the NumPy routine matrixmultiply(m1,m2)
to compute (C*F) and (C*transpose(F)) and the * operator for the final
term-by-term multiplication; that will shave probably another factor
of 50 from the time.  Result: 10 hours -> 14 seconds.

    # !!! Untested code follows !!!

    import Numeric
    from Numeric import matrixmultiply, transpose

    def f(C,F):
	return matrixmultiply(C,F) * matrixmultiply(C, transpose(F))

Michael

-- 
Michael Haggerty
mhagger at blizzard.harvard.edu