Provided by: librsb-dev_1.2.0-rc7-5_amd64 
NAME
librsb - rsb-examples - Examples of usage of librsb.
DESCRIPTION
SYNOPSIS
Detailed Description
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-2015 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 This is a first "hello RSB" example program.
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(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;
rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
rsb_err_t errval = RSB_ERR_NO_ERROR;
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!0);
printf("Initializing the library...0);
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) !=
RSB_ERR_NO_ERROR)
{
printf("Error initializing the library!0);
goto err;
}
printf("Correctly initialized the library.0);
printf("Attempting to set the"
" RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE library option.0);
{
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:3s).0,errbuf);
if(errval&RSB_ERRS_UNSUPPORTED_FEATURES)
{
printf("This error may be safely ignored.0);
}
else
{
printf("Some unexpected error occurred!0);
goto err;
}
}
else
{
printf("Setting back the "
"RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE"
" library option.0);
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!0);
goto err;
}
printf("Correctly allocated a matrix.0);
printf("Summary information of the matrix:0);
/* 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("0);
if((errval =
rsb_spmv(RSB_TRANSPOSITION_N,&one,mtxAp,B,1,&one,X,1))
!= RSB_ERR_NO_ERROR )
{
printf("Error performing a multiplication!0);
goto err;
}
printf("Correctly performed a SPMV.0);
rsb_mtx_free(mtxAp);
printf("Correctly freed the matrix.0);
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
printf("Error finalizing the library!0);
goto err;
}
printf("Correctly finalized the library.0);
printf("Program terminating with no error.0);
return 0;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.0);
return -1;
}
/*
Copyright (C) 2008-2015 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 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 */
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.0);
return 0;
#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: */
rsb_err_t errval = RSB_ERR_NO_ERROR;
int i;
printf("Hello, RSB!0);
/* initialize the library */
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly initialized the library.0);
/* 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 ?!0);
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.0);
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.0);
goto err;
}
printf("Correctly performed a SPMV.0);
/* deallocate matrix A */
if( BLAS_usds(A) )
{
goto err;
}
printf("Correctly freed the matrix.0);
/* finalize the library */
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly finalized the library.0);
printf("Program terminating with no error.0);
return 0;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.0);
return -1;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}
/*
Copyright (C) 2008-2015 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 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 */
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.0);
return 0;
#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: */
rsb_err_t errval = RSB_ERR_NO_ERROR;
int i;
printf("Hello, RSB!0);
/* initialize the library */
if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly initialized the library.0);
/* 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 ?!0);
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.0);
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.0);
goto err;
}
printf("Correctly performed a SPMV.0);
/* deallocate matrix A */
if( BLAS_usds(A) )
{
goto err;
}
printf("Correctly freed the matrix.0);
/* finalize the library */
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!= RSB_ERR_NO_ERROR)
{
goto err;
}
printf("Correctly finalized the library.0);
printf("Program terminating with no error.0);
return 0;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.0);
return -1;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}
/*
Copyright (C) 2008-2015 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 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;
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_SPBLAS_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
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.0,typecode);
errval = rsb_mtx_clone(&mtxAp, typecode, transA, NULL, mtxMp,
flagsA);
if( (errval) != RSB_ERR_NO_ERROR )
goto err;
printf("Base matrix:0);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s0,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.0,dt);
else
printf("After %lfs, thread autotuning declared speedup of %lg x,"
" when using threads count of %d.0,dt,sf,tn);
printf("0);
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 threads0,dt,sf,tn);
}
else
{
printf("After %lfs, global autotuning declared speedup of %lg x,"
" when using threads count of %d and a new matrix:0,dt,sf,tn);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s0,ib);
}
printf("0);
/* 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;
goto ret;
ret:
return 0;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.0);
return -1;
}
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_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
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;
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 = nrA;
rsb_nnz_idx_t ldC = ncA;
/* misc variables */
rsb_err_t errval = RSB_ERR_NO_ERROR;
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_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 */
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]) )
return tune_from_file(argv[1],wvat);
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.0,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++;
}
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 );
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:0,(size_t)nnzA);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s0,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.0,tt,dt);
printf("Threads autotuning (may take more than %lfs)...0,
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.0,dt);
else
printf("After %lfs, autotuning routine declared speedup of %lg x,"
" when using threads count of %d.0,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("%s0,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 x0,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 )...0, 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.0,dt);
else
printf("After %lfs, autotuning routine declared speedup of %lg x,"
" when using threads count of %d.0,dt,sf,tn);
rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
printf("%s0,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 x0,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.0);
errval = RSB_ERR_NO_ERROR;
}
else
if(etime) /* This will only work if enabled at configure time. */
printf("Elapsed program time is %5.2lfs0,etime);
if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
!=RSB_ERR_NO_ERROR)
goto err;
return 0;
err:
rsb_perror(NULL,errval);
printf("Program terminating with error.0);
return -1;
}
/*
Copyright (C) 2008-2015 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 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>
int main(const int argc, char * const argv[])
{
#ifndef RSB_NUMERICAL_TYPE_DOUBLE
printf("Skipping a test because of 'double' type opted out.0);
return 0;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
blas_sparse_matrix A = blas_invalid_handle;
rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
rsb_char_t * filename = argc > 1 ? argv[1] : "pd.mtx";
printf("Hello, RSB!0);
if((rsb_perror(NULL,
rsb_lib_init(RSB_NULL_INIT_OPTIONS)))!=RSB_ERR_NO_ERROR)
{
printf("Error while initializing the library.0);
goto err;
}
printf("Correctly initialized the library.0);
A = rsb_load_spblas_matrix_file_as_matrix_market(filename,
typecode );
if( A == blas_invalid_handle )
{
printf("Error while loading matrix %s from file.0,
filename);
goto err;
}
printf("Correctly loaded and allocated a matrix"
" from file %s.0,filename);
if( BLAS_usgp(A,blas_symmetric) == 1 )
printf("Matrix is symmetric0);
if( BLAS_usgp(A,blas_hermitian) == 1 )
printf("Matrix is hermitian0);
printf("Now SPMV with NULL vectors will be attempted,"
" resulting in an error (so don't worry).0);
if(BLAS_dusmv(blas_no_trans,-1,A,NULL,1,NULL,1))
{
printf("Correctly detected an error condition.0);
goto okerr;
}
printf("No error detected ?0f you see this line printed out,"
" please report as a bug, because the above NULL pointers"
" should have been detected0);
return -1;
okerr:
printf("Program correctly recovered from intentional"
" error condition.0);
if(BLAS_usds(A))
{
printf("Error while freeing the matrix!0);
goto err;
}
printf("Correctly freed the matrix.0);
err:
if(rsb_perror(NULL,
rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))!=RSB_ERR_NO_ERROR)
{
printf("Failed finalizing the library.0);
goto ferr;
}
printf("Correctly finalized the library.0);
return 0;
ferr:
return -1;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}
/*
Copyright (C) 2008-2015 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-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;
rsb_blk_idx_t brA = RSB_DEFAULT_BLOCKING, bcA=RSB_DEFAULT_BLOCKING;
rsb_err_t errval = RSB_ERR_NO_ERROR;
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;
rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
/* library initialization */
if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
{
return -1;
}
/* allocation */
mtxAp = rsb_mtx_alloc_from_coo_const(
VA,IA,JA,nnzA,typecode,nrA,ncA,
brA,bcA,RSB_FLAG_NOFLAGS,NULL);
if(!mtxAp)
{
return -1;
}
/* 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 -1;
}
/* 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 expecy 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 0;
err:
rsb_perror(NULL,errval);
return -1;
}
/*
Copyright (C) 2008-2015 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-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[])
{
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 };
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE ZERO = 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 -1;
/* allocation */
mtxAp = rsb_mtx_alloc_from_coo_const(VA,IA,JA,nnzA,
typecode,nrA,ncA,br,bc,RSB_FLAG_NOFLAGS,NULL);
if(!mtxAp)
return -1;
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:%lg0,it,norm,norm-oldnorm);
{void *tmp=bnow;bnow=bnext;bnext=tmp;/* pointers swap */}
if(WANT_VERBOSE)
{
printf("norm:%lg0,norm);
if(isinf(norm))
/* isinf is a C99 feature (need correct
* compilation flags) */
goto err;
for(i=0;i<2;++i)
printf("x[%d]=%lg0,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 0;
err:
rsb_perror(NULL,errval);
return -1;
}
!
! Copyright (C) 2008-2016 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 blas_sparse
USE rsb ! For the second part of the example
IMPLICIT NONE
INTEGER :: res, istat = 0, i
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 :: 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) = (/((1), ir=1,nc)/) ! reference x ! (/1, 1/)
REAL(KIND=8),PARAMETER :: cy(nr) = (/((alpha+alpha*nr), ir=1,nr)/) ! reference cy after ! (/9, 9/)
REAL(KIND=8) :: y(nr) = (/((alpha), ir=1,nr)/) ! 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
CALL usmv(transn,alpha,a,x,incx,y,incy,istat)
IF (istat.NE.0) GOTO 9997
!
DO i = 1, nr
IF (y(i).NE.cy(i)) print *, "first check results are not ok"
IF (y(i).NE.cy(i)) GOTO 9997
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
CALL usmv(transn,alpha,a,x,incx,y,incy,istat)
CALL usmm(blas_colmajor,transn,1, alpha,a,x,nr,y,nc,istat) ! Equivalent to the above (as long as incx=incy=1).
CALL usmm(blas_colmajor,transn,1,-alpha,a,x,nr,y,nc,istat) ! Subtract the last usmm call contribution.
IF (istat.NE.0) GOTO 9997
!
DO i = 1, nr
IF (y(i).NE.cy(i)) print *,"second check results are not ok"
IF (y(i).NE.cy(i)) GOTO 9997
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,&
& rsb_io_want_extra_verbose_interface,rsb_io_want_verbose_tuning,&
& rsb_lib_set_opt
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)
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
res = rsb_lib_exit(eo)
print *, "FAILED:", failed
print *, "PASSED:", passed
IF (failed .GT. 0) THEN
stop 1
END IF
END PROGRAM
!
! Copyright (C) 2008-2016 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
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
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
Author
librsb was written by Michele Martone; this documentation has been generated by Doxygen.
SEE ALSO
rsb-examples rsb-spblas.h rsb.h