1function sstest
2%SSTEST exhaustive performance test for SSMULT.
3%
4% Example
5%   sstest
6%
7% See also ssmult, ssmultsym, ssmult_install, sstest2, mtimes.
8
9% Copyright 2007-2009, Timothy A. Davis, http://www.suitesparse.com
10
11N = [500:50:1000 1100:100:3000 3200:200:5000 ] ;
12
13% warmup for more accurate timings
14A = sparse (1) ;
15B = sparse (1) ;
16C = A*B ;
17D = ssmult(A,B) ;
18err = norm (C-D,1) ;
19if (err > 0)
20    error ('test failure') ;
21end
22clear C D
23
24titles = { ...
25    'C=A*B blue, C=B*A red, both real', ...
26    'A real, B complex', 'A complex, B real', 'both complex' } ;
27
28xlabels = { '(A random, B diagonal)', '(A random, B permutation)', ...
29    '(A random, B tridiagonal)' } ;
30
31fprintf ('\nIn the next plots, speedup is the time for MATLAB C=A*B divided\n');
32fprintf ('by the time for C=ssmult(A,B).  The X-axis is n, the dimension\n') ;
33fprintf ('of the square matrices A and B.  A is a sparse random matrix with\n');
34fprintf ('1%% nonzero values.  B is diagonal in the first row of plots,\n') ;
35fprintf ('a permutation in the 2nd row, and tridiagonal in the third.\n') ;
36fprintf ('C=A*B is in blue, C=B*A is in red.  A and B are both real in the\n') ;
37fprintf ('first column of plots, B is complex in the 2nd, A in the 3rd, and\n');
38fprintf ('both are complex in the 4th column of plots.  You will want to\n') ;
39fprintf ('maximize the figure; otherwise the text is too hard to read.\n') ;
40fprintf ('\nBe aware that in MATLAB 7.6 and later, C=A*B in MATLAB uses\n') ;
41fprintf('SSMULT (but with some additional MathWorks-specific optimizations)\n');
42fprintf ('so you are comparing nearly identical codes.\n');
43% input ('Hit enter to continue: ', 's') ;
44
45tlim = 0.1 ;
46clf ;
47
48for fig = 1:3
49
50    fprintf ('Testing C=A*B and C=B*A %s\n', xlabels {fig}) ;
51
52    T = zeros (length(N),4,4) ;
53
54    for k = 1:length(N)
55
56        n = N (k) ;
57        try
58
59            A = sprand (n,n,0.01) ;
60            if (fig == 1)
61                % B diagonal
62                B = spdiags (rand (n,1), 0, n, n) ;
63            elseif (fig == 2)
64                % B permutation
65                B = spdiags (rand (n,1), 0, n, n) ;
66                B = B (:,randperm(n)) ;
67            else
68                % B tridiagonal
69                B = spdiags (rand (n,3), -1:1, n, n) ;
70            end
71
72            for kind = 1:4
73
74                if (kind == 2)
75                    % A complex, B real
76                    A = A + 1i*sprand (A) ;
77                elseif (kind == 3)
78                    % A real, B complex
79                    A = real (A) ;
80                    B = B + 1i*sprand (B) ;
81                elseif (kind == 4)
82                    % both complex
83                    A = A + 1i*sprand (A) ;
84                    B = B + 1i*sprand (B) ;
85                end
86
87                %---------------------------------------------------------------
88                % C = A*B
89                %---------------------------------------------------------------
90
91                t1 = 0 ;
92                trials = 0 ;
93                tic
94                while (t1 < tlim)
95                    C = A*B ;
96                    trials = trials + 1 ;
97                    t1 = toc ;
98                end
99                t1 = t1 / trials ;
100
101                t2 = 0 ;
102                trials = 0 ;
103                tic
104                while (t2 < tlim)
105                    D = ssmult (A,B) ;
106                    trials = trials + 1 ;
107                    t2 = toc ;
108                end
109                t2 = t2 / trials ;
110
111                err = norm (C-D,1) ;
112                if (err > 0)
113                    error ('test failure') ;
114                end
115                clear C
116                clear D
117
118                %---------------------------------------------------------------
119                % C = B*A
120                %---------------------------------------------------------------
121
122                t3 = 0 ;
123                trials = 0 ;
124                tic
125                while (t3 < tlim)
126                    C = B*A ;
127                    trials = trials + 1 ;
128                    t3 = toc ;
129                end
130                t3 = t3 / trials ;
131
132                t4 = 0 ;
133                trials = 0 ;
134                tic
135                while (t4 < tlim)
136                    D = ssmult (B,A) ;
137                    trials = trials + 1 ;
138                    t4 = toc ;
139                end
140                t4 = t4 / trials ;
141
142                err = norm (C-D,1) ;
143                if (err > 0)
144                    error ('test failure') ;
145                end
146                clear C
147                clear D
148
149                %---------------------------------------------------------------
150
151                T (k,kind,1) = t1 ;
152                T (k,kind,2) = t2 ;
153                T (k,kind,3) = t3 ;
154                T (k,kind,4) = t4 ;
155                subplot (3,4,kind + 4*(fig-1)) ;
156                plot (N(1:k), T (1:k,kind,1) ./ T (1:k,kind,2), 'o', ...
157                      N(1:k), T (1:k,kind,3) ./ T (1:k,kind,4), 'rx', ...
158                      [N(1) n], [1 1], 'k') ;
159                xlabel (['n ' xlabels{fig}]) ;
160                ylabel ('speedup') ;
161                axis tight
162                title (titles {kind}) ;
163                drawnow
164
165            end
166
167        catch me
168            % probably because we ran out of memory ...
169            disp (me.message) ;
170            break ;
171        end
172    end
173end
174
175