Matrix multiplication

If A. is an M by N matrix, B. is an N by R matrix, and C. is an M by R matrix, then C.=A.*B. means that:

 do I = 1 to M   do K = 1 to R     C.I.K = 0     do J = 1 to N       C.I.K = C.I.K + A.I.J * B.J.K     end   end end

This calculation would involve M*N*R multiplications. The algorithm, discovered by Shmuel Winograd, uses FORMAT(N/2,,0)*M*R+(N%2)*(M+R) multiplications. It saves about half of the multiplications.

ALGORITHM

 Odd = N // 2 do K = 1 to R   Y.K = 0   do J = 1 to N % 2     JpJ = J + J; JpJm1 = JpJ - 1     Y.K = Y.K + B.JpJm1.K * B.JpJ.K   end end do I = 1 to M   X = 0   do J = 1 to N % 2     JpJ = J + J; JpJm1 = JpJ - 1     X = X + A.I.JpJ * A.I.JpJm1   end   do K = 1 to R     C.I.K = 0     if Odd then Z = A.I.N * B.N.K            else Z = 0     do J = 1 to N % 2       JpJ = J + J; JpJm1 = JpJ - 1       C.I.K = C.I.K + (A.I.JpJ + B.JpJm1.K) *,                       (A.I.JpJm1 + B.JpJ.K)     end J     C.I.K = C.I.K - X - Y.K + Z   end K end I

COMPARISON
I used the following program for Nd=10,100,200

 parse arg Nd; numeric digits Nd M = 20; N = 20; R = 20 A. = 1 / 3; B. = 1 / 3 /* fragments of code follow */ ...

Results of the tests

Running time in seconds for M = 20, N = 20, R = 20
Algorithm Nd = 10 Nd = 100 Nd = 200
definition  0.93 23.13 90.13