Provided by: librsb-dev_1.3.0.2+dfsg-6.1build1_amd64 
NAME
librsb - rsb-examples - Example programs and code
DESCRIPTION
- Examples of usage of librsb in C (<rsb.h> and <blas_sparse.h>), C++ (optional
<rsb.hpp>), Fortran (<rsb.F90> and <rsb_blas_sparse.F90>).
SYNOPSIS
Files
file assemble.cpp
C++ example based on <rsb.hpp> assembling RsbMatrix by pieces.
file autotune.cpp
C++ example based on <rsb.hpp> for performance-benchmarking a matrix from file using
RsbMatrix.tune_spmm(), for various right-hand side counts.
file bench.cpp
C++ example based on <rsb.hpp> for performance-benchmarking a matrix from file using
RsbMatrix.spmm(), for various right-hand side counts.
file build.cpp
C++ example based on <rsb.hpp> using timings for common matrix operations on
RsbMatrix: RsbMatrix.get_coo(), rsb_coo_sort(), rsb_time().
file example.cpp
C++ example based on <rsb.hpp> using RsbMatrix.spmm().
file misc.cpp
C++ example based on <rsb.hpp> showing various RsbMatrix operations.
file mtx2bin.cpp
C++ example based on <rsb.hpp> converting a Matrix Market file into a custom format
using RsbMatrix.get_coo().
file render.cpp
C++ example based on <rsb.hpp> of invoking the autotuner on an RsbMatrix matrix with
RsbMatrix.tune_spmm() and rendering it with RsbMatrix.rndr().
file span.cpp
C++ example based on <rsb.hpp> illustrating use of RsbMatrix.file_save(), and
std::span-based versions of RsbMatrix.tune_spmm(), RsbMatrix.spmv().
file twonnz.cpp
C++ example based on <rsb.hpp> measuring RsbMatrix.spmm() performance of a matrix with
only two elements; this is is effectively measuring performance of result vector
scaling.
file autotune.c
C 'RSB autotune' example program based on <rsb.h>. uses rsb_file_mtx_load(),
rsb_spmm(), rsb_tune_spmm().
file backsolve.c
C triangular solve example program. Uses rsb_spsv(),rsb_tune_spsm(). Based on <rsb.h>.
file cplusplus.cpp
C++ program using the C <rsb.h> interface. Uses rsb_tune_spmm(), rsb_spmv().
file hello-spblas.c
A first C 'hello RSB' example program using a Sparse BLAS interface and <rsb.h>. Uses
BLAS_duscr_begin(), BLAS_ussp(), BLAS_usgp(), BLAS_duscr_insert_entries(),
BLAS_duscr_end(), BLAS_dusget_element(),BLAS_dusmv(),BLAS_usds().
file hello.c
A first 'hello RSB' example C program based on <rsb.h>. Uses rsb_lib_set_opt(),
rsb_mtx_get_info_str().
file io-spblas.c
Example C program using the Sparse BLAS interface and reading from file using
rsb_blas_file_mtx_load(), BLAS_usgp(), BLAS_dusmv(), BLAS_usds().
file power.c
A toy <rsb.h>-based C program implementing the power method for computing matrix
eigenvalues. Uses rsb_spmv().
file snippets.c
Collection of C snippets of other examples. Used piecewise the documentation. Not
intended to be read as example.
file transpose.c
A toy <rsb.h>-based C program showing instantiation, transposition and other
operations on a single matrix. Uses rsb_mtx_clone(), rsb_file_mtx_save(),
rsb_file_mtx_get_dims(), rsb_file_mtx_load().
Detailed Description
Examples of usage of librsb in C (<rsb.h> and <blas_sparse.h>), C++ (optional <rsb.hpp>),
Fortran (<rsb.F90> and <rsb_blas_sparse.F90>).
The following fully working example programs illustrate correct ways of using the library.
• First example in C, using <rsb.h>: examples_hello_c.
• First example in C, using <blas_sparse.h>: examples_hello_spblas_c.
• Autotuning example in C, using <rsb.h>: examples_autotune_c.
• I/O example in C, using <blas_sparse.h>: examples_io_spblas_c.
• Example transposing a matrix in C, using <rsb.h> in C: examples_transpose_c.
• Example showing the power method in C, using <rsb.h> in C: examples_power_c.
• Example in Fortran, using <rsb_blas_sparse.F90>: examples_fortran_F90.
• Example in Fortran, using <rsb.F90>: examples_fortran_rsb_fi_F90.
• Example in C, using <rsb.h>: examples_backsolve_c.
• Misc example snippets in C, using <rsb.h>: examples_snippets_c.
• Benchmark invocation from shell script: examples_bench_sh.
Once installed librsb, the script displayed here (examples/make.sh) should be sufficient
to build these examples:
#!/bin/bash
# Example script to build the librsb example programs.
# Uses librsb-config in the $PATH for build flags.
# Environment-provided $LIBRSB_CONFIG override that.
set -e
set -x
srcdir=${srcdir:-`pwd`}
builddir=${builddir:-`pwd`}
prefix="/usr"
exec_prefix="${prefix}"
bindir="${exec_prefix}/bin"
if test -z "${LIBRSB_CONFIG}"; then
export PATH="${bindir}:$PATH"
fi
LIBRSB_CONFIG=${LIBRSB_CONFIG:-librsb-config}
PKG_CONFIG=pkg-config
WANT_PKGCONFIG="yes"
WANT_CLEANUP=${WANT_CLEANUP:-false}
if test x"yes" == x"yes" ; then
CXX="`${LIBRSB_CONFIG} --cxx`"
if test x"rsblib" != x"" -a x"${CXX}" != x"" ; then
for s in ${srcdir}/*.cpp
do
p=${builddir}/`basename ${s/.cpp/}`
rm -f $p
CXXFLAGS=`${LIBRSB_CONFIG} --cxxflags --I_opts`
LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
LINK=`${LIBRSB_CONFIG} --link`
o="${p}.o"
ccmd="$CXX $CXXFLAGS -c $s -o $o"
lcmd="$LINK $o $LDFLAGS -o $p"
echo "$ccmd && $lcmd"
( $ccmd && $lcmd )
${WANT_CLEANUP} && rm -f "$p"
# one may use a single command, but that's error-prone (may miss libraries):
#cmd="$CXX $CXXFLAGS $s $LDFLAGS -o $p"
#echo $cmd
#$cmd
done
fi
fi
if test x"yes" == x"yes" ; then
for s in ${srcdir}/*.c
do
p=`basename ${s/.c/}`
if test $p == hello-spblas -a x"yes" != x"yes" ; then continue; fi
if test $p == io-spblas -a x"yes" != x"yes" ; then continue; fi
rm -f $p
CFLAGS=`${LIBRSB_CONFIG} --I_opts --cppflags`
LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
CC=`${LIBRSB_CONFIG} --cc`
LINK=`${LIBRSB_CONFIG} --link`
o="${p}.o"
ccmd="$CC $CFLAGS -c $s -o $o"
lcmd="$LINK $o $LDFLAGS -o $p"
echo "$ccmd && $lcmd"
( $ccmd && $lcmd )
${WANT_CLEANUP} && rm -f "$p"
# one may use a single command, but that's error-prone (may miss libraries):
#cmd="$CC $CFLAGS $s $LDFLAGS -o $p"
#echo $cmd
#$cmd
if test x"${WANT_PKGCONFIG}" != x"no" ; then
CFLAGS=`${PKG_CONFIG} --cflags librsb`
LIBS=`${PKG_CONFIG} --libs --static librsb`
ccmd="$CC $CFLAGS -c $s -o $o"
lcmd="$LINK $o $LIBS -o $p"
${WANT_CLEANUP} && rm -f "$p"
echo "$ccmd && $lcmd"
( $ccmd && $lcmd )
fi
done
fi
if test x"yes" == x"yes" ; then
if test x"yes" = x"yes" ; then
FP=${srcdir}/fortran_rsb_fi.F90
if test x"yes" = x"yes" ; then
FP+=\ ${srcdir}/fortran.F90
fi
# activated if you have built the Fortran modules and installed them in the right path.
for s in $FP
do
p=`basename ${s/.F90/}`
rm -f $p
FCFLAGS=`${LIBRSB_CONFIG} --I_opts --fcflags`
LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
FC=`${LIBRSB_CONFIG} --fc`
LINK=`${LIBRSB_CONFIG} --link`
o="${p}.o"
FCLIBS=`${LIBRSB_CONFIG} --fclibs`
ccmd="$FC $FCFLAGS -c $s -o $o"
lcmd="$LINK $o $LDFLAGS $FCLIBS -o $p"
echo "$ccmd && $lcmd"
( $ccmd && $lcmd )
${WANT_CLEANUP} && rm -f "$p"
done
fi
fi
echo " [*] done building examples!"
examples/hello.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
\ingroup rsb-examples
@file
@author Michele Martone
@brief A first "hello RSB" example C program based on <rsb.h>.
Uses #rsb_lib_set_opt(), #rsb_mtx_get_info_str().
\include hello.c
*/
#include <rsb.h> /* librsb header to include */
#include <stdio.h> /* printf() */
int main(const int argc, char * const argv[])
{
/*!
A Hello-RSB program.
This program shows how to use the rsb.h interface correctly to:
- initialize the library using #rsb_lib_init()
- set library options using #rsb_lib_set_opt()
- revert such changes
- allocate (build) a single sparse matrix in the RSB format
using #rsb_mtx_alloc_from_coo_const()
- prints information obtained via #rsb_mtx_get_info_str()
- multiply the matrix times a vector using #rsb_spmv()
- deallocate the matrix using #rsb_mtx_free()
- finalize the library using #rsb_lib_exit()
In this example, we use #RSB_DEFAULT_TYPE as matrix type.
This type depends on what was configured at library build time.
* */
const rsb_blk_idx_t bs = RSB_DEFAULT_BLOCKING;
const rsb_blk_idx_t brA = bs, bcA = bs;
const RSB_DEFAULT_TYPE one = 1;
const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
const rsb_nnz_idx_t nnzA = 4; /* matrix nonzeroes count */
const rsb_coo_idx_t nrA = 3; /* matrix rows count */
const rsb_coo_idx_t ncA = 3; /* matrix columns count */
/* nonzero row indices coordinates: */
const rsb_coo_idx_t IA[] = {0,1,2,2};
/* nonzero column indices coordinates: */
const rsb_coo_idx_t JA[] = {0,1,2,2};
const RSB_DEFAULT_TYPE VA[] = {11,22,32,1};/* values of nonzeroes */
RSB_DEFAULT_TYPE X[] = { 0, 0, 0 }; /* X vector's array */
const RSB_DEFAULT_TYPE B[] = { -1, -2, -5 }; /* B vector's array */
char ib[200];
struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
rsb_err_t errval = RSB_ERR_NO_ERROR;
printf("Hello, RSB!\n");
printf("Initializing the library...\n");
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) !=
RSB_ERR_NO_ERROR)
{
printf("Error initializing the library!\n");
goto err;
}
printf("Correctly initialized the library.\n");
printf("Attempting to set the"
" RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE library option.\n");
{
rsb_int_t evi=1;
/* Setting a single optional library parameter. */
errval = rsb_lib_set_opt(
RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE, &evi);
if(errval != RSB_ERR_NO_ERROR)
{
char errbuf[256];
rsb_strerror_r(errval,&errbuf[0],sizeof(errbuf));
printf("Failed setting the"
" RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE"
" library option (reason string:\n%s).\n",errbuf);
if(errval&RSB_ERRS_UNSUPPORTED_FEATURES)
{
printf("This error may be safely ignored.\n");
}
else
{
printf("Some unexpected error occurred!\n");
goto err;
}
}
else
{
printf("Setting back the "
"RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE"
" library option.\n");
evi = 0;
errval = rsb_lib_set_opt(RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE,
&evi);
errval = RSB_ERR_NO_ERROR;
}
}
mtxAp = rsb_mtx_alloc_from_coo_const(
VA,IA,JA,nnzA,typecode,nrA,ncA,brA,bcA,
RSB_FLAG_NOFLAGS /* default format will be chosen */
|RSB_FLAG_DUPLICATES_SUM/* duplicates will be summed */
,&errval);
if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
{
printf("Error while allocating the matrix!\n");
goto err;
}
printf("Correctly allocated a matrix.\n");
printf("Summary information of the matrix:\n");
/* print out the matrix summary information */
rsb_mtx_get_info_str(mtxAp,"RSB_MIF_MATRIX_INFO__TO__CHAR_P",
ib,sizeof(ib));
printf("%s",ib);
printf("\n");
if((errval =
rsb_spmv(RSB_TRANSPOSITION_N,&one,mtxAp,B,1,&one,X,1))
!= RSB_ERR_NO_ERROR )
{
printf("Error performing a multiplication!\n");
goto err;
}
printf("Correctly performed a SPMV.\n");
rsb_mtx_free(mtxAp);
printf("Correctly freed the matrix.\n");
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
printf("Error finalizing the library!\n");
goto err;
}
printf("Correctly finalized the library.\n");
printf("Program terminating with no error.\n");
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.\n");
return EXIT_FAILURE;
}
examples/hello-spblas.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
\ingroup rsb-examples
@file
@author Michele Martone
@brief A first C "hello RSB" example program using
a Sparse BLAS interface and <rsb.h>.
Uses #BLAS_duscr_begin(), #BLAS_ussp(), #BLAS_usgp(),
#BLAS_duscr_insert_entries(), #BLAS_duscr_end(),
#BLAS_dusget_element(),#BLAS_dusmv(),#BLAS_usds().
\include hello-spblas.c
*/
#include <rsb.h> /* for rsb_lib_init */
#include <blas_sparse.h> /* Sparse BLAS on the top of librsb */
#include <stdio.h> /* printf */
int main(const int argc, char * const argv[])
{
/*!
* A Hello/Sparse BLAS program.
*
* This program shows how to use the blas_sparse.h
* interface correctly to:
*
* - initialize the library using #rsb_lib_init()
* - allocate (build) a single sparse matrix in the RSB
* format using #BLAS_duscr_begin()/#BLAS_duscr_insert_entries()
* /#BLAS_duscr_end()
* - extract one matrix element with #BLAS_dusget_element()
* - multiply the matrix times a vector using #BLAS_dusmv()
* - deallocate the matrix using #BLAS_usds()
* - finalize the library using
* #rsb_lib_exit(#RSB_NULL_EXIT_OPTIONS)
*/
#ifndef RSB_NUMERICAL_TYPE_DOUBLE
printf("'double' type configured out."
" Please reconfigure the library with it and recompile.\n");
return EXIT_SUCCESS;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
blas_sparse_matrix A = blas_invalid_handle; /* handle for A */
const int nnz = 4; /* number of nonzeroes of matrix A */
const int nr = 3; /* number of A's rows */
const int nc = 3; /* number of A's columns */
/* A's nonzero elements row indices (coordinates): */
#ifdef RSB_WANT_LONG_IDX_TYPE
const int64_t IA[] = { 0, 1, 2, 2 };
#else /* RSB_WANT_LONG_IDX_TYPE */
const int IA[] = { 0, 1, 2, 2 };
#endif /* RSB_WANT_LONG_IDX_TYPE */
/* A's nonzero elements column indices (coordinates): */
#ifdef RSB_WANT_LONG_IDX_TYPE
const int64_t JA[] = { 0, 1, 0, 2 };
#else /* RSB_WANT_LONG_IDX_TYPE */
const int JA[] = { 0, 1, 0, 2 };
#endif /* RSB_WANT_LONG_IDX_TYPE */
/* A's nonzero values (matrix coefficients): */
double VA[] = { 11.0, 22.0, 13.0, 33.0 };
/* the X vector's array: */
double X[] = { 0.0, 0.0, 0.0 };
/* the B vector's array: */
const double B[] = { -1.0, -2.0, -2.0 };
/* the (known) result array: */
const double AB[] = { 11.0+26.0, 44.0, 66.0+13.0 };
/* rsb error variable: */
rsb_err_t errval = RSB_ERR_NO_ERROR;
int i;
printf("Hello, RSB!\n");
/* initialize the library */
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly initialized the library.\n");
/* initialize a matrix descriptor */
A = BLAS_duscr_begin(nr,nc);
if( A == blas_invalid_handle )
{
goto err;
}
/* specify properties (e.g.: symmetry)*/
if( BLAS_ussp(A,blas_lower_symmetric) != 0 )
{
goto err;
}
/* get properties (e.g.: symmetry) */
if( BLAS_usgp(A,blas_lower_symmetric) != 1 )
{
printf("Symmetry property non set ?!\n");
goto err;
}
/* insert the nonzeroes (here, all at once) */
if( BLAS_duscr_insert_entries(A, nnz, VA, IA, JA)
== blas_invalid_handle)
{
goto err;
}
/* finalize (allocate) the matrix build */
if( BLAS_duscr_end(A) == blas_invalid_handle )
{
goto err;
}
printf("Correctly allocated a matrix.\n");
VA[0] = 0.0;
if( BLAS_dusget_element(A, IA[0], JA[0], &VA[0]) )
{
goto err;
}
/* a check */
if( VA[0] != 11.0 )
{
goto err;
}
/* compute X = X + (-1) * A * B */
if(BLAS_dusmv(blas_no_trans,-1,A,B,1,X,1))
{
goto err;
}
for( i = 0 ; i < nc; ++i )
if( X[i] != AB[i] )
{
printf("Computed SPMV result seems wrong. Terminating.\n");
goto err;
}
printf("Correctly performed a SPMV.\n");
/* deallocate matrix A */
if( BLAS_usds(A) )
{
goto err;
}
printf("Correctly freed the matrix.\n");
/* finalize the library */
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly finalized the library.\n");
printf("Program terminating with no error.\n");
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.\n");
return EXIT_FAILURE;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}
examples/autotune.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
\ingroup rsb-examples
@file
@author Michele Martone
@brief C "RSB autotune" example program based on <rsb.h>.
uses #rsb_file_mtx_load(), rsb_spmm(), rsb_tune_spmm().
\include autotune.c
*/
#include <rsb.h> /* librsb header to include */
#include <stdio.h> /* printf() */
#include <ctype.h> /* isdigit() */
#include <stdlib.h> /* atoi() */
/* #include "rsb_internals.h" */
static int tune_from_file(char * const filename, rsb_int_t wvat)
{
struct rsb_mtx_t *mtxMp = NULL;
/* spmv specific variables */
const void * alphap = NULL; // equivalent to 1
const void * betap = NULL; // equivalent to 1
rsb_flags_t order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER;
const rsb_coo_idx_t nrhs = 2; /* number of right hand sides */
rsb_trans_t transA = RSB_TRANSPOSITION_N; /* transposition */
rsb_nnz_idx_t ldB = 0;
rsb_nnz_idx_t ldC = 0;
/* misc variables */
rsb_err_t errval = RSB_ERR_NO_ERROR;
rsb_time_t dt;
char ib[200];
const char*is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
/* misc variables */
/* input autotuning variables */
rsb_int_t oitmax = 1 /*15*/; /* auto-tune iterations */
rsb_time_t tmax = 0.1; /* time per autotune operation */
/* output autotuning variables */
rsb_flags_t flagsA = RSB_FLAG_NOFLAGS;
/* int ione = 1; */
rsb_type_t typecodea [] = RSB_MATRIX_TYPE_CODES_ARRAY;
int typecodei;
errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS);
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
/*
errval = rsb_lib_set_opt(RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE, &ione);
*/
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
printf("Loading matrix from file \"%s\".\n",filename);
mtxMp = rsb_file_mtx_load(filename, flagsA, typecodea[0], &errval);
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
for( typecodei = 0 ; typecodei < RSB_IMPLEMENTED_TYPES; ++typecodei )
{
rsb_type_t typecode = typecodea[typecodei];
struct rsb_mtx_t *mtxAp = NULL;
struct rsb_mtx_t *mtxOp = NULL;
rsb_real_t sf = 0.0;
rsb_int_t tn = 0;
sf = 0.0;
tn = 0;
printf("Considering %c clone.\n",typecode);
errval = rsb_mtx_clone(&mtxAp, typecode, transA, NULL, mtxMp,
flagsA);
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
printf("Base matrix:\n");
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s\n\n",ib);
dt = -rsb_time();
errval = rsb_tune_spmm(NULL, &sf, &tn, oitmax, tmax, transA,
alphap, mtxAp, nrhs, order, NULL, ldB, betap, NULL, ldC);
dt += rsb_time();
if(tn == 0)
printf("After %lfs, autotuning routine did not find a better"
" threads count configuration.\n",dt);
else
printf("After %lfs, thread autotuning declared speedup of %lg x,"
" when using threads count of %d.\n",dt,sf,tn);
printf("\n");
dt = -rsb_time();
mtxOp = mtxAp;
errval = rsb_tune_spmm(&mtxAp, &sf, &tn, oitmax, tmax, transA,
alphap, NULL, nrhs, order, NULL, ldB, betap, NULL, ldC);
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
dt += rsb_time();
if( mtxOp == mtxAp )
{
printf("After %lfs, global autotuning found old matrix optimal,"
" with declared speedup %lg x when using %d threads\n",dt,sf,tn);
}
else
{
printf("After %lfs, global autotuning declared speedup of %lg x,"
" when using threads count of %d and a new matrix:\n",dt,sf,tn);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s\n",ib);
}
printf("\n");
/* user is expected to:
errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
and use mtxAp in SpMV.
*/
rsb_mtx_free(mtxAp);
mtxAp = NULL;
}
rsb_mtx_free(mtxMp);
mtxMp = NULL;
err:
rsb_perror(NULL,errval);
if ( errval != RSB_ERR_NO_ERROR )
printf("Program terminating with error.\n");
return errval;
}
int main(const int argc, char * const argv[])
{
/*!
Autotuning example.
*/
/* matrix variables */
struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
const int bs = RSB_DEFAULT_BLOCKING;
rsb_coo_idx_t nrA = 500; /* number of rows */
rsb_coo_idx_t ncA = 500; /* number of cols */
const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
const rsb_coo_idx_t rd = 1;/* every rd rows one is non empty */
const rsb_coo_idx_t cd = 4;/* every cd cols one is non empty */
rsb_nnz_idx_t nnzA = (nrA/rd)*(ncA/cd); /* nonzeroes */
rsb_coo_idx_t*IA = NULL;
rsb_coo_idx_t*JA = NULL;
RSB_DEFAULT_TYPE*VA = NULL;
/* spmv specific variables */
const RSB_DEFAULT_TYPE alpha = 1;
const RSB_DEFAULT_TYPE beta = 1;
RSB_DEFAULT_TYPE*Cp = NULL;
RSB_DEFAULT_TYPE*Bp = NULL;
rsb_flags_t order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER;
const rsb_coo_idx_t nrhs = 2; /* number of right hand sides */
const rsb_trans_t transA = RSB_TRANSPOSITION_N;
rsb_nnz_idx_t ldB = nrA;
rsb_nnz_idx_t ldC = ncA;
/* misc variables */
rsb_err_t errval = RSB_ERR_NO_ERROR;
const size_t so = sizeof(RSB_DEFAULT_TYPE);
const size_t si = sizeof(rsb_coo_idx_t);
rsb_time_t dt,odt;
rsb_int_t t;
const rsb_int_t tt = 100; /* will repeat spmv tt times */
char ib[200];
const char*is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
/* misc counters */
rsb_coo_idx_t ci;
rsb_coo_idx_t ri;
rsb_coo_idx_t ni;
rsb_int_t nrhsi;
/* misc variables */
rsb_time_t etime = 0.0;
/* input autotuning variables */
const rsb_int_t oitmax = 15; /* auto-tune iterations */
const rsb_time_t tmax = 0.1; /* time per autotune operation */
/* input/output autotuning variables */
rsb_int_t tn = 0; /* threads number */
/* output autotuning variables */
rsb_real_t sf = 0.0; /* speedup factor obtained from auto tuning */
const rsb_int_t wvat = 1; /* want verbose autotuning; see
documentation of RSB_IO_WANT_VERBOSE_TUNING */
if(argc > 1 && !isdigit(argv[1][0]) )
{
errval = tune_from_file(argv[1],wvat);
goto ret;
}
if(argc > 1)
{
nrA = ncA = atoi(argv[1]);
if ( nrA < RSB_MIN_MATRIX_DIM || (nrA > (RSB_MAX_MATRIX_DIM) ))
goto err;
nnzA = (nrA/rd)*(ncA/cd);
ldB = nrA;
ldC = ncA;
}
printf("Creating %d x %d matrix with %d nonzeroes.\n",(int)nrA,
(int)ncA, (int)nnzA);
IA = calloc(nnzA, si);
JA = calloc(nnzA, si);
VA = calloc(nnzA, so);
Bp = calloc(nrhs*ncA ,so);
Cp = calloc(nrhs*nrA ,so);
if( ! ( VA && IA && JA && Bp && Cp ) )
goto err;
for(nrhsi=0;nrhsi<nrhs;++nrhsi)
for(ci=0;ci<ncA/cd;++ci)
Bp[nrhsi*ldC+ci] = 1.0;
for(nrhsi=0;nrhsi<nrhs;++nrhsi)
for(ri=0;ri<nrA/rd;++ri)
Cp[nrhsi*ldC+ri] = 1.0;
ni = 0;
for(ci=0;ci<ncA/cd;++ci)
for(ri=0;ri<nrA/rd;++ri)
{
VA[ni] = nrA * ri + ci,
IA[ni] = ri;
JA[ni] = ci;
ni++;
}
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS))
!= RSB_ERR_NO_ERROR) goto err;
errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
if( (errval) != RSB_ERR_NO_ERROR )
{
printf("Error setting option!\n");
goto err;
}
mtxAp = rsb_mtx_alloc_from_coo_const(
VA,IA,JA,nnzA,typecode,nrA,ncA,bs,bs,
RSB_FLAG_NOFLAGS,&errval);
/* VA, IA, JA are not necessary anymore */
free(VA);
free(IA);
free(JA);
VA = NULL;
IA = NULL;
JA = NULL;
if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
goto err;
printf("Allocated matrix of %zd nonzeroes:\n",(size_t)nnzA);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s\n\n",ib);
dt = - rsb_time();
for(t=0;t<tt;++t)
/*
If nrhs == 1, the following is equivalent to
rsb_spmv(transA,alphap,mtxAp,Bp,1,betap,Cp,1);
*/
rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
dt += rsb_time();
odt = dt;
printf("Before auto-tuning, %d multiplications took %lfs.\n",tt,dt);
printf("Threads autotuning (may take more than %lfs)...\n",
oitmax*tmax);
dt = -rsb_time();
errval = rsb_tune_spmm(NULL, &sf, &tn, oitmax, tmax, transA,
&alpha, mtxAp, nrhs, order, Bp, ldB, &beta, Cp, ldC);
dt += rsb_time();
if(errval != RSB_ERR_NO_ERROR)
goto err;
if(tn == 0)
printf("After %lfs, autotuning routine did not find a better"
" threads count configuration.\n",dt);
else
printf("After %lfs, autotuning routine declared speedup of %lg x,"
" when using threads count of %d.\n",dt,sf,tn);
errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
if(errval != RSB_ERR_NO_ERROR)
goto err;
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s\n",ib);
dt = -rsb_time();
for(t=0;t<tt;++t)
/*rsb_spmv(transA,&alpha,mtxAp,Bp,1,&beta,Cp,1);*/
rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
dt += rsb_time();
printf("After threads auto-tuning, %d multiplications took %lfs"
" -- effective speedup of %lg x\n",tt,dt,odt/dt);
odt = dt;
tn = 0; /* this will restore default threads count */
errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
if(errval != RSB_ERR_NO_ERROR)
goto err;
errval = rsb_lib_get_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
if(errval != RSB_ERR_NO_ERROR)
goto err;
printf("Matrix autotuning (may take more than %lfs; using %d"
" threads )...\n", oitmax*tmax, tn);
/* A negative tn will request also threads autotuning: */
/* tn = -tn; */
dt = -rsb_time();
errval = rsb_tune_spmm(&mtxAp, &sf, &tn, oitmax, tmax, transA,
&alpha, NULL, nrhs, order, Bp, ldB, &beta, Cp, ldC);
dt += rsb_time();
if(errval != RSB_ERR_NO_ERROR)
goto err;
if(tn == 0)
printf("After %lfs, autotuning routine did not find a better"
" threads count configuration.\n",dt);
else
printf("After %lfs, autotuning routine declared speedup of %lg x,"
" when using threads count of %d.\n",dt,sf,tn);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s\n",ib);
dt = -rsb_time();
for(t=0;t<tt;++t)
/*rsb_spmv(transA,&alpha,mtxAp,Bp,1,&beta,Cp,1);*/
rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
dt += rsb_time();
printf("After threads auto-tuning, %d multiplications took %lfs"
" -- further speedup of %lg x\n",tt,dt,odt/dt);
rsb_mtx_free(mtxAp);
free(Cp);
free(Bp);
errval = rsb_lib_get_opt(RSB_IO_WANT_LIBRSB_ETIME,&etime);
if(errval == RSB_ERR_UNSUPPORTED_FEATURE)
{
printf("librsb timer-based profiling is not supported in "
"this build. If you wish to have it, re-configure librsb "
"with its support. So you can safely ignore the error you"
" might just have seen printed out on screen.\n");
errval = RSB_ERR_NO_ERROR;
}
else
if(etime) /* This will only work if enabled at configure time. */
printf("Elapsed program time is %5.2lfs\n",etime);
ret:
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!=RSB_ERR_NO_ERROR)
goto err;
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.\n");
return EXIT_FAILURE;
}
examples/io-spblas.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
\ingroup rsb-examples
@file
@author Michele Martone
@brief Example C program using the Sparse BLAS interface
and reading from file using \ref rsb_blas_file_mtx_load(),
\ref BLAS_usgp(), \ref BLAS_dusmv(), \ref BLAS_usds().
\include io-spblas.c
*/
#include <rsb.h> /* for rsb_lib_init */
#include <blas_sparse.h>
#include <stdio.h>
int main(const int argc, char * const argv[])
{
#ifndef RSB_NUMERICAL_TYPE_DOUBLE
printf("Skipping a test because of 'double' type opted out.\n");
return EXIT_SUCCESS;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
blas_sparse_matrix A = blas_invalid_handle;
const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
const rsb_char_t * filename = argc > 1 ? argv[1] : "pd.mtx";
printf("Hello, RSB!\n");
if((rsb_perror(NULL,
rsb_lib_init(RSB_NULL_INIT_OPTIONS)))!=RSB_ERR_NO_ERROR)
{
printf("Error while initializing the library.\n");
goto err;
}
printf("Correctly initialized the library.\n");
A = rsb_blas_file_mtx_load(filename, typecode );
if( A == blas_invalid_handle )
{
printf("Error while loading matrix %s from file.\n",
filename);
goto err;
}
printf("Correctly loaded and allocated a matrix"
" from file %s.\n",filename);
if( BLAS_usgp(A,blas_symmetric) == 1 )
printf("Matrix is symmetric\n");
if( BLAS_usgp(A,blas_hermitian) == 1 )
printf("Matrix is hermitian\n");
printf("Now SPMV with NULL vectors will be attempted,"
" resulting in an error (so don't worry).\n");
if(BLAS_dusmv(blas_no_trans,-1,A,NULL,1,NULL,1))
{
printf("Correctly detected an error condition.\n");
goto okerr;
}
printf("No error detected ?\nIf you see this line printed out,"
" please report as a bug, because the above NULL pointers"
" should have been detected\n");
return EXIT_FAILURE;
okerr:
printf("Program correctly recovered from intentional"
" error condition.\n");
if(BLAS_usds(A))
{
printf("Error while freeing the matrix!\n");
goto err;
}
printf("Correctly freed the matrix.\n");
err:
if(rsb_perror(NULL,
rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))!=RSB_ERR_NO_ERROR)
{
printf("Failed finalizing the library.\n");
goto ferr;
}
printf("Correctly finalized the library.\n");
return EXIT_SUCCESS;
ferr:
return EXIT_FAILURE;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}
examples/transpose.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
@file
@author Michele Martone
@brief A toy <rsb.h>-based C program showing instantiation,
transposition and other operations on a single matrix.
Uses \ref rsb_mtx_clone(), \ref rsb_file_mtx_save(),
\ref rsb_file_mtx_get_dims(), \ref rsb_file_mtx_load().
\ingroup rsb-examples
\include transpose.c
*/
#include <rsb.h>
#include <stdio.h> /* printf */
int main(const int argc, char * const argv[])
{
struct rsb_mtx_t *mtxAp = NULL;
const rsb_blk_idx_t brA = RSB_DEFAULT_BLOCKING,
bcA = RSB_DEFAULT_BLOCKING;
rsb_nnz_idx_t nnzA = 4;
rsb_coo_idx_t nrA = 3;
rsb_coo_idx_t ncA = 3;
const rsb_coo_idx_t IA[] = { 0, 1, 2, 0 };
const rsb_coo_idx_t JA[] = { 0, 1, 2, 2 };
const RSB_DEFAULT_TYPE VA[] = { 11, 22, 33, 13 };
RSB_DEFAULT_TYPE XV[] = { 0,0,0,0,0,0 };
rsb_coo_idx_t vl = 0;
const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
rsb_err_t errval = RSB_ERR_NO_ERROR;
/* library initialization */
if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
{
return EXIT_FAILURE;
}
/* allocation */
mtxAp = rsb_mtx_alloc_from_coo_const(
VA,IA,JA,nnzA,typecode,nrA,ncA,
brA,bcA,RSB_FLAG_NOFLAGS,NULL);
if(!mtxAp)
{
return EXIT_FAILURE;
}
/* printout */
if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
{
if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
goto err;
}
/* matrix transposition */
if( RSB_ERR_NO_ERROR != (errval =
rsb_mtx_clone(&mtxAp,RSB_NUMERICAL_TYPE_SAME_TYPE,
RSB_TRANSPOSITION_T,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS)))
{
goto err;
}
/* printout */
if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
{
if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
goto err;
}
rsb_mtx_free(mtxAp);
/* doing the same after load from file */
mtxAp = rsb_file_mtx_load("pd.mtx",
RSB_FLAG_NOFLAGS,typecode,NULL);
if(!mtxAp)
{
return EXIT_FAILURE;
}
/* printout */
if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
{
if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
goto err;
}
/* one can see dimensions in advance, also */
if(RSB_ERR_NO_ERROR!=(errval =
rsb_file_mtx_get_dims("pd.mtx",&nrA,&ncA,&nnzA,NULL)))
{
if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
goto err;
}
/* A matrix can be rendered to Postscript. */
{
if(RSB_ERR_NO_ERROR!=(errval =
rsb_mtx_rndr("pd.eps",mtxAp,512,512,RSB_MARF_EPS_B)))
goto err;
}
rsb_mtx_free(mtxAp);
/* also vectors can be loaded */
if(RSB_ERR_NO_ERROR!=(errval =
rsb_file_vec_load("vf.mtx",typecode,NULL,&vl )))
goto err;
/* we expect vf.mtx to be 6 rows long */
if( vl != 6 )
{
goto err;
}
if(RSB_ERR_NO_ERROR!=(errval =
rsb_file_vec_load("vf.mtx",typecode,XV, NULL )))
goto err;
/* matrices can be rendered from file to a pixelmap as well */
{
unsigned char pixmap[3*2*2];
if(RSB_ERR_NO_ERROR!=(errval =
rsb_file_mtx_rndr(pixmap,"pd.mtx",2,2,2,RSB_MARF_RGB)))
goto err;
}
if(RSB_ERR_NO_ERROR != rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
{
goto err;
}
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
return EXIT_FAILURE;
}
examples/power.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
@file
@author Michele Martone
@brief A toy <rsb.h>-based C program implementing the power
method for computing matrix eigenvalues. Uses #rsb_spmv().
\ingroup rsb-examples
\include power.c
*/
#include <stdio.h> // printf
#include <math.h> // sqrt
#include <stdlib.h> // calloc
#include <rsb.h>
int main(const int argc, char * const argv[])
{
const int WANT_VERBOSE = 0;
struct rsb_mtx_t *mtxAp = NULL;
const int bs = RSB_DEFAULT_BLOCKING;
const int br = bs, bc = bs; /* bs x bs blocked */
const rsb_nnz_idx_t nnzA = 4;
const rsb_coo_idx_t nrA = 3;
const rsb_coo_idx_t ncA = 3;
rsb_int_t it = 0;
const rsb_int_t maxit = 100;
const rsb_coo_idx_t IA[] = { 0, 1, 2, 0 };
const rsb_coo_idx_t JA[] = { 0, 1, 2, 2 };
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[] = { 11, 22, 33, 13 };
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE ZERO = 0;
int i;
rsb_err_t errval = 0;
RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE norm = 0.0, /* nu */
oldnorm = 1.0, /* oldnorm */
*b1 = NULL, *b2 = NULL,
*bnow = NULL, *bnext = NULL;/* b1 and b2 aliases */
rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
size_t ds = 0;
/* tolerance */
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE tol = 1e-14;
/* library initialization */
if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
return EXIT_FAILURE;
/* allocation */
mtxAp = rsb_mtx_alloc_from_coo_const(VA,IA,JA,nnzA,
typecode,nrA,ncA,br,bc,RSB_FLAG_NOFLAGS,NULL);
if(!mtxAp)
return EXIT_FAILURE;
ds = (nrA)*sizeof(RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE);
b1 = calloc(1,ds);
b2 = calloc(1,ds);
if(! (b1 && b2))
{
errval = RSB_ERR_ENOMEM;
goto err;
}
for( i = 0; i < nrA; ++i )
b1[i] = 1;
bnow = b1, bnext = b2;/* b,b' */
while( fabs(norm-oldnorm) > tol && it<maxit )
{
++ it;
oldnorm = norm;
/* b'<-Ab */
if(( rsb_spmv(RSB_TRANSPOSITION_N,NULL,mtxAp,bnow,
1,&ZERO,bnext,1)) != RSB_ERR_NO_ERROR )
goto err;
/* nu<-||Ab||^2 */
norm = 0;
for(i=0;i<nrA;++i)
norm += bnext[i]*bnext[i];
/* nu<-||Ab|| */
norm = sqrt(norm);
norm = 1.0/norm;
/* b'<- Ab / ||Ab|| */
for(i=0;i<nrA;++i)
bnext[i] *= norm;
norm = 1.0/norm;
printf("it:%d norm:%lg norm diff:%lg\n",it,norm,norm-oldnorm);
{void *tmp=bnow;bnow=bnext;bnext=tmp;/* pointers swap */}
if(WANT_VERBOSE)
{
printf("norm:%lg\n",norm);
if(isinf(norm))
/* isinf is a C99 feature (need correct
* compilation flags) */
goto err;
for(i=0;i<2;++i)
printf("x[%d]=%lg\n",i,((double*)bnext)[i]);
}
}
/* the biggest eigenvalue should be in bnow */
rsb_mtx_free(mtxAp);
free(b1);
free(b2);
if(rsb_lib_exit(RSB_NULL_EXIT_OPTIONS)!=RSB_ERR_NO_ERROR)
goto err;
if( it == maxit )
{
printf("ERROR: hit iterations limit without convergence!");
errval=RSB_ERR_GENERIC_ERROR;
}
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
return EXIT_FAILURE;
}
examples/fortran.F90:
!
! Copyright (C) 2008-2022 Michele Martone
!
! This file is part of librsb.
!
! librsb is free software; you can redistribute it and/or modify it
! under the terms of the GNU Lesser General Public License as published
! by the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! librsb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with librsb; see the file COPYING.
! If not, see <http://www.gnu.org/licenses/>.
!
!> @file
!! @brief Sparse BLAS-based usage:
!! uscr_begin(), usgp(), ussp(), usmv(), ...
!! \include fortran.F90
SUBROUTINE blas_sparse_mod_example(res)
USE blas_sparse
USE rsb ! For the second part of the example and RSB_IDX_KIND
USE iso_c_binding
IMPLICIT NONE
INTEGER :: i, j
INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0, res
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr ! matrix pointer
INTEGER(C_INT) :: A
INTEGER,PARAMETER :: transn = blas_no_trans
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
REAL(KIND=8),parameter :: alpha = 3
! Symmetric (declared via lower triangle) matrix based example, e.g.:
! 1 0
! 1 1
! declaration of VA,IA,JA
!INTEGER,PARAMETER :: nr = 100
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nr = 20
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nc = nr
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nnz = (nr*(nr+1))/2 ! half the square
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 1
INTEGER(KIND=RSB_IDX_KIND) :: nt = 0
INTEGER :: ic, ir
! For compatibility with (or compassion for) flang-13, we turn these three lines off:
! REAL(KIND=8),PARAMETER :: VA(nnz) = (/ ((1, ic=1,ir), ir=1,nr ) /) ! (/1, 1, 1/)
! INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: IA(nnz) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
! INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: JA(nnz) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
! and use these other ones instead, initializing VA,IA,JA later:
REAL(KIND=8) :: va(nnz)
INTEGER(KIND=RSB_IDX_KIND) :: IA(nnz)
INTEGER(KIND=RSB_IDX_KIND) :: JA(nnz)
REAL(KIND=8) :: x(nc,nrhs) = reshape((/((1), ic=1,nc*nrhs)/),[nc,nrhs]) ! reference x ! (/1, 1/)
REAL(KIND=8),parameter :: cy(nr,nrhs) = reshape((/((alpha+alpha*nr), ir=1,nr*nrhs)/),[nr,nrhs]) ! reference cy after ! (/9, 9/)
REAL(KIND=8) :: y(nr,nrhs) = reshape((/((alpha), ir=1,nr*nrhs)/),[nr,nrhs]) ! y will be overwritten ! (/3, 3/)
va(:) = (/ ((1, ic=1,ir), ir=1,nr ) /) ! (/1, 1, 1/)
ia(:) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
ja(:) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
! First example part: pure blas_sparse code.
res = 0
CALL duscr_begin(nr,nc,a,res)
IF (res.NE.0) GOTO 9999
CALL ussp(a,blas_lower_symmetric,istat)
IF (istat.NE.0) GOTO 9997
CALL ussp(a,blas_rsb_spmv_autotuning_on,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) print *,"autotuning returned nonzero:", istat &
&," ...did you enable autotuning ?"
!
! First style example
CALL uscr_insert_entries(a,nnz,va,ia,ja,istat)
IF (istat.NE.0) GOTO 9997
CALL uscr_end(a,istat)
IF (istat.NE.0) GOTO 9997
! CALL ussp(A,blas_rsb_duplicates_sum,istat)
! CALL uscr_insert_entries(A,nnz,VA,IA,JA,istat) ! uncomment this to activate add of coefficients to pattern
CALL usgp(a,blas_rsb_spmv_autotuning_on,nt) ! (experimental)
IF (nt.NE.0) print*,"autotuner chose ",nt," threads"
CALL ussp(a,blas_rsb_spmv_autotuning_off,istat) ! (experimental) turns auto-tuning + thread setting off
IF (istat.NE.0) GOTO 9997
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
IF (istat.NE.0) GOTO 9997
!
DO j = 1, nrhs
DO i = 1, nr
IF (y(i,j).NE.cy(i,j)) print *, "first check results are not ok"
IF (y(i,j).NE.cy(i,j)) GOTO 9997
END DO
END DO
!
y(:,:) = alpha ! reset
!
! Second style example
CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) GOTO 9997
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat) ! Equivalent to the above (as long as incx=incy=1).
CALL usmm(blas_colmajor,transn,nrhs,-alpha,a,x,nr,y,nc,istat) ! Subtract the last usmm call contribution.
IF (istat.NE.0) GOTO 9997
!
DO j = 1, nrhs
DO i = 1, nr
IF (y(i,j).NE.cy(i,j)) print *,"second check results are not ok"
IF (y(i,j).NE.cy(i,j)) GOTO 9997
END DO
END DO
!
print *, "check results are ok"
! Second part of the example: access to the rsb.h interface via
! the ISO C Binding interface.
mtxap = rsb_blas_get_mtx(a) ! get pointer to rsb structure (as in the rsb.h API)
IF(nr.LT.5) istat = rsb_file_mtx_save(mtxap,c_null_char) ! write to stdout (only if matrix small enough)
GOTO 9998
9997 res = -1
9998 CONTINUE
CALL usds(a,istat)
IF (istat.NE.0) res = -1
9999 CONTINUE
end SUBROUTINE blas_sparse_mod_example
! Example loading a matrix from file and measuring SPMM.
SUBROUTINE blas_sparse_io_example(res)
USE blas_sparse
USE rsb ! For rsb_blas_file_mtx_load
USE iso_c_binding
IMPLICIT NONE
INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
INTEGER :: j
INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0
INTEGER :: A
INTEGER,PARAMETER :: transn = blas_no_trans
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
COMPLEX(KIND=8),PARAMETER :: alpha = 3
INTEGER(KIND=RSB_IDX_KIND) :: nr
INTEGER(KIND=RSB_IDX_KIND) :: nc
INTEGER(KIND=RSB_IDX_KIND) :: nz
INTEGER(KIND=RSB_IDX_KIND) :: st
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 4
COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: x
COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: y
CHARACTER(KIND=C_CHAR,LEN=7),TARGET :: filename = 'pd.mtx'//c_null_char
INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double_complex
REAL(KIND=c_double) :: mvt,mmt,omt
INTEGER(KIND=C_INT),TARGET::IZERO=0
res = 0
a = rsb_blas_file_mtx_load(filename, typecode);
IF (a.EQ.blas_invalid_handle) GOTO 9997
CALL usgp(a,blas_num_rows,nr)
CALL usgp(a,blas_num_cols,nc)
CALL usgp(a,blas_num_nonzeros,nz)
print*,"Read matrix ",filename(1:6)," ",nr,"x",nc,":",nz
CALL usgp(a,blas_general,st)
IF (st .EQ. 1) print*,"Matrix has no symmetry"
CALL usgp(a,blas_upper_symmetric,st)
IF (st .EQ. 1) print*,"Matrix is upper symmetric"
CALL usgp(a,blas_upper_hermitian,st)
IF (st .EQ. 1) print*,"Matrix is upper hermitian"
! ...
IF (istat.NE.0) GOTO 9997
WRITE(*,'(a,i0)') "Using NRHS=",nrhs
ALLOCATE( x(nc,nrhs))
ALLOCATE( y(nr,nrhs))
x = 1.0
y = 0.0
mvt = -rsb_time()
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
IF (istat.NE.0) GOTO 9997
mvt = mvt + rsb_time()
WRITE(*,'(a,e12.4,a)') "Repeated USMV took ",mvt," s"
y = 0.0
mmt = -rsb_time()
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
mmt = mmt + rsb_time()
WRITE(*,'(a,e12.4,a)') "A single USMM took ",mmt," s"
WRITE(*,'(a,g11.4,a)')"USMM-to-USMV speed ratio is is ", mvt/mmt, "x"
print*,"Call auto-tuning routine.."
! Change IZERO value to 1 to get verbose tuning again.
res = rsb_lib_set_opt(rsb_io_want_verbose_tuning,c_loc(izero))
IF (res.NE.0) GOTO 9997
CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) GOTO 9997
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
print*,"Repeat measurement."
y = 0.0
omt = -rsb_time()
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
omt = omt + rsb_time()
WRITE(*,'(a,e12.4,a)') "Tuned USMM took ",omt," s"
WRITE(*,'(a,g11.4,a)')"Tuned-to-untuned speed ratio is is ",mmt/omt,"x"
GOTO 9998
9997 res = -1
9998 CONTINUE
CALL usds(a,istat)
IF (istat.NE.0) res = -1
end SUBROUTINE blas_sparse_io_example
PROGRAM main
USE rsb, ONLY: rsb_lib_init, rsb_lib_exit, c_ptr, c_null_ptr,&
& rsb_io_want_extra_verbose_interface,rsb_io_want_verbose_tuning,&
& rsb_lib_set_opt,rsb_idx_kind
! USE blas_sparse, only: RSB_BLAS_IDX_KIND ! only if using long indices
USE iso_c_binding
IMPLICIT NONE
INTEGER :: passed = 0, failed = 0
INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
!TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
!TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
TYPE(C_PTR),PARAMETER :: EO = c_null_ptr
TYPE(C_PTR),PARAMETER :: IO = c_null_ptr
INTEGER(KIND=C_INT),TARGET::IONE=1
res = rsb_lib_init(io)
res = rsb_lib_set_opt(rsb_io_want_verbose_tuning,c_loc(ione))
CALL blas_sparse_mod_example(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
CALL blas_sparse_io_example(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
res = rsb_lib_exit(eo)
print *, "FAILED:", failed
print *, "PASSED:", passed
IF (failed .GT. 0) THEN
stop 1
END IF
END PROGRAM
examples/fortran_rsb_fi.F90:
!
! Copyright (C) 2008-2022 Michele Martone
!
! This file is part of librsb.
!
! librsb is free software; you can redistribute it and/or modify it
! under the terms of the GNU Lesser General Public License as published
! by the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! librsb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with librsb; see the file COPYING.
! If not, see <http://www.gnu.org/licenses/>.
!
!> @file.
!! @brief RSB.F90-based usage:
!! rsb_mtx_alloc_from_coo_const(),
!! rsb_tune_spmm(),
!! rsb_file_mtx_save(),
!! rsb_spmv(),
!! ...
!! \include fortran_rsb_fi.F90
SUBROUTINE rsb_mod_example1(res)
! Example using an unsymmetric matrix specified as COO, and autotuning, built at once.
USE rsb
USE iso_c_binding
IMPLICIT NONE
INTEGER ::res
INTEGER,TARGET :: istat = 0, i
INTEGER :: transt = rsb_transposition_n !
INTEGER(KIND=RSB_IDX_KIND) :: incx = 1, incy = 1
REAL(KIND=8),target :: alpha = 3, beta = 1
! 1 1
! 1 1
! declaration of VA,IA,JA
INTEGER(KIND=RSB_IDX_KIND) :: nnz = 4
INTEGER(KIND=RSB_IDX_KIND) :: nr = 2
INTEGER(KIND=RSB_IDX_KIND) :: nc = 2
INTEGER(KIND=RSB_IDX_KIND) :: nrhs = 1
INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
INTEGER :: flags = rsb_flag_noflags
INTEGER(KIND=RSB_IDX_KIND),TARGET :: IA(4) = (/0, 1, 1,0/)
INTEGER(KIND=RSB_IDX_KIND),TARGET :: JA(4) = (/0, 0, 1,1/)
REAL(KIND=8),target :: va(4) = (/1,1,1,1/)
REAL(KIND=8),target :: x(2) = (/1, 1/)! reference x
REAL(KIND=8),target :: cy(2) = (/9, 9/)! reference cy after
REAL(KIND=8),target :: y(2) = (/3, 3/)! y will be overwritten
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr ! matrix pointer
REAL(KIND=8) :: tmax = 2.0 ! tuning max time
INTEGER :: titmax = 2 ! tuning max iterations
INTEGER,TARGET :: ont = 0 ! optimal number of threads
!TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
!TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
TYPE(C_PTR),PARAMETER :: EO = c_null_ptr
TYPE(C_PTR),PARAMETER :: IO = c_null_ptr
INTEGER,TARGET :: errval
res = 0
errval = rsb_lib_init(io)
IF (errval.NE.rsb_err_no_error) GOTO 9997
mtxap = rsb_mtx_alloc_from_coo_const(c_loc(va),c_loc(ia),c_loc(ja)&
&,nnz,&
& rsb_numerical_type_double,nr,nc,1,1,flags,c_loc(istat))
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_file_mtx_save(mtxap,c_null_char)
! Structure autotuning:
istat = rsb_tune_spmm(c_loc(mtxap),c_null_ptr,c_null_ptr,titmax,&
& tmax,&
& transt,c_loc(alpha),c_null_ptr,nrhs,order,c_loc(x),nr,&
& c_loc(beta),c_loc(y),nc)
IF (istat.NE.rsb_err_no_error) GOTO 9997
! Thread count autotuning:
istat = rsb_tune_spmm(c_null_ptr,c_null_ptr,c_loc(ont),titmax,&
& tmax,&
& transt,c_loc(alpha),mtxap,nrhs,order,c_loc(x),nr,c_loc(beta),&
& c_loc(y),nc)
print *, "Optimal number of threads:", ont
y(:) = (/3, 3/)! restoring reference y (rsb_tune_spmm has overwritten it)
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_file_mtx_save(mtxap,c_null_char)
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_spmv(transt,c_loc(alpha),mtxap,c_loc(x),incx,&
& c_loc(beta),c_loc(y),incy)
IF (istat.NE.rsb_err_no_error) GOTO 9997
DO i = 1, 2
IF (y(i).NE.cy(i)) print *, "type=d dims=2x2 sym=g diag=g &
&blocks=1x1 usmv alpha= 3 beta= 1 incx=1 incy=1 trans=n is not ok"
IF (y(i).NE.cy(i)) GOTO 9997
END DO
print*,"type=d dims=2x2 sym=g diag=g blocks=1x1 usmv alpha= 3&
& beta= 1 incx=1 incy=1 trans=n is ok"
IF ( res .NE.rsb_err_no_error) GOTO 9997
GOTO 9998
9997 res = -1
9998 CONTINUE
mtxap = rsb_mtx_free(mtxap)
IF (istat.NE.rsb_err_no_error) res = -1
! 9999 CONTINUE
errval=rsb_lib_exit(eo) ! librsb finalization
IF (errval.NE.rsb_err_no_error)&
& stop "error calling rsb_lib_exit"
print *, "rsb module fortran test is ok"
istat = rsb_perror(c_null_ptr,istat)
end SUBROUTINE rsb_mod_example1
SUBROUTINE rsb_mod_example2(res)
! Example using an unsymmetric matrix specified as COO, and autotuning, built piecewise.
USE rsb
USE iso_c_binding
IMPLICIT NONE
INTEGER,TARGET :: errval
INTEGER :: res
INTEGER :: transt = rsb_transposition_n ! no transposition
INTEGER(KIND=RSB_IDX_KIND) :: incX = 1, incb = 1 ! X, B vectors increment
REAL(KIND=8),target :: alpha = 3,beta = 1
INTEGER(KIND=RSB_IDX_KIND) :: nnzA = 4, nra = 3, nca = 3 ! nonzeroes, rows, columns of matrix A
INTEGER(KIND=RSB_IDX_KIND),TARGET :: IA(4) = (/1, 2, 3, 3/) ! row indices
INTEGER(KIND=RSB_IDX_KIND),TARGET :: JA(4) = (/1, 2, 1, 3/) ! column indices
INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double
INTEGER :: flags =rsb_flag_default_matrix_flags+rsb_flag_symmetric
REAL(KIND=8),target :: va(4) = (/11.0, 22.0, 13.0, 33.0/) ! coefficients
REAL(KIND=8),target :: x(3) = (/ 0, 0, 0/)
REAL(KIND=8),target :: b(3) = (/-1.0, -2.0, -2.0/)
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr
TYPE(C_PTR) :: mtxApp = c_null_ptr
REAL(KIND=8),target :: etime = 0.0
!TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
!TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
TYPE(C_PTR),PARAMETER :: EO = c_null_ptr
TYPE(C_PTR),PARAMETER :: IO = c_null_ptr
errval = rsb_lib_init(io) ! librsb initialization
IF (errval.NE.rsb_err_no_error) &
& stop "error calling rsb_lib_init"
#if defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ < 5)
#define RSB_SKIP_BECAUSE_OLD_COMPILER 1
#endif
#ifndef RSB_SKIP_BECAUSE_OLD_COMPILER
mtxap = rsb_mtx_alloc_from_coo_begin(nnza,typecode,nra,nca,flags,&
& c_loc(errval)) ! begin matrix creation
errval = rsb_mtx_set_vals(mtxap,&
& c_loc(va),c_loc(ia),c_loc(ja),nnza,flags) ! insert some nonzeroes
mtxapp = c_loc(mtxap) ! Old compilers like e.g.: Gfortran 4.4.7 will NOT compile this.
IF (errval.NE.rsb_err_no_error) &
& stop "error calling rsb_mtx_set_vals"
errval = rsb_mtx_alloc_from_coo_end(mtxapp) ! end matrix creation
IF (errval.NE.rsb_err_no_error) &
& stop "error calling rsb_mtx_alloc_from_coo_end"
errval = rsb_spmv(transt,c_loc(alpha),mtxap,c_loc(x),&
& incx,c_loc(beta),c_loc(b),incb) ! X := X + (3) * A * B
IF (errval.NE.rsb_err_no_error)&
& stop "error calling rsb_spmv"
mtxap = rsb_mtx_free(mtxap) ! destroy matrix
! The following is optional and depends on configure options, so it is allowed to fail
errval = rsb_lib_get_opt(rsb_io_want_librsb_etime,c_loc(etime))
IF (errval.EQ.rsb_err_no_error)&
& print*,"Time spent in librsb is:",etime
! IF (errval.NE.0)STOP "error calling rsb_lib_get_opt"
errval = rsb_err_no_error
IF (errval.NE.rsb_err_no_error) &
& stop "error calling rsb_mtx_free"
#else
print*,"You have an old Fortran compiler not supporting C_LOC."
print*,"Skipping a part of the test"
#endif
errval=rsb_lib_exit(eo) ! librsb finalization
IF (errval.NE.rsb_err_no_error)&
& stop "error calling rsb_lib_exit"
print *, "rsb module fortran test is ok"
res = errval
end SUBROUTINE rsb_mod_example2
SUBROUTINE rsb_mod_example3(res)
! Example using a symmetric matrix specified as CSR, and autotuning, built at once.
USE rsb
USE iso_c_binding
IMPLICIT NONE
INTEGER ::res
INTEGER,TARGET :: istat = 0, i
INTEGER :: transt = rsb_transposition_n !
INTEGER(KIND=RSB_IDX_KIND) :: incx = 1, incy = 1
REAL(KIND=8),target :: alpha = 4, beta = 1
! 11 21
! 21 22
! declaration of VA,IP,JA
INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nnz = 3
INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nr = 2
INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nc = 2
INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nrhs = 1
INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
INTEGER :: flags = rsb_flag_noflags + rsb_flag_symmetric + &
& rsb_flag_fortran_indices_interface ! optional flags
INTEGER(KIND=RSB_IDX_KIND),TARGET :: IP(3) = (/1, 2, 4/)
INTEGER(KIND=RSB_IDX_KIND),TARGET :: JA(3) = (/1, 1, 2/)
REAL(KIND=8),target :: va(3) = (/11,21,22/) ! lower triangle coefficients
REAL(KIND=8),target :: x(2) = (/1, 2/)! reference x
REAL(KIND=8),target :: cy(2) = (/215.0, 264.0/)! reference cy after
REAL(KIND=8),target :: y(2) = (/3, 4/)! y will be overwritten
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr ! matrix pointer
REAL(KIND=8) :: tmax = 2.0 ! tuning max time
INTEGER :: titmax = 2 ! tuning max iterations
INTEGER,TARGET :: ont = 0 ! optimal number of threads
TYPE(C_PTR),PARAMETER :: EO = c_null_ptr
TYPE(C_PTR),PARAMETER :: IO = c_null_ptr
INTEGER,TARGET :: errval
errval = rsb_lib_init(io) ! librsb initialization
IF (errval.NE.rsb_err_no_error) &
& stop "error calling rsb_lib_init"
res = 0
mtxap = rsb_mtx_alloc_from_csr_const(c_loc(va),c_loc(ip),c_loc(ja)&
&,nnz,rsb_numerical_type_double,nr,nc,1,1,flags,c_loc(istat))
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_file_mtx_save(mtxap,c_null_char)
! Structure autotuning:
istat = rsb_tune_spmm(c_loc(mtxap),c_null_ptr,c_null_ptr,titmax,&
& tmax,&
& transt,c_loc(alpha),c_null_ptr,nrhs,order,c_loc(x),nr,&
& c_loc(beta),c_loc(y),nc)
IF (istat.NE.rsb_err_no_error) GOTO 9997
! Thread count autotuning:
istat = rsb_tune_spmm(c_null_ptr,c_null_ptr,c_loc(ont),titmax,&
& tmax,&
& transt,c_loc(alpha),mtxap,nrhs,order,c_loc(x),nr,c_loc(beta),&
& c_loc(y),nc)
print *, "Optimal number of threads:", ont
y(:) = (/3, 4/)! restoring reference y (rsb_tune_spmm has overwritten it)
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_file_mtx_save(mtxap,c_null_char)
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_spmv(transt,c_loc(alpha),mtxap,c_loc(x),incx,&
& c_loc(beta),c_loc(y),incy)
print *, y
IF (istat.NE.rsb_err_no_error) GOTO 9997
DO i = 1, 2
IF (y(i).NE.cy(i)) print *, "type=d dims=2x2 sym=s diag=g &
&blocks=1x1 usmv alpha= 4 beta= 1 incx=1 incy=1 trans=n is not ok"
IF (y(i).NE.cy(i)) GOTO 9997
END DO
print*,"type=d dims=2x2 sym=s diag=g blocks=1x1 usmv alpha= 4&
& beta= 1 incx=1 incy=1 trans=n is ok"
GOTO 9998
errval=rsb_lib_exit(eo) ! librsb finalization
IF (errval.NE.rsb_err_no_error)&
& stop "error calling rsb_lib_exit"
print *, "rsb module fortran test is ok"
9997 res = -1
9998 CONTINUE
mtxap = rsb_mtx_free(mtxap)
IF (istat.NE.rsb_err_no_error) res = -1
! 9999 CONTINUE
istat = rsb_perror(c_null_ptr,istat)
end SUBROUTINE rsb_mod_example3
PROGRAM main
USE rsb
IMPLICIT NONE
INTEGER :: res = rsb_err_no_error, passed = 0, failed = 0
CALL rsb_mod_example1(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
CALL rsb_mod_example2(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
CALL rsb_mod_example3(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
print *, "FAILED:", failed
print *, "PASSED:", passed
IF (failed.GT.0) THEN
stop 1
END IF
END PROGRAM
examples/backsolve.c:
/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.
*/
/*!
\ingroup rsb-examples
@file
@author Michele Martone
@brief C triangular solve example program.
Uses #rsb_spsv(),#rsb_tune_spsm().
Based on <rsb.h>.
\include hello.c
*/
#include <rsb.h> /* librsb header to include */
#include <stdio.h> /* printf() */
int main(const int argc, char * const argv[])
{
/*!
A Hello-RSB program.
This program shows how to use the rsb.h interface correctly to:
- initialize the library using #rsb_lib_init()
- allocate (build) a single sparse matrix in the RSB format
using #rsb_mtx_alloc_from_coo_const(), with implicit diagonal
- print information obtained via #rsb_mtx_get_info_str()
- multiply the triangular matrix using #rsb_spmv()
- autotune the matrix for rsb_spsv with #rsb_tune_spsm()
- solve the triangular system using #rsb_spsv()
- deallocate the matrix using #rsb_mtx_free()
- finalize the library using #rsb_lib_exit()
In this example, we use #RSB_DEFAULT_TYPE as matrix type.
This type depends on what was configured at library build time.
* */
const int bs = RSB_DEFAULT_BLOCKING;
const int brA = bs, bcA = bs;
const RSB_DEFAULT_TYPE one = 1;
const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
const rsb_nnz_idx_t nnzA = 7; /* matrix nonzeroes count */
const rsb_coo_idx_t nrA = 6; /* matrix rows count */
const rsb_coo_idx_t ncA = 6; /* matrix columns count */
/* nonzero row indices coordinates: */
const rsb_coo_idx_t IA[] = {1,2,3,4,5,6,1};
/* nonzero column indices coordinates: */
const rsb_coo_idx_t JA[] = {1,2,3,4,5,6,6};
const RSB_DEFAULT_TYPE VA[] = {1,1,1,1,1,1,1};/*values of nonzeroes*/
RSB_DEFAULT_TYPE X[] = { 0,0,0,0,0,0 }; /* X vector's array */
const RSB_DEFAULT_TYPE B[] = { 1,1,1,1,1,1 }; /* B */
struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
char ib[200];
int i;
rsb_err_t errval = RSB_ERR_NO_ERROR;
const rsb_int_t wvat = 1; /* want verbose autotuning; see
documentation of RSB_IO_WANT_VERBOSE_TUNING */
printf("Hello, RSB!\n");
printf("Initializing the library...\n");
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) !=
RSB_ERR_NO_ERROR)
{
printf("Error initializing the library!\n");
goto err;
}
printf("Correctly initialized the library.\n");
errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
if( (errval) != RSB_ERR_NO_ERROR )
{
printf("Error setting option!\n");
goto err;
}
mtxAp = rsb_mtx_alloc_from_coo_const(
VA,IA,JA,nnzA,typecode,nrA,ncA,brA,bcA,
RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS /* force rsb */
| RSB_FLAG_DUPLICATES_SUM/* sum dups */
| RSB_FLAG_UNIT_DIAG_IMPLICIT/* ask diagonal implicit */
| RSB_FLAG_TRIANGULAR /* need triangle for spsv */
| RSB_FLAG_FORTRAN_INDICES_INTERFACE /* treat indices as 1-based */
, &errval);
if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
{
printf("Error while allocating the matrix!\n");
goto err;
}
printf("Correctly allocated a matrix with %ld nonzeroes.\n",
(long int)nnzA);
printf("Summary information of the matrix:\n");
/* print out the matrix summary information */
rsb_mtx_get_info_str(mtxAp,"RSB_MIF_MATRIX_INFO__TO__CHAR_P",
ib,sizeof(ib));
printf("%s",ib);
printf("\nMatrix printout:\n");
rsb_file_mtx_save(mtxAp, NULL);
if((errval =
rsb_spmv(RSB_TRANSPOSITION_N,&one,mtxAp,B,1,&one,X,1))
!= RSB_ERR_NO_ERROR )
{
printf("Error performing a multiplication!\n");
goto err;
}
printf("\nWe have a unitary vector:\n");
rsb_file_vec_save(NULL, typecode, B, nrA);
printf("\nMultiplying matrix by unitary vector we get:\n");
rsb_file_vec_save(NULL, typecode, X, nrA);
errval = rsb_tune_spsm(&mtxAp, NULL, NULL, 0, 0, RSB_TRANSPOSITION_N,
&one, NULL, 1, RSB_FLAG_WANT_COLUMN_MAJOR_ORDER, NULL, nrA,
NULL, NULL, nrA);
if( (errval) != RSB_ERR_NO_ERROR )
{
printf("Error performing autotuning!\n");
goto err;
}
if((errval = rsb_spsv(RSB_TRANSPOSITION_N,&one,mtxAp,X,1,X,1))
!= RSB_ERR_NO_ERROR )
{
printf("Error performing triangular solve!\n");
goto err;
}
printf("\nBacksolving we should get a unitary vector:\n");
rsb_file_vec_save(NULL, typecode, X, nrA);
for(i=0;i<nrA;++i)
if(X[i]!=one)
{
printf("Warning! Result vector not unitary!:\n");
errval = RSB_ERR_INTERNAL_ERROR;
goto err;
}
printf("All done.\n");
rsb_mtx_free(mtxAp);
printf("Correctly freed the matrix.\n");
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
printf("Error finalizing the library!\n");
goto err;
}
printf("Correctly finalized the library.\n");
printf("Program terminating with no error.\n");
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.\n");
return EXIT_FAILURE;
}
examples/bench.sh:
#!/bin/bash
#
# Copyright (C) 2008-2022 Michele Martone
#
# This file is part of librsb.
#
# librsb is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# librsb is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
# License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with librsb; see the file COPYING.
# If not, see <http://www.gnu.org/licenses/>.
# Benchmark rsbench and postprocess results.
# \ingroup rsb-examples
# @file
# @author Michele Martone
# @brief Benchmark invocation from shell script.
#
set -e
set -x
which rsbench
BRF=test.rpr
# invoke rsbench and produce a performance record, using all types and one thread
rsbench -oa -Ob --bench --lower 100 --as-symmetric --types ':' -n 1 --notranspose --compare-competitors --verbose --verbose --write-performance-record=${BRF}
# examine tuning renderings (produced by --verbose --verbose)
ls -ltr ${BRF/.rpr/}-tuning*
# convert the performance record to text form
rsbench --read-performance-record ${BRF} > ${BRF/rpr/txt}
ls -ltr ${BRF/rpr/txt}
# convert the performance record to LaTeX table document form
RSB_PR_WLTC=2 RSB_PR_SR=0 rsbench --read-performance-record ${BRF} > ${BRF/rpr/tex}
which latex || exit 0 # to compile LaTeX document
which kpsepath || exit 0 # to check LaTeX packages
find `kpsepath tex | sed 's/!!//g;s/:/\n/g;'` -name sciposter.cls || exit 0 # need sciposter class, usually in texlive-science
find `kpsepath tex | sed 's/!!//g;s/:/\n/g;'` -name translator.sty || exit 0 # need sciposter class, usually in texlive-latex-recommended
# convert the LaTeX table into a DVI (may as well use pdflatex for PDF)
latex -interaction=batchmode -file-line-error ${BRF/rpr/tex}
# convert the performance record to GNUPLOT plots
which gnuplot || exit 0
RSB_PRD_STYLE_PLT_PFN=${BRF/rpr/} RSB_PRD_STYLE_PLT_FMT=1 RSB_PR_SR=2 rsbench --read-performance-record ${BRF} > ${BRF/rpr/gnu}
# convert the GNUPLOT plots into PDF
ls -ltr ${BRF/rpr/gnu}
gnuplot ${BRF/rpr/gnu}
# convert the performance record to GNUPLOT plots, different way
RSB_PRD_STYLE_PLT_PFN=${BRF/rpr/} RSB_PRD_STYLE_PLT_FMT= RSB_PR_SR=2 rsbench --read-performance-record ${BRF} > ${BRF/rpr/gnu}
gnuplot ${BRF/rpr/gnu}
ls -ltr ${BRF/.rpr/}*.png
ls -ltr ${BRF/.rpr/}*.eps
ls -ltr ${BRF} ${BRF/rpr/tex} ${BRF/rpr/dvi} ${BRF/rpr/txt}
# clean up
rm ${BRF/rpr/}{aux,log,out,gnu}
rm ${BRF/.rpr/}*{.png,.eps}
rm ${BRF} ${BRF/rpr/tex} ${BRF/rpr/dvi} ${BRF/rpr/txt}
exit
Most of the snippets in the documentation come from examples/snippets.c.
Author
librsb was written by Michele Martone; this documentation has been generated by Doxygen.
SEE ALSO
rsb-examples rsb-spblas.h rsb.h rsb.hpp