Example programs and code

Examples of usage of librsb in C (<rsb.h> and <blas_sparse.h>), C++ (optional <rsb.hpp>), Fortran (<rsb.F90> and <rsb_blas_sparse.F90>). More...

Files
file	assemble.cpp
	C++ example based on <rsb.hpp> assembling RsbMatrix by pieces.
file	autotune.cpp
	C++ example based on <rsb.hpp> for performance-benchmarking a matrix from file using RsbMatrix.tune_spmm(), for various right-hand side counts.
file	bench.cpp
	C++ example based on <rsb.hpp> for performance-benchmarking a matrix from file using RsbMatrix.spmm(), for various right-hand side counts.
file	build.cpp
	C++ example based on <rsb.hpp> using timings for common matrix operations on RsbMatrix: RsbMatrix.get_coo(), rsb_coo_sort(), rsb_time().
file	example.cpp
	C++ example based on <rsb.hpp> using RsbMatrix.spmm().
file	misc.cpp
	C++ example based on <rsb.hpp> showing various RsbMatrix operations.
file	mtx2bin.cpp
	C++ example based on <rsb.hpp> converting a Matrix Market file into a custom format using RsbMatrix.get_coo().
file	render.cpp
	C++ example based on <rsb.hpp> of invoking the autotuner on an RsbMatrix matrix with RsbMatrix.tune_spmm() and rendering it with RsbMatrix.rndr().
file	span.cpp
	C++ example based on <rsb.hpp> illustrating use of RsbMatrix.file_save(), and std::span-based versions of RsbMatrix.tune_spmm(), RsbMatrix.spmv().
file	twonnz.cpp
	C++ example based on <rsb.hpp> measuring RsbMatrix.spmm() performance of a matrix with only two elements; this is is effectively measuring performance of result vector scaling.
file	autotune.c
	C "RSB autotune" example program based on <rsb.h>. uses rsb_file_mtx_load(), rsb_spmm(), rsb_tune_spmm().
file	backsolve.c
	C triangular solve example program. Uses rsb_spsv(),rsb_tune_spsm(). Based on <rsb.h>.
file	cplusplus.cpp
	C++ program using the C <rsb.h> interface. Uses rsb_tune_spmm(), rsb_spmv().
file	hello-spblas.c
	A first C "hello RSB" example program using a Sparse BLAS interface and <rsb.h>. Uses BLAS_duscr_begin(), BLAS_ussp(), BLAS_usgp(), BLAS_duscr_insert_entries(), BLAS_duscr_end(), BLAS_dusget_element(),BLAS_dusmv(),BLAS_usds().
file	hello.c
	A first "hello RSB" example C program based on <rsb.h>. Uses rsb_lib_set_opt(), rsb_mtx_get_info_str().
file	io-spblas.c
	Example C program using the Sparse BLAS interface and reading from file using rsb_blas_file_mtx_load(), BLAS_usgp(), BLAS_dusmv(), BLAS_usds().
file	power.c
	A toy <rsb.h>-based C program implementing the power method for computing matrix eigenvalues. Uses rsb_spmv().
file	snippets.c
	Collection of C snippets of other examples. Used piecewise the documentation. Not intended to be read as example.
file	transpose.c
	A toy <rsb.h>-based C program showing instantiation, transposition and other operations on a single matrix. Uses rsb_mtx_clone(), rsb_file_mtx_save(), rsb_file_mtx_get_dims(), rsb_file_mtx_load().

Detailed Description

Examples of usage of librsb in C (<rsb.h> and <blas_sparse.h>), C++ (optional <rsb.hpp>), Fortran (<rsb.F90> and <rsb_blas_sparse.F90>).

The following fully working example programs illustrate correct ways of using the library.

• First example in C, using <rsb.h>: examples_hello_c.
• First example in C, using <blas_sparse.h>: examples_hello_spblas_c.
• Autotuning example in C, using <rsb.h>: examples_autotune_c.
• I/O example in C, using <blas_sparse.h>: examples_io_spblas_c.
• Example transposing a matrix in C, using <rsb.h> in C: examples_transpose_c.
• Example showing the power method in C, using <rsb.h> in C: examples_power_c.
• Example in Fortran, using <rsb_blas_sparse.F90>: examples_fortran_F90.
• Example in Fortran, using <rsb.F90>: examples_fortran_rsb_fi_F90.
• Example in C, using <rsb.h>: examples_backsolve_c.
• Misc example snippets in C, using <rsb.h>: examples_snippets_c.
• Benchmark invocation from shell script: examples_bench_sh.

Once installed librsb, the script displayed here (examples/make.sh) should be sufficient to build these examples:

#!/bin/bash
# Example script to build the librsb example programs.
set -e
set -x

srcdir=${srcdir:-`pwd`}
builddir=${builddir:-`pwd`}
LIBRSB_CONFIG=${LIBRSB_CONFIG:-librsb-config}
PKG_CONFIG=pkg-config
WANT_PKGCONFIG="no"

if test x"rsblib" != x"" ; then
for s in `ls ${srcdir}/../rsblib/examples/*.cpp`
do
        p=${builddir}/../rsblib/examples/`basename ${s/.cpp/}`
        rm -f $p
        CXXFLAGS=`${LIBRSB_CONFIG} --cppflags --I_opts`
        LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
        CXX=`${LIBRSB_CONFIG} --cxx`
        LINK=`${LIBRSB_CONFIG} --link`
        o="${p}.o"
        ccmd="$CXX $CXXFLAGS -c $s -o $o"
        lcmd="$LINK $o $LDFLAGS -o $p"
        echo "$ccmd && $lcmd"
        ( $ccmd && $lcmd )
        # one may use a single command, but that's error-prone (may miss libraries):
        #cmd="$CXX $CXXFLAGS $s $LDFLAGS -o $p"
        #echo $cmd
        #$cmd
done
fi

for s in ${srcdir}/*.c
do
        p=`basename ${s/.c/}`
        if test $p == hello-spblas -a x"yes" != x"yes" ; then continue; fi
        if test $p ==    io-spblas -a x"yes" != x"yes" ; then continue; fi
        rm -f $p 
        CFLAGS=`${LIBRSB_CONFIG} --I_opts --cppflags`
        LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
        CC=`${LIBRSB_CONFIG} --cc`
        LINK=`${LIBRSB_CONFIG} --link`
        o="${p}.o"
        ccmd="$CC $CFLAGS -c $s -o $o"
        lcmd="$LINK $o $LDFLAGS -o $p"
        echo "$ccmd && $lcmd"
        ( $ccmd && $lcmd )
        # one may use a single command, but that's error-prone (may miss libraries):
        #cmd="$CC $CFLAGS $s $LDFLAGS -o $p"
        #echo $cmd
        #$cmd
        if test x"${WANT_PKGCONFIG}" != x"no" ; then
                CFLAGS=`${PKG_CONFIG} --cflags librsb`
                LIBS=`${PKG_CONFIG} --libs --static librsb`
                ccmd="$CC $CFLAGS -c $s -o $o"
                lcmd="$LINK $o $LIBS -o $p"
                echo "$ccmd && $lcmd"
                ( $ccmd && $lcmd )
        fi
done

if test x"yes" = x"yes" ; then
if test x"yes" = x"yes" ; then
FP=${srcdir}/fortran_rsb_fi.F90
if test x"yes" = x"yes" ; then
        FP+=\ ${srcdir}/fortran.F90
fi
# activated if you have built the Fortran modules and installed them in the right path.
for s in $FP
do
        p=`basename ${s/.F90/}`
        rm -f $p 
        FCFLAGS=`${LIBRSB_CONFIG} --I_opts --fcflags`
        LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
        FC=`${LIBRSB_CONFIG} --fc`
        LINK=`${LIBRSB_CONFIG} --link`
        o="${p}.o"
        FCLIBS=`${LIBRSB_CONFIG} --fclibs`
        ccmd="$FC $FCFLAGS -c $s -o $o"
        lcmd="$LINK $o $LDFLAGS $FCLIBS -o $p"
        echo "$ccmd && $lcmd"
        ( $ccmd && $lcmd )
done
fi
fi

echo " [*] done building examples!"

examples/hello.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone

 @brief A first "hello RSB" example C program based on <rsb.h>.
   Uses #rsb_lib_set_opt(), #rsb_mtx_get_info_str().

 \include hello.c
*/
#include <rsb.h>        /* librsb header to include */
#include <stdio.h>      /* printf() */

int main(const int argc, char * const argv[])
{
        /*!
          A Hello-RSB program.
         
          This program shows how to use the rsb.h interface correctly to:
         
          - initialize the library using #rsb_lib_init()
          - set library options using #rsb_lib_set_opt()
          - revert such changes 
          - allocate (build) a single sparse matrix in the RSB format
            using #rsb_mtx_alloc_from_coo_const()
          - prints information obtained via #rsb_mtx_get_info_str()
          - multiply the matrix times a vector using #rsb_spmv()
          - deallocate the matrix using #rsb_mtx_free() 
          - finalize the library using #rsb_lib_exit()
         
          In this example, we use #RSB_DEFAULT_TYPE as matrix type.
          This type depends on what was configured at library build time.
         * */
        const rsb_blk_idx_t bs = RSB_DEFAULT_BLOCKING;
        const rsb_blk_idx_t brA = bs, bcA = bs;
        const RSB_DEFAULT_TYPE one = 1;
        const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
        const rsb_nnz_idx_t nnzA = 4;           /* matrix nonzeroes count */
        const rsb_coo_idx_t nrA = 3;            /* matrix rows count */
        const rsb_coo_idx_t ncA = 3;            /* matrix columns count */
        /* nonzero row indices coordinates: */
        const rsb_coo_idx_t IA[] = {0,1,2,2};
        /* nonzero column indices coordinates: */
        const rsb_coo_idx_t JA[] = {0,1,2,2};
        const RSB_DEFAULT_TYPE VA[] = {11,22,32,1};/* values of nonzeroes */
        RSB_DEFAULT_TYPE X[] = { 0, 0, 0 };     /* X vector's array */
        const RSB_DEFAULT_TYPE B[] = { -1, -2, -5 }; /* B vector's array */
        char ib[200];
        struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
        rsb_err_t errval = RSB_ERR_NO_ERROR;

        printf("Hello, RSB!\n");
        printf("Initializing the library...\n");
        if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) != 
                        RSB_ERR_NO_ERROR)
        {
                printf("Error initializing the library!\n");
                goto err;
        }
        printf("Correctly initialized the library.\n");

        printf("Attempting to set the"
               " RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE library option.\n");
        {
                rsb_int_t evi=1; 

                /* Setting a single optional library parameter. */
                errval = rsb_lib_set_opt(
                        RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE, &evi);
                if(errval != RSB_ERR_NO_ERROR)
                {
                        char errbuf[256];
                        rsb_strerror_r(errval,&errbuf[0],sizeof(errbuf));
                        printf("Failed setting the"
                        " RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE"
                        " library option (reason string:\n%s).\n",errbuf);
                        if(errval&RSB_ERRS_UNSUPPORTED_FEATURES)
                        {
                          printf("This error may be safely ignored.\n");
                        }
                        else
                        {
                          printf("Some unexpected error occurred!\n");
                          goto err;
                        }
                }
                else
                {
                        printf("Setting back the "
                                "RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE"
                                " library option.\n");
                        evi = 0;
                        errval = rsb_lib_set_opt(RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE,
                                        &evi);
                        errval = RSB_ERR_NO_ERROR;
                }
        }

        mtxAp = rsb_mtx_alloc_from_coo_const(
                VA,IA,JA,nnzA,typecode,nrA,ncA,brA,bcA,
                RSB_FLAG_NOFLAGS    /* default format will be chosen */
                |RSB_FLAG_DUPLICATES_SUM/* duplicates will be summed */
                        ,&errval);
        if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
        {
                printf("Error while allocating the matrix!\n");
                goto err;
        }
        printf("Correctly allocated a matrix.\n");
        printf("Summary information of the matrix:\n");
        /* print out the matrix summary information  */
        rsb_mtx_get_info_str(mtxAp,"RSB_MIF_MATRIX_INFO__TO__CHAR_P",
                        ib,sizeof(ib));
        printf("%s",ib);
        printf("\n");

        if((errval = 
                rsb_spmv(RSB_TRANSPOSITION_N,&one,mtxAp,B,1,&one,X,1))
                        != RSB_ERR_NO_ERROR )
        {
                printf("Error performing a multiplication!\n");
                goto err;
        }
        printf("Correctly performed a SPMV.\n");
        rsb_mtx_free(mtxAp);
        printf("Correctly freed the matrix.\n");
        if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
                        != RSB_ERR_NO_ERROR)
        {
                printf("Error finalizing the library!\n");
                goto err;
        }
        printf("Correctly finalized the library.\n");
        printf("Program terminating with no error.\n");
        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        printf("Program terminating with error.\n");
        return EXIT_FAILURE;
}

examples/hello-spblas.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone

 @brief A first C "hello RSB" example program using
        a Sparse BLAS interface and <rsb.h>.
        Uses #BLAS_duscr_begin(), #BLAS_ussp(), #BLAS_usgp(),
        #BLAS_duscr_insert_entries(), #BLAS_duscr_end(),
        #BLAS_dusget_element(),#BLAS_dusmv(),#BLAS_usds().

 \include hello-spblas.c
*/
#include <rsb.h>        /* for rsb_lib_init */
#include <blas_sparse.h>        /* Sparse BLAS on the top of librsb */
#include <stdio.h>      /* printf */

int main(const int argc, char * const argv[])
{
        /*!
         * A Hello/Sparse BLAS program.
         *
         * This program shows how to use the blas_sparse.h
         * interface correctly to:
         *
         * - initialize the library using #rsb_lib_init()
         * - allocate (build) a single sparse matrix in the RSB
         *   format using #BLAS_duscr_begin()/#BLAS_duscr_insert_entries()
         *   /#BLAS_duscr_end()
         * - extract one matrix element with #BLAS_dusget_element()
         * - multiply the matrix times a vector using #BLAS_dusmv()
         * - deallocate the matrix using #BLAS_usds() 
         * - finalize the library using
         *   #rsb_lib_exit(#RSB_NULL_EXIT_OPTIONS) 
        */
#ifndef RSB_NUMERICAL_TYPE_DOUBLE   
        printf("'double' type configured out."
        " Please reconfigure the library with it and recompile.\n");
        return EXIT_SUCCESS;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
        blas_sparse_matrix A = blas_invalid_handle; /* handle for A */
        const int nnz = 4;      /* number of nonzeroes of matrix A */
        const int  nr = 3;      /* number of A's rows */
        const int  nc = 3;      /* number of A's columns */
        /* A's nonzero elements row indices (coordinates): */
#ifdef RSB_WANT_LONG_IDX_TYPE 
        const int64_t IA[] = { 0, 1, 2, 2 };
#else /* RSB_WANT_LONG_IDX_TYPE */
        const int   IA[] = { 0, 1, 2, 2 };
#endif /* RSB_WANT_LONG_IDX_TYPE */
        /* A's nonzero elements column indices (coordinates): */
#ifdef RSB_WANT_LONG_IDX_TYPE 
        const int64_t JA[] = { 0, 1, 0, 2 };
#else /* RSB_WANT_LONG_IDX_TYPE */
        const int   JA[] = { 0, 1, 0, 2 };
#endif /* RSB_WANT_LONG_IDX_TYPE */
        /* A's nonzero values (matrix coefficients): */
        double VA[] = { 11.0, 22.0, 13.0, 33.0  };
        /* the X vector's array: */
        double X[] = { 0.0, 0.0, 0.0 };
        /* the B vector's array: */
        const double B[] = { -1.0, -2.0, -2.0 };
        /* the (known) result array: */
        const double AB[] = { 11.0+26.0, 44.0, 66.0+13.0 };
        /* rsb error variable: */
        rsb_err_t errval = RSB_ERR_NO_ERROR;
        int i;

        printf("Hello, RSB!\n");
        /* initialize the library */
        if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) 
                        != RSB_ERR_NO_ERROR)
        {
                goto err;
        }
        printf("Correctly initialized the library.\n");

        /* initialize a matrix descriptor */
        A = BLAS_duscr_begin(nr,nc);
        if( A == blas_invalid_handle )
        {
                goto err;
        }
        
        /* specify properties (e.g.: symmetry)*/
        if( BLAS_ussp(A,blas_lower_symmetric) != 0 )
        {
                goto err;
        }

        /* get properties (e.g.: symmetry) */
        if( BLAS_usgp(A,blas_lower_symmetric) != 1 )
        {
                printf("Symmetry property non set ?!\n");
                goto err;
        }

        /* insert the nonzeroes (here, all at once) */
        if( BLAS_duscr_insert_entries(A, nnz, VA, IA, JA)
                        == blas_invalid_handle)
        {
                goto err;
        }

        /* finalize (allocate) the matrix build  */
        if( BLAS_duscr_end(A) == blas_invalid_handle )
        {
                goto err;
        }
        printf("Correctly allocated a matrix.\n");

        VA[0] = 0.0;
        if( BLAS_dusget_element(A, IA[0], JA[0], &VA[0]) )
        {
                goto err;
        }

        /* a check */
        if( VA[0] != 11.0 )
        {
                goto err;
        }

        /* compute X = X + (-1) * A * B   */
        if(BLAS_dusmv(blas_no_trans,-1,A,B,1,X,1))
        {
                goto err;
        }

        for( i = 0 ; i < nc; ++i )
                if( X[i] != AB[i] )
                {
                        printf("Computed SPMV result seems wrong. Terminating.\n");
                        goto err;
                }
        printf("Correctly performed a SPMV.\n");

        /* deallocate matrix A */
        if( BLAS_usds(A) )
        {
                goto err;
        }
        printf("Correctly freed the matrix.\n");

        /* finalize the library */
        if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
                        != RSB_ERR_NO_ERROR)
        {
                goto err;
        }
        printf("Correctly finalized the library.\n");
        printf("Program terminating with no error.\n");

        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        printf("Program terminating with error.\n");
        return EXIT_FAILURE;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}

examples/autotune.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone

 @brief C "RSB autotune" example program based on <rsb.h>.
   uses #rsb_file_mtx_load(), rsb_spmm(), rsb_tune_spmm().

 \include autotune.c
*/
#include <rsb.h>        /* librsb header to include */
#include <stdio.h>      /* printf() */
#include <ctype.h>      /* isdigit() */
#include <stdlib.h>     /* atoi() */
/* #include "rsb_internals.h" */

static int tune_from_file(char * const filename, rsb_int_t wvat)
{
        struct rsb_mtx_t *mtxMp = NULL;
        /* spmv specific variables */
        const void * alphap = NULL; // equivalent to 1
        const void * betap = NULL; // equivalent to 1
        rsb_flags_t order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER;
        const rsb_coo_idx_t nrhs = 2;  /* number of right hand sides */
        rsb_trans_t transA = RSB_TRANSPOSITION_N; /* transposition */
        rsb_nnz_idx_t ldB = 0;
        rsb_nnz_idx_t ldC = 0;
        /* misc variables */
        rsb_err_t errval = RSB_ERR_NO_ERROR;
        rsb_time_t dt;
        char ib[200];
        const char*is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
        /* misc variables */
        /* input autotuning variables */
        rsb_int_t oitmax = 1 /*15*/;    /* auto-tune iterations */
        rsb_time_t tmax = 0.1;  /* time per autotune operation */
        /* output autotuning variables */
        rsb_flags_t flagsA = RSB_FLAG_NOFLAGS;
        /* int ione = 1; */
        rsb_type_t typecodea [] = RSB_MATRIX_TYPE_CODES_ARRAY;
        int typecodei;

        errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS);

        if( (errval) != RSB_ERR_NO_ERROR )
                goto err;

        errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
        
        /*
        errval = rsb_lib_set_opt(RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE, &ione);
        */

        if( (errval) != RSB_ERR_NO_ERROR )
                goto err;

        printf("Loading matrix from file \"%s\".\n",filename);

        mtxMp = rsb_file_mtx_load(filename, flagsA, typecodea[0], &errval);

        if( (errval) != RSB_ERR_NO_ERROR )
                goto err;

        for( typecodei = 0 ; typecodei < RSB_IMPLEMENTED_TYPES; ++typecodei )
        {
                rsb_type_t typecode = typecodea[typecodei];
                struct rsb_mtx_t *mtxAp = NULL;
                struct rsb_mtx_t *mtxOp = NULL;
                rsb_real_t sf = 0.0;
                rsb_int_t tn = 0;

                sf = 0.0;
                tn = 0;

                printf("Considering %c clone.\n",typecode);
                
                errval = rsb_mtx_clone(&mtxAp, typecode, transA, NULL, mtxMp,
                                flagsA);

                if( (errval) != RSB_ERR_NO_ERROR )
                        goto err;

                printf("Base matrix:\n");
                rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
                printf("%s\n\n",ib);

                dt = -rsb_time();
                errval = rsb_tune_spmm(NULL, &sf, &tn, oitmax, tmax, transA,
                     alphap, mtxAp, nrhs, order, NULL, ldB, betap, NULL, ldC);

                dt += rsb_time();
                if(tn == 0)
                printf("After %lfs, autotuning routine did not find a better"
                        " threads count configuration.\n",dt);
                else
                printf("After %lfs, thread autotuning declared speedup of %lg x,"
                        " when using threads count of %d.\n",dt,sf,tn);
                printf("\n");


                dt = -rsb_time();

                mtxOp = mtxAp;
                errval = rsb_tune_spmm(&mtxAp, &sf, &tn, oitmax, tmax, transA,
                        alphap, NULL, nrhs, order, NULL, ldB, betap, NULL, ldC);
                if( (errval) != RSB_ERR_NO_ERROR )
                        goto err;

                dt += rsb_time();
                if( mtxOp == mtxAp )
                {
                        printf("After %lfs, global autotuning found old matrix optimal,"
                        " with declared speedup %lg x when using %d threads\n",dt,sf,tn);
                }
                else
                {
                        printf("After %lfs, global autotuning declared speedup of %lg x,"
                        " when using threads count of %d and a new matrix:\n",dt,sf,tn);
                        rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
                        printf("%s\n",ib);
                }
                printf("\n");

                /* user is expected to:
                errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
                and use mtxAp in SpMV.
                */
                rsb_mtx_free(mtxAp);
                mtxAp = NULL;
        }
        rsb_mtx_free(mtxMp);
        mtxMp = NULL;

err:
        rsb_perror(NULL,errval);
        if ( errval != RSB_ERR_NO_ERROR )
                printf("Program terminating with error.\n");
        return errval;
}

int main(const int argc, char * const argv[])
{
        /*!
         Autotuning example.
         */
        /* matrix variables */
        struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
        const int bs = RSB_DEFAULT_BLOCKING;
        rsb_coo_idx_t nrA = 500; /* number of rows */
        rsb_coo_idx_t ncA = 500; /* number of cols */
        const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
        const rsb_coo_idx_t rd = 1;/* every rd rows one is non empty */
        const rsb_coo_idx_t cd = 4;/* every cd cols one is non empty */
        rsb_nnz_idx_t nnzA = (nrA/rd)*(ncA/cd); /* nonzeroes */
        rsb_coo_idx_t*IA = NULL;
        rsb_coo_idx_t*JA = NULL;
        RSB_DEFAULT_TYPE*VA = NULL;
        /* spmv specific variables */
        const RSB_DEFAULT_TYPE alpha = 1;
        const RSB_DEFAULT_TYPE beta = 1;
        RSB_DEFAULT_TYPE*Cp = NULL;
        RSB_DEFAULT_TYPE*Bp = NULL;
        rsb_flags_t order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER;
        const rsb_coo_idx_t nrhs = 2;  /* number of right hand sides */
        const rsb_trans_t transA = RSB_TRANSPOSITION_N;
        rsb_nnz_idx_t ldB = nrA;
        rsb_nnz_idx_t ldC = ncA;
        /* misc variables */
        rsb_err_t errval = RSB_ERR_NO_ERROR;
        const size_t so = sizeof(RSB_DEFAULT_TYPE);
        const size_t si = sizeof(rsb_coo_idx_t);
        rsb_time_t dt,odt;
        rsb_int_t t;
        const rsb_int_t tt = 100;       /* will repeat spmv tt times */
        char ib[200];
        const char*is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
        /* misc counters */
        rsb_coo_idx_t ci;
        rsb_coo_idx_t ri;
        rsb_coo_idx_t ni;
        rsb_int_t nrhsi;
        /* misc variables */
        rsb_time_t etime = 0.0;
        /* input autotuning variables */
        const rsb_int_t oitmax = 15;    /* auto-tune iterations */
        const rsb_time_t tmax = 0.1;    /* time per autotune operation */
        /* input/output autotuning variables */
        rsb_int_t tn = 0;       /* threads number */
        /* output autotuning variables */
        rsb_real_t sf = 0.0;    /* speedup factor obtained from auto tuning */
        const rsb_int_t wvat = 1; /* want verbose autotuning; see 
                documentation of RSB_IO_WANT_VERBOSE_TUNING */

        if(argc > 1 && !isdigit(argv[1][0]) )
        {
                errval = tune_from_file(argv[1],wvat);
                goto ret;
        }

        if(argc > 1)
        {
                nrA = ncA = atoi(argv[1]);
                if ( nrA < RSB_MIN_MATRIX_DIM || (nrA > (RSB_MAX_MATRIX_DIM) ))
                        goto err;

                nnzA = (nrA/rd)*(ncA/cd);
                ldB = nrA;
                ldC = ncA;
        }

        printf("Creating %d x %d matrix with %d nonzeroes.\n",(int)nrA,
                (int)ncA, (int)nnzA);

        IA = calloc(nnzA, si);
        JA = calloc(nnzA, si);
        VA = calloc(nnzA, so);
        Bp = calloc(nrhs*ncA ,so);
        Cp = calloc(nrhs*nrA ,so);

        if( ! ( VA && IA && JA && Bp && Cp ) )
                goto err;

        for(nrhsi=0;nrhsi<nrhs;++nrhsi)
                for(ci=0;ci<ncA/cd;++ci)
                        Bp[nrhsi*ldC+ci] = 1.0;

        for(nrhsi=0;nrhsi<nrhs;++nrhsi)
                for(ri=0;ri<nrA/rd;++ri)
                        Cp[nrhsi*ldC+ri] = 1.0;

        ni = 0;

        for(ci=0;ci<ncA/cd;++ci)
                for(ri=0;ri<nrA/rd;++ri)
                {
                        VA[ni] = nrA * ri + ci,
                        IA[ni] = ri;
                        JA[ni] = ci;
                        ni++;
                }

        if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS))
                        != RSB_ERR_NO_ERROR) goto err;

        errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
        if( (errval) != RSB_ERR_NO_ERROR )
        {
                printf("Error setting option!\n");
                goto err;
        }

        mtxAp = rsb_mtx_alloc_from_coo_const(
                VA,IA,JA,nnzA,typecode,nrA,ncA,bs,bs,
                RSB_FLAG_NOFLAGS,&errval);

        /* VA, IA, JA are not necessary anymore */
        free(VA);
        free(IA);
        free(JA);
        VA = NULL;
        IA = NULL;
        JA = NULL;

        if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
                goto err;

        printf("Allocated matrix of %zd nonzeroes:\n",(size_t)nnzA);
        rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
        printf("%s\n\n",ib);

        dt = - rsb_time();
        for(t=0;t<tt;++t)
                /* 
                   If nrhs == 1, the following is equivalent to
                   rsb_spmv(transA,alphap,mtxAp,Bp,1,betap,Cp,1);
                */
                rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
        dt += rsb_time();
        odt = dt;
        printf("Before auto-tuning, %d multiplications took %lfs.\n",tt,dt);

        printf("Threads autotuning (may take more than %lfs)...\n",
                        oitmax*tmax);
        dt = -rsb_time();
        errval = rsb_tune_spmm(NULL, &sf, &tn, oitmax, tmax, transA,
                        &alpha, mtxAp, nrhs, order, Bp, ldB, &beta, Cp, ldC);
        dt += rsb_time();
        if(errval != RSB_ERR_NO_ERROR)
                goto err;

        if(tn == 0)
        printf("After %lfs, autotuning routine did not find a better"
                        " threads count configuration.\n",dt);
        else
        printf("After %lfs, autotuning routine declared speedup of %lg x,"
                        " when using threads count of %d.\n",dt,sf,tn);

        errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
        if(errval != RSB_ERR_NO_ERROR)
                goto err;

        rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
        printf("%s\n",ib);

        dt = -rsb_time();
        for(t=0;t<tt;++t)
                /*rsb_spmv(transA,&alpha,mtxAp,Bp,1,&beta,Cp,1);*/
                rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
        dt += rsb_time();
        printf("After threads auto-tuning, %d multiplications took %lfs"
                        "  --  effective speedup of %lg x\n",tt,dt,odt/dt);
        odt = dt;


        tn = 0; /* this will restore default threads count */
        errval = rsb_lib_set_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
        if(errval != RSB_ERR_NO_ERROR)
                goto err;
        errval = rsb_lib_get_opt(RSB_IO_WANT_EXECUTING_THREADS,&tn);
        if(errval != RSB_ERR_NO_ERROR)
                goto err;

        printf("Matrix autotuning (may take more than %lfs; using %d"
                        " threads )...\n", oitmax*tmax, tn);

        /* A negative tn will request also threads autotuning: */
        /* tn = -tn; */

        dt = -rsb_time();
        errval = rsb_tune_spmm(&mtxAp, &sf, &tn, oitmax, tmax, transA,
                &alpha,  NULL, nrhs, order, Bp, ldB, &beta, Cp, ldC);
        dt += rsb_time();

        if(errval != RSB_ERR_NO_ERROR)
                goto err;

        if(tn == 0)
        printf("After %lfs, autotuning routine did not find a better"
                        " threads count configuration.\n",dt);
        else
        printf("After %lfs, autotuning routine declared speedup of %lg x,"
                        " when using threads count of %d.\n",dt,sf,tn);

        rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
        printf("%s\n",ib);

        dt = -rsb_time();
        for(t=0;t<tt;++t)
                /*rsb_spmv(transA,&alpha,mtxAp,Bp,1,&beta,Cp,1);*/
                rsb_spmm(transA,&alpha,mtxAp,nrhs,order,Bp,ldB,&beta,Cp,ldC);
        dt += rsb_time();
        printf("After threads auto-tuning, %d multiplications took %lfs"
                        "  --  further speedup of %lg x\n",tt,dt,odt/dt);

        rsb_mtx_free(mtxAp);
        free(Cp);
        free(Bp);


        errval = rsb_lib_get_opt(RSB_IO_WANT_LIBRSB_ETIME,&etime);
        if(errval == RSB_ERR_UNSUPPORTED_FEATURE)
        {
                printf("librsb timer-based profiling is not supported in "
                "this build. If you wish to have it, re-configure librsb "
                "with its support. So you can safely ignore the error you"
                " might just have seen printed out on screen.\n");
                errval = RSB_ERR_NO_ERROR;
        }
        else
        if(etime) /* This will only work if enabled at configure time. */
                printf("Elapsed program time is %5.2lfs\n",etime);

ret:
        if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
                        !=RSB_ERR_NO_ERROR)
                goto err;
        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        printf("Program terminating with error.\n");
        return EXIT_FAILURE;
}

examples/io-spblas.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone

 @brief Example C program using the Sparse BLAS interface
        and reading from file using \ref rsb_blas_file_mtx_load(),
        \ref BLAS_usgp(), \ref BLAS_dusmv(), \ref BLAS_usds().

 \include io-spblas.c
*/
#include <rsb.h>        /* for rsb_lib_init */
#include <blas_sparse.h>
#include <stdio.h>
        
int main(const int argc, char * const argv[])
{
#ifndef RSB_NUMERICAL_TYPE_DOUBLE   
        printf("Skipping a test because of 'double' type opted out.\n");
        return EXIT_SUCCESS;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
        blas_sparse_matrix A = blas_invalid_handle;
        const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
        const rsb_char_t * filename = argc > 1 ? argv[1] : "pd.mtx";

        printf("Hello, RSB!\n");
        if((rsb_perror(NULL,
                rsb_lib_init(RSB_NULL_INIT_OPTIONS)))!=RSB_ERR_NO_ERROR)
        {
                printf("Error while initializing the library.\n");
                goto err;
        }

        printf("Correctly initialized the library.\n");

        A = rsb_blas_file_mtx_load(filename, typecode );
        if( A == blas_invalid_handle )
        {
                printf("Error while loading matrix %s from file.\n",
                                filename);
                goto err;
        }

        printf("Correctly loaded and allocated a matrix"
                        " from file %s.\n",filename);

        if( BLAS_usgp(A,blas_symmetric) == 1 )
                printf("Matrix is symmetric\n");

        if( BLAS_usgp(A,blas_hermitian) == 1 )
                printf("Matrix is hermitian\n");

        printf("Now SPMV with NULL vectors will be attempted,"
                        " resulting in an error (so don't worry).\n");

        if(BLAS_dusmv(blas_no_trans,-1,A,NULL,1,NULL,1))
        {
                printf("Correctly detected an error condition.\n");
                goto okerr;
        }

        printf("No error detected ?\nIf you see this line printed out,"
                " please report as a bug, because the above NULL pointers"
                " should have been detected\n");
        return EXIT_FAILURE;
okerr:
        printf("Program correctly recovered from intentional"
                        " error condition.\n");
        if(BLAS_usds(A))
        {
                printf("Error while freeing the matrix!\n");
                goto err;
        }

        printf("Correctly freed the matrix.\n");
err:
        if(rsb_perror(NULL,
                rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))!=RSB_ERR_NO_ERROR)
        {
                printf("Failed finalizing the library.\n");
                goto ferr;
        }

        printf("Correctly finalized the library.\n");
        return EXIT_SUCCESS;
ferr:
        return EXIT_FAILURE;
#endif /* RSB_NUMERICAL_TYPE_DOUBLE */
}

examples/transpose.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 @file
 @author Michele Martone

 @brief A toy <rsb.h>-based C program showing instantiation,
 transposition and other operations on a single matrix.
 Uses \ref rsb_mtx_clone(), \ref rsb_file_mtx_save(),
  \ref rsb_file_mtx_get_dims(), \ref rsb_file_mtx_load().

 \ingroup rsb_doc_examples

 \include transpose.c
*/
#include <rsb.h>
#include <stdio.h>      /* printf */

int main(const int argc, char * const argv[])
{
        struct rsb_mtx_t *mtxAp = NULL;
        const rsb_blk_idx_t brA = RSB_DEFAULT_BLOCKING,
                            bcA = RSB_DEFAULT_BLOCKING;
        rsb_nnz_idx_t nnzA = 4;
        rsb_coo_idx_t  nrA = 3;
        rsb_coo_idx_t  ncA = 3;
        const rsb_coo_idx_t    IA[] = { 0, 1, 2, 0 };
        const rsb_coo_idx_t    JA[] = { 0, 1, 2, 2 };
        const RSB_DEFAULT_TYPE VA[] = { 11, 22, 33, 13 };
        RSB_DEFAULT_TYPE XV[] = { 0,0,0,0,0,0 };
        rsb_coo_idx_t  vl = 0;
        const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
        rsb_err_t errval = RSB_ERR_NO_ERROR;

        /* library initialization */
        if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
        {
                return EXIT_FAILURE;
        }

        /* allocation */
        mtxAp = rsb_mtx_alloc_from_coo_const(
                        VA,IA,JA,nnzA,typecode,nrA,ncA,
                        brA,bcA,RSB_FLAG_NOFLAGS,NULL);
        if(!mtxAp)
        {
                return EXIT_FAILURE;
        }

        /* printout */
        if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
        {
                if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
                        goto err;
        }
        
        /* matrix transposition */
        if( RSB_ERR_NO_ERROR != (errval =
                rsb_mtx_clone(&mtxAp,RSB_NUMERICAL_TYPE_SAME_TYPE,
                RSB_TRANSPOSITION_T,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS)))
        {
                goto err;
        }

        /* printout */
        if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
        {
                if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
                        goto err;
        }

        rsb_mtx_free(mtxAp);

        /* doing the same after load from file */
        mtxAp = rsb_file_mtx_load("pd.mtx",
                RSB_FLAG_NOFLAGS,typecode,NULL);
        if(!mtxAp)
        {
                return EXIT_FAILURE;
        }

        /* printout */
        if(RSB_ERR_NO_ERROR!=(errval = rsb_file_mtx_save(mtxAp,NULL)))
        {
                if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
                        goto err;
        }

        /* one can see dimensions in advance, also */
        if(RSB_ERR_NO_ERROR!=(errval =
                rsb_file_mtx_get_dims("pd.mtx",&nrA,&ncA,&nnzA,NULL)))
        {
                if(errval != RSB_ERR_UNSUPPORTED_FEATURE)
                        goto err;
        }

        /* A matrix can be rendered to Postscript. */
        {
                if(RSB_ERR_NO_ERROR!=(errval =
                rsb_mtx_rndr("pd.eps",mtxAp,512,512,RSB_MARF_EPS_B)))
                        goto err;
        }

        rsb_mtx_free(mtxAp);

        /* also vectors can be loaded */
        if(RSB_ERR_NO_ERROR!=(errval = 
                rsb_file_vec_load("vf.mtx",typecode,NULL,&vl )))
                goto err;
        /* we expect vf.mtx to be 6 rows long */
        if( vl != 6 )
        {
                goto err;
        }

        if(RSB_ERR_NO_ERROR!=(errval = 
                rsb_file_vec_load("vf.mtx",typecode,XV, NULL )))
                goto err;

        /* matrices can be rendered from file to a pixelmap as well */
        {
                unsigned char pixmap[3*2*2];

                if(RSB_ERR_NO_ERROR!=(errval =
                rsb_file_mtx_rndr(pixmap,"pd.mtx",2,2,2,RSB_MARF_RGB)))
                        goto err;
        }

        if(RSB_ERR_NO_ERROR != rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
        {
                goto err;
        }
        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        return EXIT_FAILURE;
}

examples/power.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 @file
 @author Michele Martone

 @brief A toy <rsb.h>-based C program implementing the power
        method for computing matrix eigenvalues. Uses #rsb_spmv().
 \ingroup rsb_doc_examples

 \include power.c
*/

#include <stdio.h>      // printf
#include <math.h>       // sqrt
#include <stdlib.h>     // calloc
#include <rsb.h>

int main(const int argc, char * const argv[])
{
        const int WANT_VERBOSE = 0;
        struct rsb_mtx_t *mtxAp = NULL;
        const int bs = RSB_DEFAULT_BLOCKING;
        const int br = bs, bc = bs; /* bs x bs blocked */
        const rsb_nnz_idx_t nnzA = 4;
        const rsb_coo_idx_t  nrA = 3;
        const rsb_coo_idx_t  ncA = 3;
        rsb_int_t it = 0;
        const rsb_int_t maxit = 100;
        const rsb_coo_idx_t    IA[] = { 0, 1, 2, 0 };
        const rsb_coo_idx_t    JA[] = { 0, 1, 2, 2 };
        const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[] = { 11, 22, 33, 13 };
        const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE ZERO = 0;
        int i;
        rsb_err_t errval = 0;

        RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE norm = 0.0, /* nu */
                oldnorm = 1.0, /* oldnorm */
                *b1 = NULL, *b2 = NULL,
                *bnow = NULL, *bnext = NULL;/* b1 and b2 aliases */
        rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
        size_t ds = 0;
        /* tolerance */
        const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE tol = 1e-14;

        /* library initialization */
        if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
                return EXIT_FAILURE;

        /* allocation */
        mtxAp = rsb_mtx_alloc_from_coo_const(VA,IA,JA,nnzA,
                        typecode,nrA,ncA,br,bc,RSB_FLAG_NOFLAGS,NULL);
        if(!mtxAp)
                return EXIT_FAILURE;

        ds = (nrA)*sizeof(RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE);
        b1 = calloc(1,ds);
        b2 = calloc(1,ds);

        if(! (b1 && b2))
        {
                errval = RSB_ERR_ENOMEM;
                goto err;
        }

        for( i = 0; i < nrA; ++i )
                b1[i] = 1;

        bnow = b1, bnext = b2;/* b,b' */

        while( fabs(norm-oldnorm) > tol && it<maxit )
        {
                ++ it;
                oldnorm = norm;
                /* b'<-Ab */
                if(( rsb_spmv(RSB_TRANSPOSITION_N,NULL,mtxAp,bnow,
                        1,&ZERO,bnext,1)) != RSB_ERR_NO_ERROR )
                        goto err;
                /* nu<-||Ab||^2 */
                norm = 0;
                for(i=0;i<nrA;++i) 
                        norm += bnext[i]*bnext[i];
                /* nu<-||Ab|| */
                norm = sqrt(norm);
                norm = 1.0/norm;
                /* b'<- Ab / ||Ab|| */
                for(i=0;i<nrA;++i)
                        bnext[i] *= norm;
                norm = 1.0/norm;
                printf("it:%d norm:%lg norm diff:%lg\n",it,norm,norm-oldnorm);

                {void *tmp=bnow;bnow=bnext;bnext=tmp;/* pointers swap */}
                if(WANT_VERBOSE)
                {
                        printf("norm:%lg\n",norm);
                        if(isinf(norm))
                        /* isinf is a C99 feature (need correct
                         * compilation flags) */
                                goto err;

                        for(i=0;i<2;++i)
                                printf("x[%d]=%lg\n",i,((double*)bnext)[i]);
                }
        }
        /* the biggest eigenvalue should be in bnow */

        rsb_mtx_free(mtxAp);
        free(b1);
        free(b2);
        if(rsb_lib_exit(RSB_NULL_EXIT_OPTIONS)!=RSB_ERR_NO_ERROR)
                goto err;
        if( it == maxit )
        {
                printf("ERROR: hit iterations limit without convergence!");
                errval=RSB_ERR_GENERIC_ERROR;
        }
        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        return EXIT_FAILURE;
}

examples/fortran.F90:

! 
! Copyright (C) 2008-2021 Michele Martone
! 
! This file is part of librsb.
! 
! librsb is free software; you can redistribute it and/or modify it
! under the terms of the GNU Lesser General Public License as published
! by the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
! 
! librsb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
! License for more details.
! 
! You should have received a copy of the GNU Lesser General Public
! License along with librsb; see the file COPYING.
! If not, see <http://www.gnu.org/licenses/>.
! 

!> @file
!! @brief Sparse BLAS-based usage:
!!  uscr_begin(), usgp(), ussp(), usmv(), ...
!! \include fortran.F90

      SUBROUTINE blas_sparse_mod_example(res)
      USE blas_sparse
      USE rsb ! For the second part of the example and RSB_IDX_KIND
      IMPLICIT NONE
      INTEGER :: i, j
      INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0, res
      TYPE(c_ptr),TARGET :: mtxap = c_null_ptr ! matrix pointer
      INTEGER :: a
      INTEGER,PARAMETER :: transn = blas_no_trans
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
      REAL(KIND=8),PARAMETER :: alpha = 3
! Symmetric (declared via lower triangle) matrix based example, e.g.:
! 1 0
! 1 1
      ! declaration of VA,IA,JA 
      !INTEGER,PARAMETER :: nr = 100
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nr = 20
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nc = nr
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nnz = (nr*(nr+1))/2 ! half the square
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 1
      INTEGER(KIND=RSB_IDX_KIND) :: nt = 0
      INTEGER :: ic, ir
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: ia(nnz) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: ja(nnz) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
      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

      ! Example loading a matrix from file and measuring SPMM.
      SUBROUTINE blas_sparse_io_example(res)
      USE blas_sparse
      USE rsb ! For rsb_blas_file_mtx_load
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
      INTEGER :: j
      INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0
      INTEGER :: a
      INTEGER,PARAMETER :: transn = blas_no_trans
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
      COMPLEX(KIND=8),PARAMETER :: alpha = 3
      INTEGER(KIND=RSB_IDX_KIND) :: nr
      INTEGER(KIND=RSB_IDX_KIND) :: nc
      INTEGER(KIND=RSB_IDX_KIND) :: nz
      INTEGER(KIND=RSB_IDX_KIND) :: st
      INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 4
      COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: x
      COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: y
      CHARACTER(KIND=C_SIGNED_CHAR,LEN=7),TARGET :: filename = "pd.mtx"//c_null_char
      INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double_complex
      REAL(KIND=C_DOUBLE) :: mvt,mmt,omt
      INTEGER(KIND=C_INT),TARGET::izero=0

      res = 0

      a = rsb_blas_file_mtx_load(c_loc(filename), typecode);
      IF (a.EQ.blas_invalid_handle) GOTO 9997

      CALL usgp(a,blas_num_rows,nr)
      CALL usgp(a,blas_num_cols,nc)
      CALL usgp(a,blas_num_nonzeros,nz)

      print*,"Read matrix ",filename(1:6)," ",nr,"x",nc,":",nz
      CALL usgp(a,blas_general,st)
      IF (st .EQ. 1) print*,"Matrix has no symmetry"
      CALL usgp(a,blas_upper_symmetric,st)
      IF (st .EQ. 1) print*,"Matrix is upper symmetric"
      CALL usgp(a,blas_upper_hermitian,st)
      IF (st .EQ. 1) print*,"Matrix is upper hermitian"
      ! ...

      IF (istat.NE.0) GOTO 9997
      WRITE(*,'(a,i0)') "Using NRHS=",nrhs
      ALLOCATE( x(nc,nrhs))
      ALLOCATE( y(nr,nrhs))
      x = 1.0

      y = 0.0
      mvt = -rsb_time()
      DO j = 1, nrhs
        CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
      END DO
      IF (istat.NE.0) GOTO 9997
      mvt = mvt + rsb_time()
      WRITE(*,'(a,e12.4,a)') "Repeated USMV took ",mvt," s"
      
      y = 0.0
      mmt = -rsb_time()
      CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
      IF (istat.NE.0) GOTO 9997
      mmt = mmt + rsb_time()
      WRITE(*,'(a,e12.4,a)') "A single USMM took ",mmt," s"

      WRITE(*,'(a,g11.4,a)')"USMM-to-USMV speed ratio is is ", mvt/mmt, "x"

      print*,"Call auto-tuning routine.."
      ! Change IZERO value to 1 to get verbose tuning again.
      res = rsb_lib_set_opt(rsb_io_want_verbose_tuning,c_loc(izero))
      IF (res.NE.0) GOTO 9997
      CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
      IF (istat.NE.0) GOTO 9997
      CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
      IF (istat.NE.0) GOTO 9997

      print*,"Repeat measurement."
      y = 0.0
      omt = -rsb_time()
      CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
      IF (istat.NE.0) GOTO 9997
      omt = omt + rsb_time()
      WRITE(*,'(a,e12.4,a)') "Tuned USMM took ",omt," s"
      WRITE(*,'(a,g11.4,a)')"Tuned-to-untuned speed ratio is is ",mmt/omt,"x"
      
      GOTO 9998
9997      res = -1
9998      CONTINUE
      CALL usds(a,istat)
      IF (istat.NE.0) res = -1
      end SUBROUTINE blas_sparse_io_example

      PROGRAM main
      USE rsb, ONLY: rsb_lib_init, rsb_lib_exit, c_ptr, c_null_ptr,&
       & rsb_io_want_extra_verbose_interface,rsb_io_want_verbose_tuning,&
       & rsb_lib_set_opt,rsb_idx_kind
      ! USE blas_sparse, only: RSB_BLAS_IDX_KIND ! only if using long indices
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER :: passed = 0, failed = 0
      INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
      !TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
      !TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
      ! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
      TYPE(c_ptr),PARAMETER :: eo = c_null_ptr
      TYPE(c_ptr),PARAMETER :: io = c_null_ptr
      INTEGER(KIND=C_INT),TARGET::ione=1
      res = rsb_lib_init(io)
      res = rsb_lib_set_opt(rsb_io_want_verbose_tuning,c_loc(ione))
      
      CALL blas_sparse_mod_example(res)
      IF (res.LT.0) failed = failed + 1
      IF (res.EQ.0) passed = passed + 1

      CALL blas_sparse_io_example(res)
      IF (res.LT.0) failed = failed + 1
      IF (res.EQ.0) passed = passed + 1

      res = rsb_lib_exit(eo)
      
      print *, "FAILED:", failed
      print *, "PASSED:", passed
      IF (failed .GT. 0) THEN
       stop 1
      END IF
      END PROGRAM

examples/fortran_rsb_fi.F90:

! 
! Copyright (C) 2008-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/>.
! 

!> @file.
!! @brief RSB.F90-based usage:
!!        rsb_mtx_alloc_from_coo_const(),
!!        rsb_tune_spmm(),
!!        rsb_file_mtx_save(),
!!        rsb_spmv(),
!!        ...
!! \include fortran_rsb_fi.F90

      SUBROUTINE rsb_mod_example1(res)
! Example using an unsymmetric matrix specified as COO, and autotuning, built at once.
      USE rsb
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER ::res
      INTEGER,TARGET :: istat = 0, i
      INTEGER :: transt = rsb_transposition_n !
      INTEGER(KIND=RSB_IDX_KIND) :: incx = 1, incy = 1
      REAL(KIND=8),TARGET :: alpha = 3, beta = 1
! 1 1
! 1 1
      ! declaration of VA,IA,JA 
      INTEGER(KIND=RSB_IDX_KIND) :: nnz = 4
      INTEGER(KIND=RSB_IDX_KIND) :: nr = 2
      INTEGER(KIND=RSB_IDX_KIND) :: nc = 2
      INTEGER(KIND=RSB_IDX_KIND) :: nrhs = 1
      INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
      INTEGER :: flags = rsb_flag_noflags 
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ia(4) = (/0, 1, 1,0/)
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ja(4) = (/0, 0, 1,1/)
      REAL(KIND=8),TARGET :: va(4) = (/1,1,1,1/)
      REAL(KIND=8),TARGET :: x(2) = (/1, 1/)! reference x 
      REAL(KIND=8),TARGET :: cy(2) = (/9, 9/)! reference cy after 
      REAL(KIND=8),TARGET :: y(2) = (/3, 3/)! y will be overwritten
      TYPE(c_ptr),TARGET :: mtxap = c_null_ptr ! matrix pointer
      REAL(KIND=8) :: tmax = 2.0 ! tuning max time
      INTEGER :: titmax = 2 ! tuning max iterations
      INTEGER,TARGET :: ont = 0     ! optimal number of threads
      !TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
      !TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
      ! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
      TYPE(c_ptr),PARAMETER :: eo = c_null_ptr
      TYPE(c_ptr),PARAMETER :: io = c_null_ptr
      INTEGER,TARGET :: errval

      res = 0

      errval = rsb_lib_init(io)
      IF (errval.NE.rsb_err_no_error) GOTO 9997

      mtxap = rsb_mtx_alloc_from_coo_const(c_loc(va),c_loc(ia),c_loc(ja)&
       &,nnz,&
       & rsb_numerical_type_double,nr,nc,1,1,flags,c_loc(istat))

      IF (istat.NE.rsb_err_no_error) GOTO 9997

      istat = rsb_file_mtx_save(mtxap,c_null_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/)! restoring reference y (rsb_tune_spmm has overwritten it)
      IF (istat.NE.rsb_err_no_error) GOTO 9997
      
      istat = rsb_file_mtx_save(mtxap,c_null_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"

      IF ( res .NE.rsb_err_no_error) GOTO 9997
      GOTO 9998
9997      res = -1
9998      CONTINUE
      mtxap = rsb_mtx_free(mtxap)
      IF (istat.NE.rsb_err_no_error) res = -1 
! 9999      CONTINUE
      errval=rsb_lib_exit(eo)                 ! librsb finalization
      IF (errval.NE.rsb_err_no_error)&
       & stop "error calling rsb_lib_exit"
      print *, "rsb module fortran test is ok"

      istat = rsb_perror(c_null_ptr,istat)
      end SUBROUTINE rsb_mod_example1

      SUBROUTINE rsb_mod_example2(res)
! Example using an unsymmetric matrix specified as COO, and autotuning, built piecewise.
      USE rsb
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER,TARGET :: errval
      INTEGER :: res
      INTEGER :: transt = rsb_transposition_n  ! no transposition
      INTEGER(KIND=RSB_IDX_KIND) :: incx = 1, incb = 1        ! X, B vectors increment
      REAL(KIND=8),TARGET :: alpha = 3,beta = 1
      INTEGER(KIND=RSB_IDX_KIND) :: nnza = 4, nra = 3, nca = 3     ! nonzeroes, rows, columns of matrix A
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ia(4) = (/1, 2, 3, 3/)  ! row    indices
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ja(4) = (/1, 2, 1, 3/)  ! column indices
      INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double
      INTEGER :: flags =rsb_flag_default_matrix_flags+rsb_flag_symmetric
      REAL(KIND=8),TARGET :: va(4) = (/11.0, 22.0, 13.0, 33.0/) ! coefficients
      REAL(KIND=8),TARGET :: x(3) = (/   0,    0,    0/)
      REAL(KIND=8),TARGET :: b(3) = (/-1.0, -2.0, -2.0/)
      TYPE(c_ptr),TARGET  :: mtxap = c_null_ptr
      TYPE(c_ptr)  :: mtxapp = c_null_ptr
      REAL(KIND=8),TARGET :: etime = 0.0
      !TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
      !TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
      ! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
      TYPE(c_ptr),PARAMETER :: eo = c_null_ptr
      TYPE(c_ptr),PARAMETER :: io = c_null_ptr

      errval = rsb_lib_init(io)                ! librsb initialization
      IF (errval.NE.rsb_err_no_error) &
       & stop "error calling rsb_lib_init"
#if defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ < 5)
#define RSB_SKIP_BECAUSE_OLD_COMPILER 1
#endif
#ifndef RSB_SKIP_BECAUSE_OLD_COMPILER
      mtxap = rsb_mtx_alloc_from_coo_begin(nnza,typecode,nra,nca,flags,&
       & c_loc(errval)) ! begin matrix creation
      errval = rsb_mtx_set_vals(mtxap,&
       & c_loc(va),c_loc(ia),c_loc(ja),nnza,flags) ! insert some nonzeroes
      mtxapp = c_loc(mtxap) ! Old compilers like e.g.: Gfortran 4.4.7 will NOT compile this.
      IF (errval.NE.rsb_err_no_error) &
       & stop "error calling rsb_mtx_set_vals"
      errval = rsb_mtx_alloc_from_coo_end(mtxapp)                   ! end matrix creation
      IF (errval.NE.rsb_err_no_error) &
       & stop "error calling rsb_mtx_alloc_from_coo_end"
      errval = rsb_spmv(transt,c_loc(alpha),mtxap,c_loc(x),&
       & incx,c_loc(beta),c_loc(b),incb) ! X := X + (3) * A * B 
      IF (errval.NE.rsb_err_no_error)&
       & stop "error calling rsb_spmv"
      mtxap = rsb_mtx_free(mtxap)                                 ! destroy matrix

      ! The following is optional and depends on configure options, so it is allowed to fail
      errval = rsb_lib_get_opt(rsb_io_want_librsb_etime,c_loc(etime))
      IF (errval.EQ.rsb_err_no_error)&
       & print*,"Time spent in librsb is:",etime
      ! IF (errval.NE.0)STOP "error calling rsb_lib_get_opt" 
      errval = rsb_err_no_error

      IF (errval.NE.rsb_err_no_error) &
       & stop "error calling rsb_mtx_free"
#else
      print*,"You have an old Fortran compiler not supporting C_LOC."
      print*,"Skipping a part of the test"
#endif
      errval=rsb_lib_exit(eo)                 ! librsb finalization
      IF (errval.NE.rsb_err_no_error)&
       & stop "error calling rsb_lib_exit"
      print *, "rsb module fortran test is ok"
      res = errval
      end SUBROUTINE rsb_mod_example2

      SUBROUTINE rsb_mod_example3(res)
! Example using a symmetric matrix specified as CSR, and autotuning, built at once.
      USE rsb
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER ::res
      INTEGER,TARGET :: istat = 0, i
      INTEGER :: transt = rsb_transposition_n !
      INTEGER(KIND=RSB_IDX_KIND) :: incx = 1, incy = 1
      REAL(KIND=8),TARGET :: alpha = 4, beta = 1
! 11 21 
! 21 22 
      ! declaration of VA,IP,JA 
      INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nnz = 3
      INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nr = 2
      INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nc = 2
      INTEGER(KIND=RSB_IDX_KIND), PARAMETER :: nrhs = 1
      INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
      INTEGER :: flags = rsb_flag_noflags + rsb_flag_symmetric + &
       &                 rsb_flag_fortran_indices_interface        !  optional flags
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ip(3) = (/1, 2, 4/)
      INTEGER(KIND=RSB_IDX_KIND),TARGET :: ja(3) = (/1, 1, 2/)
      REAL(KIND=8),TARGET :: va(3) = (/11,21,22/) ! lower triangle coefficients
      REAL(KIND=8),TARGET :: x(2) = (/1, 2/)! reference x 
      REAL(KIND=8),TARGET :: cy(2) = (/215.0, 264.0/)! reference cy after 
      REAL(KIND=8),TARGET :: y(2) = (/3, 4/)! y will be overwritten
      TYPE(c_ptr),TARGET :: mtxap = c_null_ptr ! matrix pointer
      REAL(KIND=8) :: tmax = 2.0 ! tuning max time
      INTEGER :: titmax = 2 ! tuning max iterations
      INTEGER,TARGET :: ont = 0     ! optimal number of threads
      TYPE(c_ptr),PARAMETER :: eo = c_null_ptr
      TYPE(c_ptr),PARAMETER :: io = c_null_ptr
      INTEGER,TARGET :: errval

      errval = rsb_lib_init(io)                ! librsb initialization
      IF (errval.NE.rsb_err_no_error) &
       & stop "error calling rsb_lib_init"

      res = 0
      mtxap = rsb_mtx_alloc_from_csr_const(c_loc(va),c_loc(ip),c_loc(ja)&
       &,nnz,rsb_numerical_type_double,nr,nc,1,1,flags,c_loc(istat))

      IF (istat.NE.rsb_err_no_error) GOTO 9997

      istat = rsb_file_mtx_save(mtxap,c_null_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, 4/)! restoring reference y (rsb_tune_spmm has overwritten it)
      IF (istat.NE.rsb_err_no_error) GOTO 9997
      
      istat = rsb_file_mtx_save(mtxap,c_null_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)
      print *, y
      IF (istat.NE.rsb_err_no_error) GOTO 9997
      DO i = 1, 2
            IF (y(i).NE.cy(i)) print *, "type=d dims=2x2 sym=s diag=g &
      &blocks=1x1 usmv alpha= 4 beta= 1 incx=1 incy=1 trans=n is not ok"
            IF (y(i).NE.cy(i)) GOTO 9997
      END DO
      print*,"type=d dims=2x2 sym=s diag=g blocks=1x1 usmv alpha= 4&
       & beta= 1 incx=1 incy=1 trans=n is ok"
      GOTO 9998

      errval=rsb_lib_exit(eo)                 ! librsb finalization
      IF (errval.NE.rsb_err_no_error)&
       & stop "error calling rsb_lib_exit"
      print *, "rsb module fortran test is ok"
9997      res = -1
9998      CONTINUE
      mtxap = rsb_mtx_free(mtxap)
      IF (istat.NE.rsb_err_no_error) res = -1 
! 9999      CONTINUE
      istat = rsb_perror(c_null_ptr,istat)
      end SUBROUTINE rsb_mod_example3

      PROGRAM main
      USE rsb
      IMPLICIT NONE
      INTEGER :: res = rsb_err_no_error, passed = 0, failed = 0

      CALL rsb_mod_example1(res)
      IF (res.LT.0) failed = failed + 1
      IF (res.EQ.0) passed = passed + 1

      CALL rsb_mod_example2(res)
      IF (res.LT.0) failed = failed + 1
      IF (res.EQ.0) passed = passed + 1

      CALL rsb_mod_example3(res)
      IF (res.LT.0) failed = failed + 1
      IF (res.EQ.0) passed = passed + 1
      
      print *, "FAILED:", failed
      print *, "PASSED:", passed
      IF (failed.GT.0) THEN
       stop 1
      END IF
      END PROGRAM

examples/backsolve.c:

/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone

 @brief C triangular solve example program.
   Uses #rsb_spsv(),#rsb_tune_spsm().
   Based on <rsb.h>.

 \include hello.c
*/
#include <rsb.h>        /* librsb header to include */
#include <stdio.h>      /* printf() */

int main(const int argc, char * const argv[])
{
        /*!
          A Hello-RSB program.
         
          This program shows how to use the rsb.h interface correctly to:
         
          - initialize the library using #rsb_lib_init()
          - allocate (build) a single sparse matrix in the RSB format
            using #rsb_mtx_alloc_from_coo_const(), with implicit diagonal
          - print information obtained via #rsb_mtx_get_info_str()
          - multiply the triangular matrix using #rsb_spmv()
          - autotune the matrix for rsb_spsv with #rsb_tune_spsm()
          - solve the triangular system using #rsb_spsv()
          - deallocate the matrix using #rsb_mtx_free() 
          - finalize the library using #rsb_lib_exit()
         
          In this example, we use #RSB_DEFAULT_TYPE as matrix type.
          This type depends on what was configured at library build time.
         * */
        const int bs = RSB_DEFAULT_BLOCKING;
        const int brA = bs, bcA = bs;
        const RSB_DEFAULT_TYPE one = 1;
        const rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
        const rsb_nnz_idx_t nnzA = 7;           /* matrix nonzeroes count */
        const rsb_coo_idx_t nrA = 6;            /* matrix rows count */
        const rsb_coo_idx_t ncA = 6;            /* matrix columns count */
        /* nonzero row indices coordinates: */
        const rsb_coo_idx_t IA[] = {1,2,3,4,5,6,1};
        /* nonzero column indices coordinates: */
        const rsb_coo_idx_t JA[] = {1,2,3,4,5,6,6};
        const RSB_DEFAULT_TYPE VA[] = {1,1,1,1,1,1,1};/*values of nonzeroes*/
        RSB_DEFAULT_TYPE X[] = { 0,0,0,0,0,0 }; /* X vector's array */
        const RSB_DEFAULT_TYPE B[] = { 1,1,1,1,1,1 }; /* B */
        struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
        char ib[200];
        int i;
        rsb_err_t errval = RSB_ERR_NO_ERROR;
        const rsb_int_t wvat = 1; /* want verbose autotuning; see 
                  documentation of RSB_IO_WANT_VERBOSE_TUNING */

        printf("Hello, RSB!\n");
        printf("Initializing the library...\n");
        if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) != 
                        RSB_ERR_NO_ERROR)
        {
                printf("Error initializing the library!\n");
                goto err;
        }
        printf("Correctly initialized the library.\n");

        errval = rsb_lib_set_opt(RSB_IO_WANT_VERBOSE_TUNING, &wvat );
        if( (errval) != RSB_ERR_NO_ERROR )
        {
                printf("Error setting option!\n");
                goto err;
        }

        mtxAp = rsb_mtx_alloc_from_coo_const(
                VA,IA,JA,nnzA,typecode,nrA,ncA,brA,bcA,
                RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS /* force rsb */
                | RSB_FLAG_DUPLICATES_SUM/* sum dups */
                | RSB_FLAG_UNIT_DIAG_IMPLICIT/* ask diagonal implicit */
                | RSB_FLAG_TRIANGULAR /* need triangle for spsv */
                | RSB_FLAG_FORTRAN_INDICES_INTERFACE /* treat indices as 1-based */
                , &errval);
        if((!mtxAp) || (errval != RSB_ERR_NO_ERROR))
        {
                printf("Error while allocating the matrix!\n");
                goto err;
        }
        printf("Correctly allocated a matrix with %ld nonzeroes.\n",
                (long int)nnzA);
        printf("Summary information of the matrix:\n");
        /* print out the matrix summary information  */
        rsb_mtx_get_info_str(mtxAp,"RSB_MIF_MATRIX_INFO__TO__CHAR_P",
                        ib,sizeof(ib));
        printf("%s",ib);
        printf("\nMatrix printout:\n");

        rsb_file_mtx_save(mtxAp, NULL);

        if((errval = 
                rsb_spmv(RSB_TRANSPOSITION_N,&one,mtxAp,B,1,&one,X,1))
                        != RSB_ERR_NO_ERROR )
        {
                printf("Error performing a multiplication!\n");
                goto err;
        }
        printf("\nWe have a unitary vector:\n");
        rsb_file_vec_save(NULL, typecode, B, nrA);
        
        printf("\nMultiplying matrix by unitary vector we get:\n");
        rsb_file_vec_save(NULL, typecode, X, nrA);

        errval = rsb_tune_spsm(&mtxAp, NULL, NULL, 0, 0, RSB_TRANSPOSITION_N,
                &one, NULL, 1, RSB_FLAG_WANT_COLUMN_MAJOR_ORDER, NULL, nrA,
                 NULL, NULL, nrA);
        if( (errval) != RSB_ERR_NO_ERROR )
        {
                printf("Error performing autotuning!\n");
                goto err;
        }

        if((errval = rsb_spsv(RSB_TRANSPOSITION_N,&one,mtxAp,X,1,X,1))
                        != RSB_ERR_NO_ERROR )
        {
                printf("Error performing triangular solve!\n");
                goto err;
        }

        printf("\nBacksolving we should get a unitary vector:\n");
        rsb_file_vec_save(NULL, typecode, X, nrA);

        for(i=0;i<nrA;++i)
                if(X[i]!=one)
                {
                        printf("Warning! Result vector not unitary!:\n");
                        errval = RSB_ERR_INTERNAL_ERROR;
                        goto err;
                }

        printf("All done.\n");
        rsb_mtx_free(mtxAp);
        printf("Correctly freed the matrix.\n");
        if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS))
                        != RSB_ERR_NO_ERROR)
        {
                printf("Error finalizing the library!\n");
                goto err;
        }
        printf("Correctly finalized the library.\n");
        printf("Program terminating with no error.\n");
        return EXIT_SUCCESS;
err:
        rsb_perror(NULL,errval);
        printf("Program terminating with error.\n");
        return EXIT_FAILURE;
}

examples/bench.sh:

#!/bin/bash
#
# Copyright (C) 2008-2021 Michele Martone
# 
# This file is part of librsb.
# 
# librsb is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
# 
# librsb is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
# License for more details.
# 
# You should have received a copy of the GNU Lesser General Public
# License along with librsb; see the file COPYING.
# If not, see <http://www.gnu.org/licenses/>.

# Benchmark rsbench and postprocess results.
# \ingroup rsb_doc_examples
# @file
# @author Michele Martone
# @brief Benchmark invocation from shell script.
#
set -e
set -x
which rsbench
BRF=test.rpr
# invoke rsbench and produce a performance record, using all types and one thread
rsbench -oa -Ob --bench --lower 100 --as-symmetric \
        --types ':' -n 1 \
        --notranspose --compare-competitors \
        --verbose --verbose \
        --write-performance-record=${BRF}

# examine tuning renderings (produced by --verbose --verbose)
ls -ltr ${BRF/.rpr/}-tuning*

# convert the performance record to text form
rsbench --read-performance-record ${BRF} > ${BRF/rpr/txt}
ls -ltr ${BRF/rpr/txt}

# convert the performance record to LaTeX table document form
RSB_PR_WLTC=2 RSB_PR_SR=0 \
 rsbench --read-performance-record ${BRF} > ${BRF/rpr/tex}
which latex || exit 0 # to compile LaTeX document
which kpsepath || exit 0 # to check LaTeX packages
find `kpsepath tex | sed 's/!!//g;s/:/\n/g;'` -name sciposter.cls || exit 0 # need sciposter class, usually in texlive-science

# convert the LaTeX table into a DVI (may as well use pdflatex for PDF)
latex -interaction=batchmode -file-line-error ${BRF/rpr/tex}

# convert the performance record to GNUPLOT plots
which gnuplot || exit 0
RSB_PRD_STYLE_PLT_PFN=${BRF/rpr/} RSB_PRD_STYLE_PLT_FMT=1 RSB_PR_SR=2 \
 rsbench --read-performance-record ${BRF} > ${BRF/rpr/gnu}

# convert the GNUPLOT plots into PDF
ls -ltr ${BRF/rpr/gnu}
gnuplot ${BRF/rpr/gnu}

# convert the performance record to GNUPLOT plots, different way
RSB_PRD_STYLE_PLT_PFN=${BRF/rpr/} RSB_PRD_STYLE_PLT_FMT=  RSB_PR_SR=2 \
 rsbench --read-performance-record ${BRF} > ${BRF/rpr/gnu}
gnuplot ${BRF/rpr/gnu}
ls -ltr ${BRF/.rpr/}*.png
ls -ltr ${BRF/.rpr/}*.eps
ls -ltr ${BRF} ${BRF/rpr/tex} ${BRF/rpr/dvi} ${BRF/rpr/txt}

# clean up
rm ${BRF/rpr/}{aux,log,out,gnu}
rm ${BRF/.rpr/}*{.png,.eps}
rm      ${BRF} ${BRF/rpr/tex} ${BRF/rpr/dvi} ${BRF/rpr/txt}
exit

Most of the snippets in the documentation come from examples/snippets.c.