Mercurial > matrix-functions

comparison matrixcomp/strassenw.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 = strassenw(A, B, nmin)
2 %STRASSENW Strassen's fast matrix multiplication algorithm (Winograd variant).
3 % C = STRASSENW(A, B, NMIN), where A and B are matrices of dimension
4 % a power of 2, computes the product C = A*B.
5 % Winograd's variant of Strassen's algorithm is
6 % used recursively until dimension <= NMIN is reached,
7 % at which point standard multiplication is used.
8 % The default is NMIN = 8 (which minimizes the total number of
9 % operations).
10
11 % Reference:
12 % N. J. Higham, Accuracy and Stability of Numerical Algorithms,
13 % Second edition, Society for Industrial and Applied Mathematics,
14 % Philadelphia, PA, 2002; chap. 23.
15
16 if nargin < 3, nmin = 8; end
17
18 n = length(A);
19 if n ~= 2^( log2(n) )
20 error('The matrix dimension must be a power of 2.')
21 end
22
23 if n <= nmin
24 C = A*B;
25 else
26 m = n/2; i = 1:m; j = m+1:n;
27
28 S1 = A(j,i) + A(j,j);
29 S2 = S1 - A(i,i);
30 S3 = A(i,i) - A(j,i);
31 S4 = A(i,j) - S2;
32 S5 = B(i,j) - B(i,i);
33 S6 = B(j,j) - S5;
34 S7 = B(j,j) - B(i,j);
35 S8 = S6 - B(j,i);
36
37 M1 = strassenw( S2, S6, nmin);
38 M2 = strassenw( A(i,i), B(i,i), nmin);
39 M3 = strassenw( A(i,j), B(j,i), nmin);
40 M4 = strassenw( S3, S7, nmin);
41 M5 = strassenw( S1, S5, nmin);
42 M6 = strassenw( S4, B(j,j), nmin);
43 M7 = strassenw( A(j,j), S8, nmin);
44
45 T1 = M1 + M2;
46 T2 = T1 + M4;
47
48 C = [ M2+M3 T1+M5+M6; T2-M7 T2+M5 ];
49
50 end