librsb  1.2.0.10
Example programs and code

Examples of usage of librsb.

Examples of usage of librsb.

    The following fully working example programs illustrate correct ways of using the library.
   The script displayed here should be sufficient to build them.
#!/bin/bash
# Script to build the librsb example programs.
LIBRSB_CONFIG=${LIBRSB_CONFIG:-librsb-config}
for s in *.c
do
p=${s/.c/}
rm -f $p
CFLAGS=`${LIBRSB_CONFIG} --I_opts`
LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
CC=`${LIBRSB_CONFIG} --cc`
cmd="$CC $CFLAGS $s $LDFLAGS -o $p"
echo $cmd
$cmd
done
if test x"yes" = x"yes" ; then
# activated if you have built the Fortran modules and installed them in the right path.
for s in *.F90
do
p=${s/.F90/}
rm -f $p
CFLAGS=`${LIBRSB_CONFIG} --I_opts`
LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
FC=`${LIBRSB_CONFIG} --fc`
cmd="$FC $CFLAGS $s $LDFLAGS -o $p"
echo $cmd
$cmd
done
fi
/*
Copyright (C) 2008-2020 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_doc_examples
@file
@author Michele Martone
@brief This is a first "hello RSB" example program.
\include hello.c
*/
#include <rsb.h> /* librsb header to include */
#include <stdio.h> /* printf() */
#include <stdlib.h> // EXIT_SUCCESS
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(RSB_NULL_EXIT_OPTIONS)
In this example, we use #RSB_DEFAULT_TYPE as matrix type.
This type depends on what was configured at library build time.
* */
struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
const int bs = RSB_DEFAULT_BLOCKING;
const int brA = bs, bcA = bs;
const RSB_DEFAULT_TYPE one = 1;
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: */
rsb_coo_idx_t IA[] = {0,1,2,2};
/* nonzero column indices coordinates: */
rsb_coo_idx_t JA[] = {0,1,2,2};
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];
printf("Hello, RSB!\n");
printf("Initializing the library...\n");
{
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(
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);
{
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;
&evi);
errval = RSB_ERR_NO_ERROR;
}
}
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))
{
printf("Error performing a multiplication!\n");
goto err;
}
printf("Correctly performed a SPMV.\n");
rsb_mtx_free(mtxAp);
printf("Correctly freed the matrix.\n");
{
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;
}
rsb_err_t rsb_strerror_r(rsb_err_t errval, rsb_char_t *buf, size_t buflen)
Definition: rsb_rsb.c:713
rsb_err_t rsb_mtx_get_info_str(const struct rsb_mtx_t *mtxAp, const rsb_char_t *mis, void *minfop, size_t buflen)
Definition: rsb_rsb.c:1425
rsb_err_t rsb_spmv(rsb_trans_t transA, const void *alphap, const struct rsb_mtx_t *mtxAp, const void *Xp, rsb_coo_idx_t incX, const void *betap, void *Yp, rsb_coo_idx_t incY)
Definition: rsb_rsb.c:467
struct rsb_mtx_t * rsb_mtx_alloc_from_coo_const(const void *VA, const rsb_coo_idx_t *IA, const rsb_coo_idx_t *JA, rsb_nnz_idx_t nnzA, rsb_type_t typecode, rsb_coo_idx_t nrA, rsb_coo_idx_t ncA, rsb_blk_idx_t brA, rsb_blk_idx_t bcA, rsb_flags_t flagsA, rsb_err_t *errvalp)
Definition: rsb_rsb.c:203
rsb_err_t rsb_lib_exit(struct rsb_initopts *iop)
Definition: rsb_rsb.c:173
struct rsb_mtx_t * rsb_mtx_free(struct rsb_mtx_t *mtxAp)
Definition: rsb_rsb.c:254
rsb_err_t rsb_perror(void *stream, rsb_err_t errval)
Definition: rsb_rsb.c:691
rsb_err_t rsb_lib_init(struct rsb_initopts *iop)
Definition: rsb_rsb.c:52
@ RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE
Definition: rsb.h:658
This file declares the user interface functions and data structures for the librsb library.
#define RSB_FLAG_NOFLAGS
Definition: rsb.h:396
#define RSB_ERR_NO_ERROR
Definition: rsb.h:534
#define RSB_NULL_EXIT_OPTIONS
Definition: rsb.h:763
rsb_err_t rsb_lib_set_opt(enum rsb_opt_t iof, const void *iop)
Definition: rsb_rsb.c:91
#define RSB_DEFAULT_BLOCKING
Definition: rsb.h:347
signed int rsb_nnz_idx_t
Definition: rsb.h:281
signed int rsb_coo_idx_t
Definition: rsb.h:272
#define RSB_ERRS_UNSUPPORTED_FEATURES
Definition: rsb.h:591
signed int rsb_err_t
Definition: rsb.h:304
char rsb_type_t
Definition: rsb.h:296
#define RSB_NULL_INIT_OPTIONS
Definition: rsb.h:762
signed int rsb_int_t
Definition: rsb.h:310
#define RSB_FLAG_DUPLICATES_SUM
Definition: rsb.h:444
#define RSB_NUMERICAL_TYPE_DEFAULT
A default numerical matrix type.
Definition: rsb_types.h:201
#define RSB_TRANSPOSITION_N
N: Non transposed flag, valid for rsb_trans_t typed variables.
Definition: rsb_types.h:145
#define RSB_DEFAULT_TYPE
The default numerical matrix type (can be used for declarations), used in the example programs.
Definition: rsb_types.h:80
/*
Copyright (C) 2008-2020 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_doc_examples
@file
@author Michele Martone
@brief This is a first "hello RSB" example program using
a Sparse BLAS interface.
\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 */
#include <stdlib.h> // EXIT_SUCCESS
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): */
int IA[] = { 0, 1, 2, 2 };
/* A's nonzero elements column indices (coordinates): */
int JA[] = { 0, 1, 0, 2 };
/* 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: */
double B[] = { -1.0, -2.0, -2.0 };
/* the (known) result array: */
double AB[] = { 11.0+26.0, 44.0, 66.0+13.0 };
/* rsb error variable: */
int i;
printf("Hello, RSB!\n");
/* initialize the library */
{
goto err;
}
printf("Correctly initialized the library.\n");
/* initialize a matrix descriptor */
A = BLAS_duscr_begin(nr,nc);
{
goto err;
}
/* specify properties (e.g.: symmetry)*/
{
goto err;
}
/* get properties (e.g.: symmetry) */
{
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)
{
goto err;
}
/* finalize (allocate) the matrix build */
{
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 */
{
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 */
}
This file specifies the Sparse BLAS interface to librsb. Supported types :(float,double,...
@ blas_no_trans
Definition: blas_sparse.h:32
#define BLAS_usgp
Definition: blas_sparse.h:561
#define BLAS_ussp
Definition: blas_sparse.h:560
int blas_sparse_matrix
Definition: blas_sparse.h:165
@ blas_invalid_handle
Definition: blas_sparse.h:133
@ blas_lower_symmetric
Definition: blas_sparse.h:113
int BLAS_duscr_end(blas_sparse_matrix A)
Definition: rsb_libspblas.c:1877
int BLAS_usds(blas_sparse_matrix A)
Definition: rsb_libspblas.c:2642
int BLAS_dusmv(enum blas_trans_type transA, double alpha, blas_sparse_matrix A, const double *x, int incx, double *y, int incy)
Definition: rsb_libspblas.c:934
blas_sparse_matrix BLAS_duscr_begin(int m, int n)
Definition: rsb_libspblas.c:1470
int BLAS_dusget_element(blas_sparse_matrix A, int i, int j, double *v)
Definition: rsb_libspblas.c:3535
int BLAS_duscr_insert_entries(blas_sparse_matrix A, int nnz, const double *val, const int *indx, const int *jndx)
Definition: rsb_libspblas.c:2088
/*
Copyright (C) 2008-2020 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_doc_examples
@file
@author Michele Martone
@brief This is a first "hello RSB" example program using
a Sparse BLAS interface.
\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 */
#include <stdlib.h> // EXIT_SUCCESS
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): */
int IA[] = { 0, 1, 2, 2 };
/* A's nonzero elements column indices (coordinates): */
int JA[] = { 0, 1, 0, 2 };
/* 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: */
double B[] = { -1.0, -2.0, -2.0 };
/* the (known) result array: */
double AB[] = { 11.0+26.0, 44.0, 66.0+13.0 };
/* rsb error variable: */
int i;
printf("Hello, RSB!\n");
/* initialize the library */
{
goto err;
}
printf("Correctly initialized the library.\n");
/* initialize a matrix descriptor */
A = BLAS_duscr_begin(nr,nc);
{
goto err;
}
/* specify properties (e.g.: symmetry)*/
{
goto err;
}
/* get properties (e.g.: symmetry) */
{
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)
{
goto err;
}
/* finalize (allocate) the matrix build */
{
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 */
{
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 */
}
/*
Copyright (C) 2008-2020 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_doc_examples
@file
@author Michele Martone
@brief This is a first "RSB autotuning" example program.
\include autotuning.c
*/
#include <rsb.h> /* librsb header to include */
#include <stdio.h> /* printf() */
#include <ctype.h> /* isdigit() */
#include <stdlib.h> /* atoi() */
/* #include "rsb_internals.h" */
int tune_from_file(char * const filename, rsb_int_t wvat)
{
struct rsb_mtx_t *mtxMp = NULL;
/* spmv specific variables */
const RSB_DEFAULT_TYPE alpha = 1;
const RSB_DEFAULT_TYPE beta = 1;
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 */
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 */
int ione = 1;
rsb_type_t typecodea [] = RSB_MATRIX_SPBLAS_TYPE_CODES_ARRAY;
int typecodei;
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
/*
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,
&alpha, mtxAp, nrhs, order, NULL, ldB, &beta, 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,
&alpha, NULL, nrhs, order, NULL, ldB, &beta, 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 */
rsb_coo_idx_t rd = 1; /* every rd rows one is non empty */
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;
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 = nrA;
rsb_nnz_idx_t ldC = ncA;
/* misc variables */
size_t so = sizeof(RSB_DEFAULT_TYPE);
size_t si = sizeof(rsb_coo_idx_t);
rsb_time_t dt,odt;
rsb_int_t t,tt = 100; /* will repeat spmv tt times */
char ib[200];
const char*is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
/* misc counters */
rsb_int_t nrhsi;
/* misc variables */
rsb_time_t etime = 0.0;
/* input autotuning variables */
rsb_int_t oitmax = 15; /* auto-tune iterations */
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 */
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",nrA,ncA,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++;
}
!= RSB_ERR_NO_ERROR) goto err;
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,&alpha,mtxAp,Bp,1,&beta,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);
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 */
if(errval != RSB_ERR_NO_ERROR)
goto err;
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);
{
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:
goto err;
return EXIT_SUCCESS;;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.\n");
return EXIT_FAILURE;
}
struct rsb_mtx_t * rsb_file_mtx_load(const rsb_char_t *filename, rsb_flags_t flagsA, rsb_type_t typecode, rsb_err_t *errvalp)
Definition: rsb_rsb.c:879
rsb_err_t rsb_spmm(rsb_trans_t transA, const void *alphap, const struct rsb_mtx_t *mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void *Bp, rsb_nnz_idx_t ldB, const void *betap, void *Cp, rsb_nnz_idx_t ldC)
Definition: rsb_rsb.c:1240
rsb_time_t rsb_time(void)
Definition: rsb_rsb.c:1489
rsb_err_t rsb_tune_spmm(struct rsb_mtx_t **mtxOpp, rsb_real_t *sfp, rsb_int_t *tnp, rsb_int_t maxr, rsb_time_t maxt, rsb_trans_t transA, const void *alphap, const struct rsb_mtx_t *mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void *Bp, rsb_nnz_idx_t ldB, const void *betap, void *Cp, rsb_nnz_idx_t ldC)
Definition: rsb_rsb.c:1567
rsb_err_t rsb_mtx_clone(struct rsb_mtx_t **mtxBpp, rsb_type_t typecode, rsb_trans_t transA, const void *alphap, const struct rsb_mtx_t *mtxAp, rsb_flags_t flags)
Definition: rsb_rsb.c:274
@ RSB_IO_WANT_VERBOSE_TUNING
Definition: rsb.h:694
@ RSB_IO_WANT_EXECUTING_THREADS
Definition: rsb.h:654
@ RSB_IO_WANT_LIBRSB_ETIME
Definition: rsb.h:690
#define RSB_MAX_MATRIX_DIM
Definition: rsb.h:360
rsb_flags_t rsb_trans_t
Definition: rsb.h:318
signed int rsb_flags_t
Definition: rsb.h:291
#define RSB_FLAG_WANT_COLUMN_MAJOR_ORDER
Definition: rsb.h:414
rsb_err_t rsb_lib_get_opt(enum rsb_opt_t iof, void *iop)
Definition: rsb_rsb.c:108
double rsb_real_t
Definition: rsb.h:321
rsb_real_t rsb_time_t
Definition: rsb.h:329
#define RSB_MIN_MATRIX_DIM
Definition: rsb.h:357
#define RSB_ERR_UNSUPPORTED_FEATURE
Definition: rsb.h:567
/*
Copyright (C) 2008-2020 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_doc_examples
@file
@author Michele Martone
@brief This is an example program using a Sparse BLAS interface
and reading from file using the RSB library.
\include io-spblas.c
*/
#include <rsb.h> /* for rsb_lib_init */
#include <blas_sparse.h>
#include <stdio.h>
#include <stdlib.h> // EXIT_SUCCESS
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 */
rsb_char_t * filename = argc > 1 ? argv[1] : "pd.mtx";
printf("Hello, RSB!\n");
if((rsb_perror(NULL,
{
printf("Error while initializing the library.\n");
goto err;
}
printf("Correctly initialized the library.\n");
typecode );
{
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,
{
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 */
}
blas_sparse_matrix rsb_load_spblas_matrix_file_as_matrix_market(const rsb_char_t *filename, rsb_type_t typecode)
Definition: rsb_libspblas_handle.c:1716
@ blas_hermitian
Definition: blas_sparse.h:109
@ blas_symmetric
Definition: blas_sparse.h:108
char rsb_char_t
Definition: rsb.h:326
#define RSB_NUMERICAL_TYPE_DOUBLE
Character code for type double.
Definition: rsb_types.h:186
/*
Copyright (C) 2008-2020 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 program showing instantiation, transposition and other
operations on a single matrix.
\ingroup rsb_doc_examples
\include transpose.c
*/
#include <rsb.h>
#include <stdio.h> /* printf */
#include <stdlib.h> // EXIT_SUCCESS
int main(const int argc, char * const argv[])
{
struct rsb_mtx_t *mtxAp = NULL;
rsb_nnz_idx_t nnzA = 4;
rsb_coo_idx_t nrA = 3;
rsb_coo_idx_t ncA = 3;
rsb_coo_idx_t IA[] = { 0, 1, 2, 0 };
rsb_coo_idx_t JA[] = { 0, 1, 2, 2 };
RSB_DEFAULT_TYPE VA[] = { 11, 22, 33, 13 };
RSB_DEFAULT_TYPE XV[] = { 0,0,0,0,0,0 };
rsb_coo_idx_t vl = 0;
/* library initialization */
{
return EXIT_FAILURE;
}
/* allocation */
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)))
{
goto err;
}
/* matrix transposition */
if( RSB_ERR_NO_ERROR != (errval =
{
goto err;
}
/* printout */
if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
{
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)))
{
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)))
{
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;
}
{
goto err;
}
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
return EXIT_FAILURE;
}
rsb_err_t rsb_file_mtx_rndr(void *pmp, const rsb_char_t *filename, rsb_coo_idx_t pmlWidth, rsb_coo_idx_t pmWidth, rsb_coo_idx_t pmHeight, rsb_marf_t rflags)
Definition: rsb_rsb.c:1328
rsb_err_t rsb_file_mtx_get_dims(const rsb_char_t *filename, rsb_coo_idx_t *nrp, rsb_coo_idx_t *ncp, rsb_coo_idx_t *nzp, rsb_flags_t *flagsp)
Definition: rsb_rsb.c:668
rsb_err_t rsb_mtx_rndr(const rsb_char_t *filename, const struct rsb_mtx_t *mtxAp, rsb_coo_idx_t pmWidth, rsb_coo_idx_t pmHeight, rsb_marf_t rflags)
Definition: rsb_rsb.c:1307
rsb_err_t rsb_file_vec_load(const rsb_char_t *filename, rsb_type_t typecode, void *Yp, rsb_coo_idx_t *yvlp)
Definition: rsb_rsb.c:856
rsb_err_t rsb_file_mtx_save(const struct rsb_mtx_t *mtxAp, const rsb_char_t *filename)
Definition: rsb_rsb.c:810
#define RSB_MARF_RGB
Definition: rsb.h:805
#define RSB_MARF_EPS_B
Definition: rsb.h:807
#define RSB_FLAG_IDENTICAL_FLAGS
Definition: rsb.h:399
signed int rsb_blk_idx_t
Definition: rsb.h:264
#define RSB_TRANSPOSITION_T
T: Transposed flag value, valid for rsb_trans_t valued variables.
Definition: rsb_types.h:146
#define RSB_NUMERICAL_TYPE_SAME_TYPE
a bogus type flag for specifying no type conversion
Definition: rsb_types.h:191
/*
Copyright (C) 2008-2020 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 program implementing the power method
for computing matrix eigenvalues.
\ingroup rsb_doc_examples
\include power.c
*/
#include <stdio.h> // printf
#include <math.h> // sqrt
#include <stdlib.h> // calloc, EXIT_SUCCESS
#include <rsb.h>
int main(const int argc, char * const argv[])
{
int WANT_VERBOSE = 0;
struct rsb_mtx_t *mtxAp = NULL;
const int bs = RSB_DEFAULT_BLOCKING;
int i;
const int br = bs, bc = bs; /* bs x bs blocked */
rsb_err_t errval = 0;
rsb_nnz_idx_t nnzA = 4;
rsb_coo_idx_t nrA = 3;
rsb_coo_idx_t ncA = 3;
rsb_int_t it = 0, 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 };
oldnorm = 1.0, /* oldnorm */
*b1 = NULL, *b2 = NULL,
*bnow = NULL, *bnext = NULL;/* b1 and b2 aliases */
size_t ds = 0;
/* tolerance */
/* library initialization */
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;
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);
goto err;
if( it == maxit )
{
printf("ERROR: hit iterations limit without convergence!");
}
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
return EXIT_FAILURE;
}
#define RSB_ERR_ENOMEM
Definition: rsb.h:555
#define RSB_ERR_GENERIC_ERROR
Definition: rsb.h:537
#define RSB_NUMERICAL_TYPE_FIRST_BLAS
A default numerical matrix type; if possible, not integer one. If no such type is configured in,...
Definition: rsb_types.h:204
#define RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE
The default, blas if possible , numerical type (can be used for declarations).
Definition: rsb_types.h:82
!
! Copyright (C) 2008-2020 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/>.
!
SUBROUTINE blas_sparse_mod_example(res)
USE rsb ! For the second part of the example
IMPLICIT NONE
INTEGER :: res, istat = 0, i, j
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr ! matrix pointer
INTEGER :: A
INTEGER,PARAMETER :: transn = blas_no_trans
INTEGER,PARAMETER :: incx = 1
INTEGER,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,PARAMETER :: nr = 20
INTEGER,PARAMETER :: nc = nr
INTEGER,PARAMETER :: nnz = (nr*(nr+1))/2 ! half the square
INTEGER :: nt = 0
INTEGER :: ic, ir
INTEGER,PARAMETER :: nrhs = 2
INTEGER,PARAMETER :: IA(nnz) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
INTEGER,PARAMETER :: JA(nnz) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
REAL(KIND=8),parameter :: va(nnz) = (/ ((1, ic=1,ir), ir=1,nr ) /) ! (/1, 1, 1/)
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/)
! 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_ptr) ! 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
PROGRAM main
USE rsb, ONLY: rsb_lib_init, rsb_lib_exit, c_ptr, c_null_ptr,&
USE iso_c_binding
IMPLICIT NONE
INTEGER :: res = 0, passed = 0, failed = 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
INTEGER,TARGET::IONE=1
res = rsb_lib_init(io)
CALL blas_sparse_mod_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
@ blas_colmajor
Definition: blas_sparse.h:28
@ blas_rsb_autotune_next_operation
Definition: blas_sparse.h:153
@ blas_rsb_spmv_autotuning_off
Definition: blas_sparse.h:148
@ blas_rsb_spmv_autotuning_on
Definition: blas_sparse.h:147
Definition: rsb_blas_sparse.F90:231
inserts multiple entries
Definition: rsb_blas_sparse.F90:49
multiplication : c <- beta c + alpha A b
Definition: rsb_blas_sparse.F90:126
multiplication : c <- beta c + alpha A b
Definition: rsb_blas_sparse.F90:104
ISO C BINDING interface to rsb_file_mtx_save.
Definition: rsb.F90:738
ISO C BINDING interface to rsb_lib_exit.
Definition: rsb.F90:100
ISO C BINDING interface to rsb_lib_init.
Definition: rsb.F90:42
ISO C BINDING interface to rsb_lib_set_opt.
Definition: rsb.F90:76
Definition: rsb_blas_sparse.F90:29
subroutine usgp(A, pname, istat)
Get a matrix property.
Definition: rsb_blas_sparse.F90:289
subroutine ussp(A, pname, istat)
Set a matrix property. Should be called just after creation, before nonzeroes insertion.
Definition: rsb_blas_sparse.F90:312
subroutine uscr_end(A, istat)
Makes an assembled matrix out of a matrix in build state. After this, it is not possible anymore to i...
Definition: rsb_blas_sparse.F90:267
subroutine usds(A, istat)
Destroys a matrix.
Definition: rsb_blas_sparse.F90:247
subroutine duscr_begin(m, n, A, istat)
Allocates an empty matrix (A) and leaves it in build state.
Definition: rsb_blas_sparse.F90:363
Definition: rsb.F90:11
integer(c_int), parameter rsb_io_want_extra_verbose_interface
See RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE.
Definition: rsb.F90:1070
integer(c_int), parameter rsb_io_want_verbose_tuning
See RSB_IO_WANT_VERBOSE_TUNING.
Definition: rsb.F90:1097
!
! Copyright (C) 2008-2019 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/>.
!
SUBROUTINE rsb_mod_example1(res)
USE rsb
USE iso_c_binding
IMPLICIT NONE
INTEGER ::res
INTEGER,TARGET :: istat = 0, i
INTEGER :: transt = rsb_transposition_n ! Please note that this interface is unfinished
INTEGER :: incx = 1, incy = 1
REAL(KIND=8),target :: alpha = 3, beta = 1
! 1 1
! 1 1
! declaration of VA,IA,JA
INTEGER :: nnz = 4
INTEGER :: nr = 2
INTEGER :: nc = 2
INTEGER :: nrhs = 1
INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
INTEGER :: flags = rsb_flag_noflags
INTEGER,TARGET :: IA(4) = (/0, 1, 1,0/)
INTEGER,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
res = 0
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_ptr)
! 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/)! reference y
IF (istat.NE.rsb_err_no_error) GOTO 9997
istat = rsb_file_mtx_save(mtxap,c_null_ptr)
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"
GOTO 9998
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_example1
SUBROUTINE rsb_mod_example2(res)
USE rsb
USE iso_c_binding
IMPLICIT NONE
INTEGER,TARGET :: errval
INTEGER :: res
INTEGER :: transt = rsb_transposition_n ! no transposition
INTEGER :: incX = 1, incb = 1 ! X, B vectors increment
REAL(KIND=8),target :: alpha = 3,beta = 1
INTEGER :: nnzA = 4, nra = 3, nca = 3 ! nonzeroes, rows, columns of matrix A
INTEGER,TARGET :: IA(4) = (/1, 2, 3, 3/) ! row indices
INTEGER,TARGET :: JA(4) = (/1, 2, 1, 3/) ! column indices
INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double
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
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"
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
PROGRAM main
USE rsb
IMPLICIT NONE
INTEGER :: res = rsb_err_no_error, passed = 0, failed = 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
res = rsb_lib_init(io)
CALL rsb_mod_example1(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
res = rsb_lib_exit(eo)
CALL rsb_mod_example2(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
ISO C BINDING interface to rsb_lib_get_opt.
Definition: rsb.F90:88
ISO C BINDING interface to rsb_mtx_alloc_from_coo_begin.
Definition: rsb.F90:111
ISO C BINDING interface to rsb_mtx_alloc_from_coo_const.
Definition: rsb.F90:204
ISO C BINDING interface to rsb_mtx_alloc_from_coo_end.
Definition: rsb.F90:127
ISO C BINDING interface to rsb_mtx_free.
Definition: rsb.F90:264
ISO C BINDING interface to rsb_mtx_set_vals.
Definition: rsb.F90:643
ISO C BINDING interface to rsb_perror.
Definition: rsb.F90:30
ISO C BINDING interface to rsb_spmv.
Definition: rsb.F90:332
ISO C BINDING interface to rsb_tune_spmm.
Definition: rsb.F90:675
integer(c_int), parameter rsb_flag_default_matrix_flags
See RSB_FLAG_DEFAULT_MATRIX_FLAGS.
Definition: rsb.F90:986
integer(c_int), parameter rsb_err_no_error
See RSB_ERR_NO_ERROR.
Definition: rsb.F90:836
integer(c_int), parameter rsb_io_want_librsb_etime
See RSB_IO_WANT_LIBRSB_ETIME.
Definition: rsb.F90:1094
integer(c_int), parameter rsb_flag_noflags
See RSB_FLAG_NOFLAGS.
Definition: rsb.F90:891
integer(c_int), parameter rsb_flag_symmetric
See RSB_FLAG_SYMMETRIC.
Definition: rsb.F90:951
integer(c_int), parameter rsb_flag_want_column_major_order
See RSB_FLAG_WANT_COLUMN_MAJOR_ORDER.
Definition: rsb.F90:906
integer(c_int), parameter rsb_transposition_n
Definition: rsb.F90:1026
integer(c_signed_char), parameter rsb_numerical_type_double
Definition: rsb.F90:1034