Mercurial > matrix-functions

comparison matrixcomp/strassen.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function C = strassen(A, B, nmin)
2 %STRASSEN Strassen's fast matrix multiplication algorithm.
3 % C = STRASSEN(A, B, NMIN), where A and B are matrices of dimension
4 % a power of 2, computes the product C = A*B.
5 % Strassen's algorithm is used recursively until dimension <= NMIN
6 % is reached, at which point standard multiplication is used.
7 % The default is NMIN = 8 (which minimizes the total number of
8 % operations).
9
10 % Reference:
11 % V. Strassen, Gaussian elimination is not optimal,
12 % Numer. Math., 13 (1969), pp. 354-356.
13
14 if nargin < 3, nmin = 8; end
15
16 n = length(A);
17 if n ~= 2^( log2(n) )
18 error('The matrix dimension must be a power of 2.')
19 end
20
21 if n <= nmin
22 C = A*B;
23 else
24 m = n/2; i = 1:m; j = m+1:n;
25 P1 = strassen( A(i,i)+A(j,j), B(i,i)+B(j,j), nmin);
26 P2 = strassen( A(j,i)+A(j,j), B(i,i), nmin);
27 P3 = strassen( A(i,i), B(i,j)-B(j,j), nmin);
28 P4 = strassen( A(j,j), B(j,i)-B(i,i), nmin);
29 P5 = strassen( A(i,i)+A(i,j), B(j,j), nmin);
30 P6 = strassen( A(j,i)-A(i,i), B(i,i)+B(i,j), nmin);
31 P7 = strassen( A(i,j)-A(j,j), B(j,i)+B(j,j), nmin);
32 C = [ P1+P4-P5+P7 P3+P5; P2+P4 P1+P3-P2+P6 ];
33 end