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