octave_packages/benchmark-1.1.1/benchmark_stmm.m

   1 % Copyright (C) 2008  Jaroslav Hajek <highegg@gmail.com>
   2 %
   3 % This file is part of OctaveForge.
   4 %
   5 % OctaveForge is free software; you can redistribute it and/or modify
   6 % it under the terms of the GNU General Public License as published by
   7 % the Free Software Foundation; either version 2 of the License, or
   8 % (at your option) any later version.
   9 %
  10 % This program is distributed in the hope that it will be useful,
  11 % but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 % GNU General Public License for more details.
  14 %
  15 % You should have received a copy of the GNU General Public License
  16 % along with this software; see the file COPYING.  If not, see
  17 % <http://www.gnu.org/licenses/>.
  18 %
  19
  20 % function benchmark_stmm (n, nvec)
  21 % description:
  22 % Sparse transposed matrix-vector multiplication benchmark.
  23 % This is to test the "compound operators" feature introduced in Octave.
  24 %
  25 % arguments:
  26 % n = dimension of matrix
  27 % nvec = number of vector op repeats
  28 %
  29 % results:
  30 % time_tmm = Time for A'*B (B n^2-by-nvec matrix)
  31 % time_tmv = Time for A'*v nvec-times (v vector)
  32 % time_mtm = Time for B*A' (B nvec-by-n^2 matrix)
  33 % time_mtv = Time for v*A' nvec-times (v vector)
  34 %
  35
  36 function results = benchmark_stmm (n, nvec)
  37
  38   benchutil_default_arg ('n', 300);
  39   benchutil_default_arg ('nvec', 100);
  40
  41   benchutil_initialize (mfilename)
  42
  43   disp ('constructing sparse matrix')
  44   n = 300; % size of the grid
  45   m = n^2; % number of points
  46   X = (n-1)*rand (m, 1); Y = (n-1)*rand (m, 1);
  47   IX = ceil (X); JY = ceil (Y);
  48
  49   A = sparse(m, n^2);
  50   A = A + sparse (1:m, sub2ind ([n, n], IX  , JY  ), (IX+1-X).*(JY+1-Y), m, n^2);
  51   A = A + sparse (1:m, sub2ind ([n, n], IX+1, JY  ), (X - IX).*(JY+1-Y), m, n^2);
  52   A = A + sparse (1:m, sub2ind ([n, n], IX  , JY+1), (IX+1-X).*(Y - JY), m, n^2);
  53   A = A + sparse (1:m, sub2ind ([n, n], IX+1, JY+1), (X - IX).*(Y - JY), m, n^2);
  54
  55   v = ones (m, nvec);
  56   tic; u = A'*v; time_tmm = toc;
  57   benchutil_set_result ('time_tmm')
  58
  59   v = ones (m, 1);
  60   tic;
  61   for i=1:nvec
  62     u = A'*v;
  63   end
  64   time_tmv = toc;
  65   benchutil_set_result ('time_tmv')
  66
  67   v = ones (nvec, m);
  68   tic; u = v*A'; time_mtm = toc;
  69   benchutil_set_result ('time_mtm')
  70
  71   v = ones (1, m);
  72   tic;
  73   for i=1:nvec
  74     u = v*A';
  75   end
  76   time_mtv = toc;
  77   benchutil_set_result ('time_mtv')
  78