Commit 56229844 authored by Rafael Laboissiere's avatar Rafael Laboissiere

New upstream version 1.0.5

parent 96153631
Name: sparsersb
Version: 1.0.2
Date: 2016-10-03
Version: 1.0.5
Date: 2017-03-29
Author: Michele Martone <michelemartone@users.sourceforge.net>
Maintainer: Michele Martone <michelemartone@users.sourceforge.net>
Title: Interface to the librsb package implementing the RSB sparse matrix format.
......
Summary of important user-visible changes for releases of the sparsersb package
===============================================================================
sparsersb-1.0.5 Release Date: 2017-03-29
===============================================================================
** Intended to be used with the latest librsb-1.2.0.
** - tests won't use quit(): might emit exception according to bug #49271
===============================================================================
sparsersb-1.0.4 Release Date: 2017-03-25
===============================================================================
** Intended to be used with the latest librsb-1.2.0.
** - test script now uses ilu() instead of obsolete luinc()
===============================================================================
sparsersb-1.0.3 Release Date: 2017-03-24
===============================================================================
** Intended to be used with the latest librsb-1.2.0.
** - "symmetric" RSB representation supported
** - improved documentation (including a few typos)
** - improved demos (demo sparsersb)
** - improved tests (test sparsersb)
** - improved error messages
===============================================================================
sparsersb-1.0.2 Release Date: 2016-10-03
===============================================================================
** Thought to be used with the latest librsb-1.2.0.
** Intended to be used with the latest librsb-1.2.0.
** - builds even if octave built with --enable-64, but limited to matrices
** which would fit when using a normal setup.
===============================================================================
sparsersb-1.0.1 Release Date: 2016-08-01
===============================================================================
** Thought to be used with the latest librsb-1.2.0.
** Intended to be used with the latest librsb-1.2.0.
** Changed sparsersb's `configure --help':
** - options to build librsb from a tarball (via configure or LIBRSB_TARBALL)
** - options to use librsb-config
......
#
# Copyright (C) 2011-2015 Michele Martone <michelemartone _AT_ users.sourceforge.net>
# Copyright (C) 2011-2017 Michele Martone <michelemartone _AT_ users.sourceforge.net>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -15,5 +15,155 @@
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
# This documentation is yet to be written.
disp "This demo is to be written"
if length(getenv("SPARSERSB_TEST")) == 0 ; pkg load sparsersb ; end
disp " ***********************************************************************"
disp "** A usage example of sparsersb. **"
disp "** A case large enough for 'sparsersb' to likely outperform 'sparse'. **"
disp "** p.s.: Invoke 'demo sparsersb' to get just a first working overview. **"
disp " ***********************************************************************"
bs=100; # block size
bo=10; # block overlap
bc=700; # block count
# bs=2; # block size
# bo=1; # block overlap
# bc=2; # block count
nr=bs+(bc-1)*(bs-bo);
nc=nr;
disp "Constructing coefficients for a sparse diagonal blocks matrix."
printf ("Will use %d blocks each wide %d and overlapping %d.\n", bc, bs, bo);
ai=[];
aj=[];
av=[];
for i=1:bc
# randomly generate block
thr=0.3;
b=rand(bs)+bc*eye(bs);
b=b+b';
[bi,bj,bv]=find(b>thr);
io=(i-1)*(bs-bo); # i offset
jo=(i-1)*(bs-bo); # j offset
ai=[ai;bi+io];
aj=[aj;bj+jo];
av=[av;bv ];
endfor
nz=length(av);
printf ("Obtained %.3e nonzeroes in a %d x %d matrix, average %.1e nonzeroes/row. \n", nz, nr, nc, nz/nr );
disp "Assembling a 'sparse' matrix..."
tic;
ao=sparse(ai,aj,av);
printf ("Assembled 'sparse' in %.1es.\n", toc);
disp "Assembling a 'sparsersb' matrix..."
tic;
ar=sparsersb(ai,aj,av);
nsb=str2num(sparsersb(ar,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
printf ("Assembled 'sparsersb' in %.1es (%d RSB blocks).\n", toc, nsb);
# Uncomment to get more information:
# printf ("RSB matrix specific info: %s.\n", ds=sparsersb(ar,"get"));
# Uncomment the following to render the RSB blocks structure to a file.
# sparsersb(ar,"render","demo_sparsersb_matrix.eps")
maxt=4;
nrhs=1;
x=ones(nc,nrhs);
y=ones(nr,nrhs);
disp " ** Testing matrix-vector multiplication ********************************"
disp "Benchmarking 'sparse' matrix-vector multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ao*x;
end
ot=dt=toc;ot/=nt;
printf ("Performed %8d 'sparse' matrix-vector multiplications in %.1es, %.2es each on average.\n", nt, dt, ot);
disp "Benchmarking 'sparsersb' matrix-vector multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ar*x;
end
rt=dt=toc;rt/=nt;
printf ("Performed %8d 'sparsersb' matrix-vector multiplications in %.1es, %.2es each on average.\n", nt, dt, rt);
printf ("So 'sparsersb' is %.2ex times as fast as 'sparse'.\n", ot/rt);
disp "Attempting autotuning 'sparsersb' matrix..."
tic;
tr=sparsersb(ar,"autotune","n",nrhs);
nnb=str2num(sparsersb(tr,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
dt=toc;
if ( nnb != nsb )
printf ("Performed autotuning in %.2es (%d -> %d RSB blocks).\n", dt, nsb, nnb);
# Uncomment to get more information:
# printf ("RSB matrix specific info: %s.\n", ds=sparsersb(tr,"get"));
disp "Benchmarking 'sparsersb' matrix-matrix multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ar*x;
end
rt=dt=toc;rt/=nt;
printf ("Performed %8d 'sparsersb' matrix-vector multiplications in %.1es, %.2es each on average.\n", nt, dt, rt);
printf ("So 'sparsersb' is %.2ex times as fast as 'sparse'.\n", ot/rt);
else
printf ("Autotuning procedure did not change the matrix.\n", dt, nsb, nnb);
end
disp " ** Testing matrix-matrix multiplication (rhs matrix is multi-vector) ***"
nrhs=5;
x=ones(nc,nrhs);
y=ones(nr,nrhs);
disp "Benchmarking 'sparse' matrix-matrix multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ao*x;
end
ot=dt=toc;ot/=nt;
printf ("Performed %8d 'sparse' matrix-matrix multiplications in %.1es, %.2es each on average.\n", nt, dt, ot);
disp "Benchmarking 'sparsersb' matrix-matrix multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ar*x;
end
rt=dt=toc;rt/=nt;
printf ("Performed %8d 'sparsersb' matrix-matrix multiplications in %.1es, %.2es each on average.\n", nt, dt, rt);
printf ("So 'sparsersb' is %.2ex times as fast as 'sparse'.\n", ot/rt);
disp "Attempting autotuning 'sparsersb' matrix..."
tic;
tr=sparsersb(ar,"autotune","n",nrhs);
nnb=str2num(sparsersb(tr,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
dt=toc;
if ( nnb != nsb )
printf ("Performed autotuning in %.2es (%d -> %d RSB blocks).\n", dt, nsb, nnb);
# Uncomment to get more information:
# printf ("RSB matrix specific info: %s.\n", ds=sparsersb(tr,"get"));
disp "Benchmarking 'sparsersb' matrix-matrix multiply..."
nt=0;tic;
while (toc < maxt)
nt++;y+=ar*x;
end
rt=dt=toc;rt/=nt;
printf ("Performed %8d 'sparsersb' matrix-matrix multiplications in %.1es, %.2es each on average.\n", nt, dt, rt);
printf ("So 'sparsersb' is %.2ex times as fast as 'sparse'.\n", ot/rt);
else
printf ("Autotuning procedure did not change the matrix.\n", dt, nsb, nnb);
end
disp " ***********************************************************************"
disp "** You can adapt this benchmark to test your matrices so to check if **"
disp "** they get multiplied faster with 'sparsersb' than with 'sparse'. **"
disp " ***********************************************************************"
#
# Copyright (C) 2011-2015 Michele Martone <michelemartone _AT_ users.sourceforge.net>
# Copyright (C) 2011-2017 Michele Martone <michelemartone _AT_ users.sourceforge.net>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -15,7 +15,7 @@
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
#
# Linear Solvers benchmarks using sparsersb.
# Linear Solvers benchmark demos using sparsersb.
#
# TODO: this file shall host some linear system solution benchmarks using sparsersb.
# It may serve as a reference point when profiling sparsersb/librsb.
......@@ -23,6 +23,12 @@
#
1; # This is a script
disp " ***********************************************************************"
disp "** Usage example of sparsersb solving linear systems with GMRES. **"
disp "** Matrices large enough for 'sparsersb' to likely outperform 'sparse'.**"
disp "** p.s.: Invoke 'demo sparsersb' to get just a first working overview. **"
disp " ***********************************************************************"
function lsb_compare(A)
n=rows(A);
maxit = n;
......@@ -30,7 +36,9 @@ b = ones (n, 1);
P = diag (diag (A));
[i,j,v]=find(sparse(A));
minres=1e-7;
printf("Solving a system of %d equations, %d nonzeroes.\n",n,nnz(A));
disp " ***********************************************************************"
printf("Solving a random system of %d equations, %d nonzeroes.\n",n,nnz(A));
disp " ***********************************************************************"
tic; Ao = sparse (i,j,v,n,n);obt=toc;
onz=nnz(Ao);
......@@ -38,7 +46,7 @@ tic; [X, FLAG, RELRES, ITER] = gmres (Ao, b, [], minres, maxit, P); odt=toc;
cs="Octave ";
onv=norm(Ao*X-b);
oRELRES=RELRES;
printf("%s took %.4f = %.4f + %.4f s and gave residual %g, flag %d, error norm %g.\n",cs,obt+odt,obt,odt,RELRES,FLAG,onv);
printf("%s took %.2e = %.2e + %.2e s and gave residual %g, flag %d, error norm %g.\n",cs,obt+odt,obt,odt,RELRES,FLAG,onv);
tic; Ar = sparsersb (i,j,v,n,n);rbt=toc;
#tic; Ar = sparsersb (Ao);rbt=toc;
......@@ -46,7 +54,7 @@ rnz=nnz(Ar);
tic; [X, FLAG, RELRES, ITER] = gmres (Ar, b, [], minres, maxit, P); rdt=toc;
cs="sparsersb";
rnv=norm(Ar*X-b);
printf("%s took %.4f = %.4f + %.4f s and gave residual %g, flag %d, error norm %g.\n",cs,rbt+rdt,rbt,rdt,RELRES,FLAG,rnv);
printf("%s took %.2e = %.2e + %.2e s and gave residual %g, flag %d, error norm %g.\n",cs,rbt+rdt,rbt,rdt,RELRES,FLAG,rnv);
if (onz != rnz)
printf("Error: seems like matrices don't match: %d vs %d nonzeroes!\n",onz,rnz);
......@@ -59,57 +67,41 @@ if (RELRES>minres ) && (oRELRES<minres )
printf("Error: sparsersb was not able to solve a system octave did (residuals: %g vs %g)!",RELRES,oRELRES);
quit(1);
else
printf("Both systems were solved, speedups for overall: %g, constructor: %g, iterations: %g.\n",(obt+odt)/(rbt+rdt),(obt)/(rbt),(odt)/(rdt));
iters=ITER(length(ITER));
printf("Both systems were solved, overall speedup: %.1ex (%.1es -> %.1es) \n (matrix construction: %.1ex, %d iterations: %.1ex).\n",(obt+odt)/(rbt+rdt),(obt+odt),(rbt+rdt),(obt)/(rbt),iters,(odt)/(rdt));
#if (obt+odt)/(rbt+rdt) > 1.0
# printf("overall: %.1ex\n",(obt+odt)/(rbt+rdt));
#end
end
printf("\n");
end
# This one is based on what Carlo De Falco posted on the octave-dev mailing list:
# This one is based on what Carlo De Falco once posted on the octave-dev mailing list:
# (he used n=1000, k=15)
n = 4;
k = 1;
A= k * eye (n) + sprandn (n, n, .8);
lsb_compare(A);
n = 100;
k = 5;
A= k * eye (n) + sprandn (n, n, .2);
lsb_compare(A);
# Toy size.
#n = 4;
#k = 1;
#A= sqrt(k) * eye (n) + sprandn (n, n, .9);
#lsb_compare(A);
n = 1000;
k = 15;
A= k * eye (n) + sprandn (n, n, .2);
# Toy size.
#n = 100;
#k = 5;
#A= sqrt(k) * eye (n) + sprandn (n, n, .8);
#lsb_compare(A);
n = 2000;
k = 1000;
A= sqrt(k) * eye (n) + sprandn (n, n, .4);
lsb_compare(A);
n = 5000;
k = 1500;
A= k * eye (n) + sprandn (n, n, .2);
lsb_compare(A);
#nx=40,ny=20;
nx0=10;ny0=10;
for xm=1:2
#for ym=1:2
for ym=1:1
nx=nx0*(2^xm),ny=ny0*(2^ym);
hx=1/(nx+1);
hy=1/(ny+1);
# a symmetric example, from andreas stahel's notes on solving linear systems
Dxx=spdiags([ones(nx,1)-2*ones(nx,1) ones(nx,1)],[-1 0 1],nx,nx)/(hx^2);
Dyy=spdiags([ones(ny,1)-2*ones(ny,1) ones(ny,1)],[-1 0 1],ny,ny)/(hy^2);
A = - kron(Dxx, speye(ny))-kron(speye(nx),Dyy);
[xx,yy]=meshgrid(linspace(hx,1-hx,nx),linspace(ny,1-hy,ny));
b=xx(:);
lsb_compare(A);
# a symmetric example, from andreas stahel's notes on solving linear systems
Dy=spdiags([ones(ny,1) ones(ny,1)],[-1 1],ny,ny)/(2*hy);
A = - kron(Dxx, speye(ny))-kron(speye(nx),Dyy) + 5*kron(speye(nx),Dy);
[xx,yy]=meshgrid(linspace(hx,1-hx,nx),linspace(ny,1-hy,ny));
b=xx(:);
A= sqrt(k) * eye (n) + sprandn (n, n, .2);
lsb_compare(A);
end
end
printf "All done."
quit(1);
disp "All done."
disp "Notice how for large matrices the matrix construction is slower..."
disp "... but multiplications are faster !"
disp " ***********************************************************************"
#
# Copyright (C) 2011-2015 Michele Martone <michelemartone _AT_ users.sourceforge.net>
# Copyright (C) 2011-2017 Michele Martone <michelemartone _AT_ users.sourceforge.net>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -15,7 +15,21 @@
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
1; # This is a script
# once complete, this program will benchmark our matrix implementation against octave's
if length(getenv("SPARSERSB_TEST")) == 0 ; pkg load sparsersb ; end
disp " ***********************************************************************"
disp "** A small 'sparse' vs 'sparsersb' benchmark. **"
disp "** On a few sample problems, tests: **"
disp "** - matrix-vector multiplication (SPMV) **"
disp "** - matrix-vector multiplication transposed (SPMV_T) **"
disp "** - sparse matrix-spares matrix multiplication (SPGEMM) **"
disp "** - and shows speedup ('RSB SPEEDUP' column) **"
disp "** p.s.: Invoke 'demo sparsersb' to get just a first working overview. **"
disp " ***********************************************************************"
disp "OP ROWS COLUMNS NONZEROES OPTIME MFLOPS RSB SPEEDUP IMPLEMENTATION"
cmt="#";
#for n_=1:6*0+1
for n_=1:6
......@@ -32,7 +46,7 @@ for ro=0:1
nz=nnz(om);
M=10^6;
if ro==1
printf("%s%s\n",cmt,"reordering with colamd");
printf("%s%s\n",cmt," reordering with colamd...");
pct=-time;
p=colamd(om);
pct+=time;
......@@ -40,9 +54,11 @@ for ro=0:1
om=om(:,p);
pat+=time;
# TODO: use an array to select/specify the different reordering algorithms
printf("%g\t%g\t(%s)\n",(nz/M)/pct,(nz/M)/pat,"mflops for pct/pat");
# printf("%g\t%g\t(%s)\n",(nz/M)/pct,(nz/M)/pat,"mflops for pct/pat");
printf("# ...colamd took %.1es (%.1e nnz/s), ",pct,nz/pct);
printf( " permutation took %.1es (%.1e nnz/s)\n",pat,nz/pat);
else
printf("%s%s\n",cmt,"no reordering");
printf("%s%s\n",cmt," testing with no reordering");
end
#bm=sparsevbr(om);
bm=sparsersb(sparse(om));
......@@ -54,26 +70,27 @@ for ro=0:1
#
bt=-time; bx=bm*b; bt+=time;
t=ot; p=["octave-",version]; mflops=(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","*",m,k,nz,mflops, p);
printf("%s\t%d\t%d\t%d\t%.1es\t%g\t%.1ex\t%s\n","SPMV ",m,k,nz,t,mflops,1 ,p);
t=bt; p=["RSB"]; mflops=(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","*",m,k,nz,mflops, p);
printf("%s\t%d\t%d\t%d\t%.1es\t%g\t%.1ex\t%s\n","SPMV ",m,k,nz,t,mflops,ot/bt,p);
## spmvt
ot=-time; ox=om.'*b; ot+=time;
#
bt=-time; bx=bm.'*b; bt+=time;
t=ot; p=["octave-",version]; mflops=(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","*",m,k,nz,mflops, p);
printf("%s\t%d\t%d\t%d\t%.1es\t%g\t%.1ex\t%s\n","SPMV_T",m,k,nz,t,mflops,1 ,p);
t=bt; p=["RSB"]; mflops=(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","*",m,k,nz,mflops, p);
printf("%s\t%d\t%d\t%d\t%.1es\t%g\t%.1ex\t%s\n","SPMV_T",m,k,nz,t,mflops,ot/bt,p);
## spgemm
ot=-time; ox=om*om; ot+=time;
#
bt=-time; bx=bm*bm; bt+=time;
t=ot; p=["octave-",version]; mflops=n*(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","OCT_SPGEMM",m,k,nz,mflops, p);
t=bt; p=["RSB"]; mflops=n*(flops/M)/t;
printf("%s\t%d\t%d\t%d\t%g\t%s\n","RSB_SPGEMM",m,k,nz,mflops, p);
t=ot; p=["octave-",version];
printf("%s\t%d\t%d\t%d\t%.1es\t\t%.1ex\t%s\n","SPGEMM",m,k,nz,t,1, p);
t=bt; p=["RSB"];
printf("%s\t%d\t%d\t%d\t%.1es\t\t%.1ex\t%s\n","SPGEMM",m,k,nz,t,ot/bt,p);
endfor
endfor
disp " ***********************************************************************"
#!/usr/bin/octave -q
#
# Copyright (C) 2011-2015 Michele Martone <michelemartone _AT_ users.sourceforge.net>
# Copyright (C) 2011-2017 Michele Martone <michelemartone _AT_ users.sourceforge.net>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -17,10 +17,21 @@
#
1; # This is a script
# a benchmark program for octave/matlab
# TODO: fix output format
# TODO: correct symmetric / hermitian matrices handling
# TODO: sound, time-and-runs-based benchmarking criteria
if length(getenv("SPARSERSB_TEST")) == 0 ; pkg load sparsersb ; end
disp " ***********************************************************************"
disp "** A small 'sparse' vs 'sparsersb' test / benchmark. **"
disp "** This is meant to be a demo, but not really an example. **"
disp "** You can invoke it supplying a Matrix Market matrix (e.g. pd.mtx). **"
disp "** Without arguments, will generate a test matrix. **"
disp "** p.s.: Invoke 'demo sparsersb' to get just a first working overview. **"
disp " ***********************************************************************"
n=10;
function printbenchline(matrixname,opname,sw,times,nnz,tottime,mxxops,bpnz,msstr)
......@@ -30,6 +41,7 @@ end
if nargin <= 0
# DGEMV benchmark
disp "** Will generate a matrix... **"
for o=1024:1024
#for o=1024:256:2048*2
m=rand(o);
......@@ -41,8 +53,8 @@ for o=1024:1024
dgemvmflops=Mflops;
printf("%d GEMV for order %d in %g secs, so %10f Mflops\n",n,o,t,n*2.0*o*o/(10^6 * t));
end
quit
endif
disp " ***********************************************************************"
else # nargin > 0
# if nargin > 0, we continue
want_sparsersb_io=1;
......@@ -55,9 +67,10 @@ end
#matrices=ls("*.mtx")';
f=1;
uc=2;
uc=2; # only 2 for the moment being.
while f<=nargin
MB=1024*1024;
printf("** Will read Matrix Market matrix file %s ...\n",f);
mmn=cell2mat(argv()(f))';
mn=strtrim(mmn');
tic();
......@@ -73,7 +86,8 @@ while f<=nargin
#if(symm=="symmetric")uc+=2;endif
if(strcmp(symm,"symmetric"))uc+=1;endif
end
wr=0 ;
disp " "
wr=0 ; # write rendering to file
if wr==1
sparsersb(sparsersb(nm),"render",[mn,"-original.eps"]);
pct=-time;
......@@ -89,14 +103,14 @@ while f<=nargin
fsz=stat(mn).size;
rt=toc();
[ia,ja,va]=find(nm);
printf("%s: %.2f MBytes read in %.4f s (%10.2f MB/s)\n",mn',fsz/MB,rt,fsz/(rt*MB));
printf("%s: %.2f MBytes read in %.4f s (%10.2f MB/s)\n",mn,fsz/MB,rt,fsz/(rt*MB));
#ia=ia'; ja=ja'; va=va';
sep=" ";
csvlstring=sprintf("#mn entries nrows ncols");
csvdstring=sprintf("%%:%s%s%d%s%d%s%d",mn,sep,entries,sep,nrows,sep,ncols);
for ski=1:uc
oppnz=1;
# FIXME: what about symmetry ?
# FIXME: what about handling symmetry ?
sparsekw="sparse";
if(ski==2)sparsekw="sparsersb";endif
if(ski==3);
......@@ -104,55 +118,64 @@ for ski=1:uc
sparsekw="sparsersb";
tic(); [nm]=sparsersb(mn); rt=toc();
sparsersb(nm,"info")
printf("%s: %.2f MBytes read by librsb in %.4f s (%10.2f MB/s)\n",mn',fsz/MB,rt,fsz/(rt*MB));
printf("%s: %.2f MBytes read by librsb in %.4f s (%10.2f MB/s)\n",mn,fsz/MB,rt,fsz/(rt*MB));
endif
if(ski==4);
nm=tril(nm);
endif
[ia,ja,va]=find(nm);
rnz=nnz(nm);
printf("benchmarking %s\n",sparsekw);
#printf("symmetry ? %s\n",issymmetric(sparse(nm)));
printf(" *** Benchmarking '%s'.\n",sparsekw);
# printf("symmetry ? %d\n",issymmetric(sparse(nm)));
mrc=rows(nm); mcc=columns(nm);
if(ski!=3);
tic();
eval(["for i=1:n; om=",sparsekw,"(ia,ja,va,mrc,mcc,\"summation\"); end"]);
printf("benchmarking %s\n",sparsekw);
printf(" *** Benchmarking '%s' instantiation from 'ia,ja,va'.\n",sparsekw);
at=toc();
#if(ski==2) tic(); nm=sparsersb(om,"autotune","N");om=nm; att=toc(); ;endif
mnz=nnz(om);
amflops=n*2.0*mnz/(10^6 * at);
printf("%s (%s) %d spBLD for %d nnz in %.4f secs, so %10.2f Mflops\n",mn',sparsekw,n,rnz,at,amflops);
amflops=n*(mnz/(10^6 * at));
printf("%s '%s' %d Instantiations for %d nnz in %.4f secs, so %10.2f nnz/s\n",mn,sparsekw,n,rnz,at,amflops);
else
mnz=rnz;
end
if(ski==2)
nsb=str2num(sparsersb(om,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
printf (" ** Assembled 'sparsersb' matrix has %d RSB blocks.\n", nsb);
endif
#rm=sparsersb(ia,ja,va);# UNFINISHED
r=linspace(1,1,size(om,1))';
v=linspace(1,1,size(om,2))';
printf(" *** Benchmarking '%s' SPMV..\n",sparsekw);
tic(); for i=1:n; r+=om *v; end; umt=toc();
UMflops=oppnz*n*2.0*mnz/(10^6 * umt);
printf("%s (%s) %d spMV for %d nnz in %.4f secs, so %10.2f Mflops\n",mn',sparsekw,n,mnz,umt, UMflops);
printf("%s '%s' %d SPMV for %d nnz in %.4f secs, so %10.2f Mflops\n",mn,sparsekw,n,mnz,umt, UMflops);
bpnz=-1; # FIXME: bytes per nonzero!
msstr="?";# FIXME: matrix structure string!
# FIXME: finish the following!
#printbenchline(mn',"spMV",sparsekw,n,mnz,umt, UMflops,bpnz,msstr);
#printbenchline(mn',"SPMV",sparsekw,n,mnz,umt, UMflops,bpnz,msstr);
#
tmp=r;r=v;v=tmp;
printf(" *** Benchmarking %s SPMV_T..\n",sparsekw);
tic(); for i=1:n; r+=om.'*v; end; tmt=toc();
TMflops=oppnz*n*2.0*mnz/(10^6 * tmt);
printf("%s (%s) %d spMVT for %d nnz in %.4f secs, so %10.2f Mflops\n",mn',sparsekw,n,mnz,tmt, TMflops);
printf("%s '%s' %d spMVT for %d nnz in %.4f secs, so %10.2f Mflops\n",mn,sparsekw,n,mnz,tmt, TMflops);
if(ski<3);
csvlstring=sprintf("%s%s",csvlstring," n at amflops umt UMflops tmt TMflops");
csvdstring=sprintf("%s%s%d%s%f%s%f%s%f%s%f%s%f%s%f",csvdstring,sep,n,sep,at,sep,amflops,sep,umt,sep,UMflops,sep,tmt,sep,TMflops);
endif
disp " "
end
++f;
printf("%s\n",csvlstring);
printf("%s\n",csvdstring);
# Uncomment following lines for benchmark-oriented output:
#printf("%s\n",csvlstring);
#printf("%s\n",csvdstring);
end
printf("benchmark terminated\n");
disp " ***********************************************************************"
endif # nargin > 0
#
# Copyright (C) 2011-2015 Michele Martone <michelemartone _AT_ users.sourceforge.net>
# Copyright (C) 2011-2017 Michele Martone <michelemartone _AT_ users.sourceforge.net>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -15,32 +15,65 @@
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
#
# This program shall attempt solution of a problem saved in the MATLAB format as for the University of Florida collection.
# This program shall attempt solution of a problem saved in the MATLAB
# format as for the University of Florida collection.
#
# e.g.: http://www.cise.ufl.edu/research/sparse/mat/Hamm/memplus.mat
# One with a problem structured as in e.g.:
# http://www.cise.ufl.edu/research/sparse/mat/Hamm/memplus.mat
# http://www.cise.ufl.edu/research/sparse/mat/Schenk_ISEI/barrier2-9.mat
#
# s=load("~/memplus.mat");
# s=load("~/barrier2-9.mat");
1; # This is a script
if length(getenv("SPARSERSB_TEST")) == 0 ; pkg load sparsersb ; end
disp " ***********************************************************************"
disp "** A usage example of sparsersb. **"
disp "** You can supply 'sparsersb' matrices to iterative method routines. **"
disp "** If the matrix is large enough, this shall secure good performance **"
disp "** of matrix-vector multiply: up to you to find method+linear system ! **"
disp "** p.s.: Invoke 'demo sparsersb' to get just a first working overview. **"
disp " ***********************************************************************"
s=load(argv(){length(argv())});
n=rows(s.Problem.A);
minres=1e-7;
maxit = n;
#maxit = 100;
minres=1e-3;
#maxit = n;
maxit = 100;
b=s.Problem.b;
#A=sparse(s.Problem.A);
A=sparsersb(s.Problem.A);
oct_A=sparse(s.Problem.A);
rsb_A=sparsersb(s.Problem.A);
printf (" **** Loaded a %d x %d matrix with %.3e nonzeroes ****\n", n, columns(s.Problem.A), nnz(s.Problem.A) );
X0=[];
RELRES=2*minres;
TOTITER=0;
M1=[]; M2=[];
M1=spdiag(A)\ones(n,1);
M2=spdiag(ones(n,1));
while RELRES >= minres ;
tic; [X1, FLAG, RELRES, ITER] = pcg (A, b, minres, maxit, M1,M2,X0); odt=toc;
RELRES
ITER;
TOTITER=TOTITER+ITER
M1=sparse(diag(s.Problem.A)\ones(n,1));
M2=sparse(diag(ones(n,1)));
nsb=str2num(sparsersb(rsb_A,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
printf (" **** The 'sparsersb' matrix consists of %d RSB blocks. ****\n", nsb);
disp " *********** Invoking pcg using a 'sparse' matrix ******************* ";
tic; [X1, FLAG, RELRES, ITER] = pcg (oct_A, b, minres, maxit, M1,M2,X0); odt=toc;
toc
disp " *********** Invoking pcg using a 'sparsersb' matrix ******************* ";
tic; [X1, FLAG, RELRES, ITER] = pcg (rsb_A, b, minres, maxit, M1,M2,X0); odt=toc;
toc
disp " ** Attempting autotuning 'sparsersb' matrix (pays off on the long run * ";
tic; rsb_A=sparsersb(rsb_A,"autotune","n",1);
toc;
nsb=str2num(sparsersb(rsb_A,"get","RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T"));
printf (" **** The 'sparsersb' matrix consists of %d RSB blocks now. **** \n", nsb);
disp " ****** Invoking pcg using a 'sparsersb' matrix (might be faster now) *** ";
tic; [X1, FLAG, RELRES, ITER] = pcg (rsb_A, b, minres, maxit, M1,M2,X0); odt=toc;
toc
X0=X1;
end
disp " *********************************************************************** ";
quit(1);
......@@ -8,13 +8,13 @@ Author: Michele MARTONE
Build / use instructions (using pkg):
To use the configure auto-detected librsb from within Octave:
> pkg -local -verbose install sparsersb-1.0.2.tar.gz
> pkg -local -verbose install sparsersb-1.0.5.tar.gz
> pkg load sparsersb
> help sparsersb
Alternatively:
tar czf sparsersb-1.0.2.tar.gz
cd sparsersb-1.0.2/src
tar czf sparsersb-1.0.5.tar.gz
cd sparsersb-1.0.5/src
./configure
make
make check
......@@ -26,7 +26,7 @@ and use it; e.g.:
or setting:
export LIBRSB_TARBALL=$HOME/librsb-1.2.0.tar.gz
before entering in Octave and building with pkg:
> pkg -local -verbose install sparsersb-1.0.2.tar.gz
> pkg -local -verbose install sparsersb-1.0.5.tar.gz
On many systems, you will have to build librsb with the PIC (-fPIC on GCC)
option or you will get link time problems.
......@@ -40,5 +40,7 @@ octave
# you can use the sparsersb function, starting with e.g.:
> help sparsersb
Check out http://librsb.sf.net for the latest librsb.
Check out http://librsb.sf.net for the latest librsb release,
http://octave.sourceforge.net/sparsersb/ for the latest sparsersb release and
http://hg.code.sf.net/p/octave/sparsersb/ for the latest repository version.
================================================================================
......@@ -10,34 +10,33 @@
-- Loadable Function: V = sparsersb (S, QS)
-- Loadable Function: sparsersb (A,"save",MTXFILENAME)
-- Loadable Function: [S, NROWS, NCOLS, NNZ, REPINFO, FIELD, SYMMETRY]
= sparsersb (MTXFILENAME, MTXTYPESTRING)
= sparsersb (MTXFILENAME[, MTXTYPESTRING])
-- Loadable Function: sparsersb (S,"render", FILENAME[, RWIDTH,
RHEIGHT])
-- Loadable Function: [O =] sparsersb (S,"autotune"[, TRANSA, NRHS,
MAXR, TMAX, TN, SF])
Create or manipulate sparse matrices using the RSB format provided
by librsb, as similarly as possible to 'sparse'.
by librsb, as similarly as possible to `sparse'.
If A is a full matrix, convert it to a sparse matrix
representation, removing all zero values in the process.
Given the integer index vectors I and J, and a 1-by-'nnz' vector of
real or complex values SV, construct the sparse matrix
'S(I(K),J(K)) = SV(K)' with overall dimensions M and N.
Given the integer index vectors I and J, and a 1-by-`nnz' vector
of real or complex values SV, construct the sparse matrix
`S(I(K),J(K)) = SV(K)' with overall dimensions M and N.
The argument 'NZMAX' is ignored but accepted for compatibility with
MATLAB and 'sparsersb'.
The argument `NZMAX' is ignored but accepted for compatibility
with MATLAB and `sparsersb'.
If M or N are not specified their values are derived from the
maximum index in the vectors I and J as given by 'M = max (I)', 'N
maximum index in the vectors I and J as given by `M = max (I)', `N
= max (J)'.