Example programs and code

Examples of usage of librsb. More...

Examples of usage of librsb.

The following fully working example programs illustrate correct ways of using the library. The script displayed here should be sufficient to build them.

#!/bin/bash
# Script to build the librsb example programs.

LIBRSB_CONFIG=${LIBRSB_CONFIG:-librsb-config}

for s in *.c
do
        p=${s/.c/}
        rm -f $p 
        CFLAGS=`${LIBRSB_CONFIG} --I_opts`
        LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
        CC=`${LIBRSB_CONFIG} --cc`
        cmd="$CC $CFLAGS $s $LDFLAGS -o $p"
        echo $cmd
        $cmd
done

if test x"yes" = x"yes" ; then
# activated if you have built the Fortran modules and installed them in the right path.
for s in *.F90
do
        p=${s/.F90/}
        rm -f $p 
        CFLAGS=`${LIBRSB_CONFIG} --I_opts`
        LDFLAGS=`${LIBRSB_CONFIG} --ldflags --extra_libs`
        FC=`${LIBRSB_CONFIG} --fc`
        cmd="$FC $CFLAGS $s $LDFLAGS -o $p"
        echo $cmd
        $cmd
done
fi

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone
 @brief This is a first "hello RSB" example program.

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

int main(const int argc, char * const argv[])
{
        /*!
          A Hello-RSB program.
         
          This program shows how to use the rsb.h interface correctly to:
         
          - initialize the library using #rsb_lib_init()
          - set library options using #rsb_lib_set_opt()
          - revert such changes 
          - allocate (build) a single sparse matrix in the RSB format
            using #rsb_mtx_alloc_from_coo_const()
          - prints information obtained via #rsb_mtx_get_info_str()
          - multiply the matrix times a vector using #rsb_spmv()
          - deallocate the matrix using #rsb_mtx_free() 
          - finalize the library using #rsb_lib_exit(RSB_NULL_EXIT_OPTIONS) 
         
          In this example, we use #RSB_DEFAULT_TYPE as matrix type.
          This type depends on what was configured at library build time.
         * */
        struct rsb_mtx_t *mtxAp = NULL; /* matrix structure pointer */
        const int bs = RSB_DEFAULT_BLOCKING;
        const int brA = bs, bcA = bs;
        const RSB_DEFAULT_TYPE one = 1;
        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!\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;
}

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone
 @brief This is a first "hello RSB" example program using 
        a Sparse BLAS interface.

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

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

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone
 @brief This is a first "hello RSB" example program using 
        a Sparse BLAS interface.

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

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

/*

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 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 void * alphap = NULL; // in rsb_spmv, NULL alphap defaults to 1
        const void * betap = NULL; // in rsb_spmv, NULL betap defaults to 1
        const 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; /* 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*const is = "RSB_MIF_MATRIX_INFO__TO__CHAR_P";
        /* misc variables */
        /* input autotuning variables */
        const rsb_int_t oitmax = 1 /*15*/;      /* auto-tune iterations */
        const rsb_time_t tmax = 0.1;    /* time per autotune operation */
        /* output autotuning variables */
        const rsb_flags_t flagsA = RSB_FLAG_NOFLAGS;
        const int ione = 1;
        const 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 \"%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;
        const 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; /* 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*const 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",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:\n",(size_t)nnzA);
        rsb_mtx_get_info_str(mtxAp,is,ib,sizeof(ib));
        printf("%s\n\n",ib);

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

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

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

        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;
}

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 \ingroup rsb_doc_examples
 @file
 @author Michele Martone
 @brief This is an example program using a Sparse BLAS interface
        and reading from file using the RSB library.

 \include io-spblas.c
*/
#include <rsb.h>        /* for rsb_lib_init */
#include <blas_sparse.h>
#include <stdio.h>
#include <stdlib.h>     // EXIT_SUCCESS
        
int main(const int argc, char * const argv[])
{
#ifndef RSB_NUMERICAL_TYPE_DOUBLE   
        printf("Skipping a test because of 'double' type opted out.\n");
        return EXIT_SUCCESS;
#else /* RSB_NUMERICAL_TYPE_DOUBLE */
        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!\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_load_spblas_matrix_file_as_matrix_market(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 */
}

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 @file
 @author Michele Martone
 @brief A toy program showing instantiation, transposition and other
 operations on a single matrix.
 \ingroup rsb_doc_examples

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

int main(const int argc, char * const argv[])
{
        struct rsb_mtx_t *mtxAp = NULL;
        rsb_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 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;
}

/*

Copyright (C) 2008-2020 Michele Martone

This file is part of librsb.

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

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

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

*/
/*!
 @file
 @author Michele Martone
 @brief A toy program implementing the power method
        for computing matrix eigenvalues.
 \ingroup rsb_doc_examples

 \include power.c
*/

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

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

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

      SUBROUTINE blas_sparse_mod_example(res)
      USE blas_sparse
      USE rsb ! For the second part of the example
      IMPLICIT NONE
      INTEGER :: res, istat = 0, i, j
      TYPE(c_ptr),TARGET :: mtxap = c_null_ptr ! matrix pointer
      INTEGER :: a
      INTEGER,PARAMETER :: transn = blas_no_trans
      INTEGER,PARAMETER :: incx = 1
      INTEGER,PARAMETER :: incy = 1
      REAL(KIND=8),PARAMETER :: alpha = 3
! Symmetric (declared via lower triangle) matrix based example, e.g.:
! 1 0
! 1 1
      ! declaration of VA,IA,JA 
      !INTEGER,PARAMETER :: nr = 100
      INTEGER,PARAMETER :: nr = 20
      INTEGER,PARAMETER :: nc = nr
      INTEGER,PARAMETER :: nnz = (nr*(nr+1))/2 ! half the square
      INTEGER :: nt = 0
      INTEGER :: ic, ir
      INTEGER,PARAMETER :: nrhs = 2
      INTEGER,PARAMETER :: ia(nnz) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
      INTEGER,PARAMETER :: ja(nnz) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
      REAL(KIND=8),PARAMETER :: va(nnz) = (/ ((1, ic=1,ir), ir=1,nr ) /) ! (/1, 1, 1/)
      REAL(KIND=8) :: x(nc,nrhs) = reshape((/((1), ic=1,nc*nrhs)/),[nc,nrhs]) ! reference x ! (/1, 1/)
      REAL(KIND=8),PARAMETER :: cy(nr,nrhs) = reshape((/((alpha+alpha*nr), ir=1,nr*nrhs)/),[nr,nrhs]) ! reference cy after ! (/9, 9/)
      REAL(KIND=8) :: y(nr,nrhs) = reshape((/((alpha), ir=1,nr*nrhs)/),[nr,nrhs]) ! y will be overwritten ! (/3, 3/)
      ! First example part: pure blas_sparse code.
      res = 0
      CALL duscr_begin(nr,nc,a,res)
      IF (res.NE.0) GOTO 9999
      CALL ussp(a,blas_lower_symmetric,istat)
      IF (istat.NE.0) GOTO 9997
      CALL ussp(a,blas_rsb_spmv_autotuning_on,istat) ! (experimental) turns auto-tuning + thread setting on
      IF (istat.NE.0) print *,"autotuning returned nonzero:", istat &
       &," ...did you enable autotuning ?"
      !
      ! First style example 
      CALL uscr_insert_entries(a,nnz,va,ia,ja,istat)
      IF (istat.NE.0) GOTO 9997
      CALL uscr_end(a,istat)
      IF (istat.NE.0) GOTO 9997
      ! CALL ussp(A,blas_rsb_duplicates_sum,istat)
      ! CALL uscr_insert_entries(A,nnz,VA,IA,JA,istat) ! uncomment this to activate add of coefficients to pattern
      CALL usgp(a,blas_rsb_spmv_autotuning_on,nt)  ! (experimental)
      IF (nt.NE.0) print*,"autotuner chose ",nt," threads"
      CALL ussp(a,blas_rsb_spmv_autotuning_off,istat) ! (experimental) turns auto-tuning + thread setting off
      IF (istat.NE.0) GOTO 9997

      DO j = 1, nrhs
        CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
      END DO
      IF (istat.NE.0) GOTO 9997
      !
      DO j = 1, nrhs
      DO i = 1, nr
        IF (y(i,j).NE.cy(i,j)) print *, "first check results are not ok"
        IF (y(i,j).NE.cy(i,j)) GOTO 9997
      END DO
      END DO
      !
      y(:,:) = alpha ! reset
      !
      ! Second style example 
      CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
      IF (istat.NE.0) GOTO 9997
      DO j = 1, nrhs
        CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
      END DO
      CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat) ! Equivalent to the above (as long as incx=incy=1).
      CALL usmm(blas_colmajor,transn,nrhs,-alpha,a,x,nr,y,nc,istat) ! Subtract the last usmm call contribution.
      IF (istat.NE.0) GOTO 9997
      !
      DO j = 1, nrhs
      DO i = 1, nr
        IF (y(i,j).NE.cy(i,j)) print *,"second check results are not ok"
        IF (y(i,j).NE.cy(i,j)) GOTO 9997
      END DO
      END DO
      !
      print *, "check results are ok"
      
      ! Second part of the example: access to the rsb.h interface via
      ! the ISO C Binding interface.
      mtxap = rsb_blas_get_mtx(a) ! get pointer to rsb structure (as in the rsb.h API)
      IF(nr.LT.5) istat = rsb_file_mtx_save(mtxap,c_null_ptr) ! write to stdout (only if matrix small enough)

      GOTO 9998
9997      res = -1
9998      CONTINUE
      CALL usds(a,istat)
      IF (istat.NE.0) res = -1
9999      CONTINUE
      end SUBROUTINE blas_sparse_mod_example

      PROGRAM main
      USE rsb, ONLY: rsb_lib_init, rsb_lib_exit, c_ptr, c_null_ptr,&
       & 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-2019 Michele Martone
! 
! This file is part of librsb.
! 
! librsb is free software; you can redistribute it and/or modify it
! under the terms of the GNU Lesser General Public License as published
! by the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
! 
! librsb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
! License for more details.
! 
! You should have received a copy of the GNU Lesser General Public
! License along with librsb; see the file COPYING.
! If not, see <http://www.gnu.org/licenses/>.
! 
      SUBROUTINE rsb_mod_example1(res)
      USE rsb
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER ::res
      INTEGER,TARGET :: istat = 0, i
      INTEGER :: transt = rsb_transposition_n ! Please note that this interface is unfinished
      INTEGER :: incx = 1, incy = 1
      REAL(KIND=8),TARGET :: alpha = 3, beta = 1
! 1 1
! 1 1
      ! declaration of VA,IA,JA 
      INTEGER :: nnz = 4
      INTEGER :: nr = 2
      INTEGER :: nc = 2
      INTEGER :: nrhs = 1
      INTEGER :: order = rsb_flag_want_column_major_order ! rhs layout
      INTEGER :: flags = rsb_flag_noflags 
      INTEGER,TARGET :: ia(4) = (/0, 1, 1,0/)
      INTEGER,TARGET :: ja(4) = (/0, 0, 1,1/)
      REAL(KIND=8),TARGET :: va(4) = (/1,1,1,1/)
      REAL(KIND=8),TARGET :: x(2) = (/1, 1/)! reference x 
      REAL(KIND=8),TARGET :: cy(2) = (/9, 9/)! reference cy after 
      REAL(KIND=8),TARGET :: y(2) = (/3, 3/)! y will be overwritten
      TYPE(c_ptr),TARGET :: mtxap = c_null_ptr ! matrix pointer
      REAL(KIND=8) :: tmax = 2.0 ! tuning max time
      INTEGER :: titmax = 2 ! tuning max iterations
      INTEGER,TARGET :: ont = 0     ! optimal number of threads

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

      IF (istat.NE.rsb_err_no_error) GOTO 9997

      istat = rsb_file_mtx_save(mtxap,c_null_ptr)

      ! Structure autotuning:
      istat = rsb_tune_spmm(c_loc(mtxap),c_null_ptr,c_null_ptr,titmax,&
       & tmax,&
       & transt,c_loc(alpha),c_null_ptr,nrhs,order,c_loc(x),nr,&
       & c_loc(beta),c_loc(y),nc)

      IF (istat.NE.rsb_err_no_error) GOTO 9997

      ! Thread count autotuning:
      istat = rsb_tune_spmm(c_null_ptr,c_null_ptr,c_loc(ont),titmax,&
       & tmax,&
       & transt,c_loc(alpha),mtxap,nrhs,order,c_loc(x),nr,c_loc(beta),&
       & c_loc(y),nc)
      print *, "Optimal number of threads:", ont

      y(:) = (/3, 3/)! reference y 
      IF (istat.NE.rsb_err_no_error) GOTO 9997
      
      istat = rsb_file_mtx_save(mtxap,c_null_ptr)
      IF (istat.NE.rsb_err_no_error) GOTO 9997

      istat = rsb_spmv(transt,c_loc(alpha),mtxap,c_loc(x),incx,&
       & c_loc(beta),c_loc(y),incy)
      IF (istat.NE.rsb_err_no_error) GOTO 9997
      DO i = 1, 2
            IF (y(i).NE.cy(i)) print *, "type=d dims=2x2 sym=g diag=g &
      &blocks=1x1 usmv alpha= 3 beta= 1 incx=1 incy=1 trans=n is not ok"
            IF (y(i).NE.cy(i)) GOTO 9997
      END DO
      print*,"type=d dims=2x2 sym=g diag=g blocks=1x1 usmv alpha= 3&
       & beta= 1 incx=1 incy=1 trans=n is ok"
      GOTO 9998
9997      res = -1
9998      CONTINUE
      mtxap = rsb_mtx_free(mtxap)
      IF (istat.NE.rsb_err_no_error) res = -1 
! 9999      CONTINUE
      istat = rsb_perror(c_null_ptr,istat)
      end SUBROUTINE rsb_mod_example1

      SUBROUTINE rsb_mod_example2(res)
      USE rsb
      USE iso_c_binding
      IMPLICIT NONE
      INTEGER,TARGET :: errval
      INTEGER :: res
      INTEGER :: transt = rsb_transposition_n  ! no transposition
      INTEGER :: incx = 1, incb = 1        ! X, B vectors increment
      REAL(KIND=8),TARGET :: alpha = 3,beta = 1
      INTEGER :: nnza = 4, nra = 3, nca = 3     ! nonzeroes, rows, columns of matrix A
      INTEGER,TARGET :: ia(4) = (/1, 2, 3, 3/)  ! row    indices
      INTEGER,TARGET :: ja(4) = (/1, 2, 1, 3/)  ! column indices
      INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double
      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