spm.c

/**
 *
 * @file spm.c
 *
 *  PaStiX spm routines
 *  PaStiX is a software package provided by Inria Bordeaux - Sud-Ouest,
 *  LaBRI, University of Bordeaux 1 and IPB.
 *
 * @version 5.1.0
 * @author Xavier Lacoste
 * @author Pierre Ramet
 * @author Mathieu Faverge
 * @date 2013-06-24
 *
 **/
#include "common.h"
#include "spm.h"

#include "z_spm.h"
#include "c_spm.h"
#include "d_spm.h"
#include "s_spm.h"
#include "p_spm.h"

static int (*conversionTable[3][3][6])(pastix_spm_t*) = {
    /* From CSC */
    {{ NULL, NULL, NULL, NULL, NULL, NULL },
     { p_spmConvertCSC2CSR,
       NULL,
       s_spmConvertCSC2CSR,
       d_spmConvertCSC2CSR,
       c_spmConvertCSC2CSR,
       z_spmConvertCSC2CSR },
     { p_spmConvertCSC2IJV,
       NULL,
       s_spmConvertCSC2IJV,
       d_spmConvertCSC2IJV,
       c_spmConvertCSC2IJV,
       z_spmConvertCSC2IJV }},
    /* From CSR */
    {{ p_spmConvertCSR2CSC,
       NULL,
       s_spmConvertCSR2CSC,
       d_spmConvertCSR2CSC,
       c_spmConvertCSR2CSC,
       z_spmConvertCSR2CSC },
     { NULL, NULL, NULL, NULL, NULL, NULL },
     { p_spmConvertCSR2IJV,
       NULL,
       s_spmConvertCSR2IJV,
       d_spmConvertCSR2IJV,
       c_spmConvertCSR2IJV,
       z_spmConvertCSR2IJV }},
    /* From IJV */
    {{ p_spmConvertIJV2CSC,
       NULL,
       s_spmConvertIJV2CSC,
       d_spmConvertIJV2CSC,
       c_spmConvertIJV2CSC,
       z_spmConvertIJV2CSC },
     { p_spmConvertIJV2CSR,
       NULL,
       s_spmConvertIJV2CSR,
       d_spmConvertIJV2CSR,
       c_spmConvertIJV2CSR,
       z_spmConvertIJV2CSR },
     { NULL, NULL, NULL, NULL, NULL, NULL }}
};


/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Init the spm structure given as parameter
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          The sparse matrix to init.
 *
 *******************************************************************************/
void
spmInit( pastix_spm_t *spm )
{
    spm->mtxtype = PastixGeneral;
    spm->flttype = PastixDouble;
    spm->fmttype = PastixCSC;

    spm->gN   = 0;
    spm->n    = 0;
    spm->gnnz = 0;
    spm->nnz  = 0;

    spm->gNexp   = 0;
    spm->nexp    = 0;
    spm->gnnzexp = 0;
    spm->nnzexp  = 0;

    spm->dof      = 1;
    spm->dofs     = NULL;
    spm->layout   = PastixColMajor;

    spm->colptr   = NULL;
    spm->rowptr   = NULL;
    spm->loc2glob = NULL;
    spm->values   = NULL;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Free the spm structure
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          The sparse matrix to free.
 *
 *******************************************************************************/
void
spmExit( pastix_spm_t *spm )
{
    if(spm->colptr != NULL)
        memFree_null(spm->colptr);
    if(spm->rowptr != NULL)
        memFree_null(spm->rowptr);
    if(spm->loc2glob != NULL)
        memFree_null(spm->loc2glob);
    if(spm->values != NULL)
        memFree_null(spm->values);
    if(spm->dofs != NULL)
        memFree_null(spm->dofs);
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Rebase the spm
 *
 * Rebase the arrays of the spm to the given value. If the value is equal to the
 * original base, then nothing is performed.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          The sparse matrix to rebase.
 *
 * @param[in] baseval
 *          The base value to use in the graph (0 or 1).
 *
 *******************************************************************************/
void
spmBase( pastix_spm_t *spm,
         int           baseval )
{
    pastix_int_t baseadj;
    pastix_int_t i, n, nnz;

    /* Parameter checks */
    if ( spm == NULL ) {
        pastix_error_print("spmBase: spm pointer is NULL");
        return;
    }
    if ( (spm->colptr == NULL) ||
         (spm->rowptr == NULL) )
    {
        pastix_error_print("spmBase: spm pointer is not correctly initialized");
        return;
    }
    if ( (baseval != 0) &&
         (baseval != 1) )
    {
        pastix_error_print("spmBase: baseval is incorrect, must be 0 or 1");
        return;
    }

    baseadj = baseval - spmFindBase( spm );
    if (baseadj == 0)
	return;

    n   = spm->n;
    nnz = spm->nnz;

    switch(spm->fmttype)
    {
    case PastixCSC:
        assert( nnz == (spm->colptr[n] - spm->colptr[0]) );

        for (i = 0; i <= n; i++) {
            spm->colptr[i] += baseadj;
        }
        for (i = 0; i < nnz; i++) {
            spm->rowptr[i] += baseadj;
        }
        break;

    case PastixCSR:
        assert( nnz == (spm->rowptr[n] - spm->rowptr[0]) );
        for (i = 0; i <= n; i++) {
            spm->rowptr[i] += baseadj;
        }
        for (i = 0; i < nnz; i++) {
            spm->colptr[i] += baseadj;
        }
        break;
    case PastixIJV:
        for (i = 0; i < nnz; i++) {
            spm->rowptr[i] += baseadj;
            spm->colptr[i] += baseadj;
        }
    }

    if (spm->loc2glob != NULL) {
        for (i = 0; i < n; i++) {
            spm->loc2glob[i] += baseadj;
        }
    }
    if (spm->dofs != NULL) {
        for (i = 0; i <= n; i++) {
            spm->dofs[i] += baseadj;
        }
    }
    return;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Search the base used in the spm structure.
 *
 *******************************************************************************
 *
 * @param[in] spm
 *          The sparse matrix structure.
 *
 ********************************************************************************
 *
 * @return  The baseval used in the given sparse matrix structure.
 *
 *******************************************************************************/
pastix_int_t
spmFindBase( const pastix_spm_t *spm )
{

    pastix_int_t i, *tmp, baseval;

    /*
     * Check the baseval, we consider that arrays are sorted by columns or rows
     */
    baseval = pastix_imin( *(spm->colptr), *(spm->rowptr) );
    /*
     * if not:
     */
    if ( ( baseval != 0 ) &&
         ( baseval != 1 ) )
    {
        baseval = spm->n;
        tmp = spm->colptr;
        for(i=0; i<spm->nnz; i++, tmp++){
            baseval = pastix_imin( *tmp, baseval );
        }
    }

    return baseval;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief  Convert the storage format of the spm to ofmttype.
 *
 * Convert the storage format of the given sparse matrix from any of the
 * following format: PastixCSC, PastixCSR, or PastixIJV to one of these.
 *
 *******************************************************************************
 *
 * @param[in] ofmttype
 *          The output format of the sparse matrix. It might be PastixCSC,
 *          PastixCSR, or PastixIJV.
 *
 * @param[in,out] spm
 *          The sparse matrix structure to convert.
 *
 ********************************************************************************
 *
 * @return
 *        \retval PASTIX_SUCCESS if the conversion happened successfuly
 *        \retval PASTIX_ERR_BADPARAMETER if one the parameter is incorrect.
 *
 *******************************************************************************/
int
spmConvert( int ofmttype, pastix_spm_t *spm )
{
    if ( conversionTable[spm->fmttype][ofmttype][spm->flttype] ) {
        if ( spm->dof != 1 ) {
            //pastix_error_print( "spmConvert: Conversion of non unique dof not yet implemented\n");
            return PASTIX_ERR_NOTIMPLEMENTED;
        }
        return conversionTable[spm->fmttype][ofmttype][spm->flttype]( spm );
    }
    else {
        return PASTIX_SUCCESS;
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Convert the spm matrix into a dense matrix for test purpose. DO NOT
 * USE with large matrices.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          The sparse matrix structure to convert.
 *
 ********************************************************************************
 *
 * @return
 *        The pointer to the allocated array storing the dense version of the
 *        matrix.
 *
 *******************************************************************************/
void *
spm2Dense( const pastix_spm_t *spm )
{
    switch (spm->flttype) {
    case PastixFloat:
        return s_spm2dense( spm );
    case PastixComplex32:
        return c_spm2dense( spm );
    case PastixComplex64:
        return z_spm2dense( spm );
    case PastixDouble:
        return d_spm2dense( spm );
    default:
        return NULL;
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Compute the norm of the spm.
 *
 * Return the ntype norm of the sparse matrix spm.
 *
 *     spmNorm = ( max(abs(spm(i,j))), NORM = PastixMaxNorm
 *               (
 *               ( norm1(spm),         NORM = PastixOneNorm
 *               (
 *               ( normI(spm),         NORM = PastixInfNorm
 *               (
 *               ( normF(spm),         NORM = PastixFrobeniusNorm
 *
 *  where norm1 denotes the one norm of a matrix (maximum column sum),
 *  normI denotes the infinity norm of a matrix (maximum row sum) and
 *  normF denotes the Frobenius norm of a matrix (square root of sum
 *  of squares). Note that max(abs(spm(i,j))) is not a consistent matrix
 *  norm.
 *
 *******************************************************************************
 *
 * @param[in] ntype
 *          = PastixMaxNorm: Max norm
 *          = PastixOneNorm: One norm
 *          = PastixInfNorm: Infinity norm
 *          = PastixFrobeniusNorm: Frobenius norm
 *
 * @param[in] spm
 *          The sparse matrix structure.
 *
 ********************************************************************************
 *
 * @return
 *          \retval the norm described above. Note that for simplicity, even if
 *          the norm of single real or single complex matrix is computed with
 *          single precision, the returned norm is stored in double precision
 *          number.
 *          \retval -1., if the floating point of the sparse matrix is
 *          undefined.
 *
 *******************************************************************************/
double
spmNorm( int ntype,
         const pastix_spm_t *spm )
{
    pastix_spm_t *spmtmp = (pastix_spm_t*)spm;
    double norm = -1.;

    if ( spm->dof != 1 ) {
        fprintf(stderr, "WARNING: spm expanded due to non implemented norm for non-expanded spm\n");
        spmtmp = spmExpand( spm );
    }
    switch (spm->flttype) {
    case PastixFloat:
        norm = (double)s_spmNorm( ntype, spmtmp );
        break;

    case PastixDouble:
        norm = d_spmNorm( ntype, spmtmp );
        break;

    case PastixComplex32:
        norm = (double)c_spmNorm( ntype, spmtmp );
        break;

    case PastixComplex64:
        norm = z_spmNorm( ntype, spmtmp );
        break;

    case PastixPattern:
    default:
        ;
    }

    if ( spmtmp != spm ) {
        spmExit( spmtmp );
        free(spmtmp);
    }
    return norm;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Sort the subarray of edges of each vertex in a CSC or CSR spm
 *
 * This routine sorts the subarray of edges of each vertex in a
 * centralized spm stored in CSC or CSR format. Nothing is performed if IJV
 * format is used.
 *
 * WARNING: This function should NOT be called if dof is greater than 1.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the same sparse matrix with subarrays of edges sorted by
 *          ascending order.
 *
 ********************************************************************************
 *
 * @return
 *          \retval PASTIX_SUCCESS if the sort was called
 *          \retval PASTIX_ERR_BADPARAMETER, if the given spm was incorrect.
 *
 *******************************************************************************/
int
spmSort( pastix_spm_t *spm )
{
    if ( spm->dof != 1 ) {
        fprintf(stderr, "WARNING: spm expanded due to non implemented sort for non-expanded spm\n");
        spm = spmExpand( spm );
    }
    switch (spm->flttype) {
    case PastixPattern:
        p_spmSort( spm );
        break;
    case PastixFloat:
        s_spmSort( spm );
        break;
    case PastixDouble:
        d_spmSort( spm );
        break;
    case PastixComplex32:
        c_spmSort( spm );
        break;
    case PastixComplex64:
        z_spmSort( spm );
        break;
    default:
        return PASTIX_ERR_BADPARAMETER;
    }
    return PASTIX_SUCCESS;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Merge mulitple entries in a spm by summing them.
 *
 * This routine merge the multiple entries in a sparse
 * matrix by suming their values together. The sparse matrix needs to be sorted
 * first (see spmSort()).
 *
 * WARNING: Not implemented for CSR and IJV format.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the reduced sparse matrix of multiple entries were present
 *          in it. The multiple values for a same vertex are sum up together.
 *
 ********************************************************************************
 *
 * @return
 *          \retval If >=0, the number of vertices that were merged
 *          \retval PASTIX_ERR_BADPARAMETER, if the given spm was incorrect.
 *
 *******************************************************************************/
pastix_int_t
spmMergeDuplicate( pastix_spm_t *spm )
{
    if ( spm->dof != 1 ) {
        fprintf(stderr, "WARNING: spm expanded due to non implemented merge for non-expanded spm\n");
        spm = spmExpand( spm );
    }
    switch (spm->flttype) {
    case PastixPattern:
        return p_spmMergeDuplicate( spm );

    case PastixFloat:
        return s_spmMergeDuplicate( spm );

    case PastixDouble:
        return d_spmMergeDuplicate( spm );

    case PastixComplex32:
        return c_spmMergeDuplicate( spm );

    case PastixComplex64:
        return z_spmMergeDuplicate( spm );

    default:
        return PASTIX_ERR_BADPARAMETER;
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Symmetrize the pattern of the spm
 *
 * Symmetrize the pattern of the input spm by edges when A(i,j) exists, but
 * A(j,i) does not. When values are associated to the edge, zeroes are added to
 * the values array.
 *
 * WARNING: Not implemented for CSR and IJV format.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the reduced sparse matrix of multiple entries were present
 *          in it. The multiple values for a same vertex are sum up together.
 *
 ********************************************************************************
 *
 * @return
 *          \retval If >=0, the number of vertices that were merged
 *          \retval PASTIX_ERR_BADPARAMETER, if the given spm was incorrect.
 *
 *******************************************************************************/
pastix_int_t
spmSymmetrize( pastix_spm_t *spm )
{
    if ( spm->dof != 1 ) {
        fprintf(stderr, "WARNING: spm expanded due to non implemented symmetrize for non-expanded spm\n");
        spm = spmExpand( spm );
    }
    switch (spm->flttype) {
    case PastixPattern:
        return p_spmSymmetrize( spm );

    case PastixFloat:
        return s_spmSymmetrize( spm );

    case PastixDouble:
        return d_spmSymmetrize( spm );

    case PastixComplex32:
        return c_spmSymmetrize( spm );

    case PastixComplex64:
        return z_spmSymmetrize( spm );

    default:
        return PASTIX_ERR_BADPARAMETER;
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Check the correctness of a spm.
 *
 * This routine initializes the sparse matrix to fit the PaStiX requirements. If
 * needed, the format is changed to CSC, the duplicated vertices are merged
 * together by summing their values; the graph is made symmetric for matrices
 * with unsymmetric pattern, new values are set to 0.; Only the lower part is
 * kept for the symmetric matrices.
 *
 * On exit, if no changes have been made, the initial sparse matrix is returned,
 * otherwise a copy of the sparse matrix structured fixed to meet the PaStiX
 * requirements is returned.
 *
 *******************************************************************************
 *
 * @param[in,out] spm
 *          The pointer to the sparse matrix structure to check, and correct.
 *          On exit, the subarrays related to each column might have been sorted
 *          by ascending order.
 *
 *******************************************************************************
 *
 * @return
 *          \retval If no modifications were made to the initial matrix
 *                  structure, the one given as parameter is returned
 *          \retval Otherwise, the news sparse matrix structure is returned. It
 *                  must be destroyed with spmExit() and a free of the returned
 *                  pointer.
 *
 *******************************************************************************/
pastix_spm_t *
spmCheckAndCorrect( pastix_spm_t *spm )
{
    pastix_spm_t *newspm = NULL;
    pastix_int_t count;

    /* Let's work on a copy */
    newspm = spmCopy( spm );

    /* PaStiX works on CSC matrices */
    spmConvert( PastixCSC, newspm );

    if ( newspm->dof != 1 ) {
        fprintf(stderr, "WARNING: newspm expanded due to missing check functions implementations\n");
        newspm = spmExpand( newspm );
    }

    /* Sort the rowptr for each column */
    spmSort( newspm );

    /* Merge the duplicated entries by summing the values */
    count = spmMergeDuplicate( newspm );
    if ( count > 0 ) {
        fprintf(stderr, "spmCheckAndCorrect: %ld entries have been merged\n", (int64_t)count );
    }

    /**
     * If the matrix is symmetric or hermitian, we keep only the upper or lower
     * part, otherwise, we symmetrize the graph to get A+A^t, new values are set
     * to 0.
     */
    if ( newspm->mtxtype == PastixGeneral ) {
        count = spmSymmetrize( newspm );
        if ( count > 0 ) {
            fprintf(stderr, "spmCheckAndCorrect: %ld entries have been added for symmetry\n", (int64_t)count );
        }
    }
    else {
        //spmToLower( newspm );
    }

    /**
     * Check if we return the new one, or the original one because no changes
     * have been made
     */
    if (( spm->fmttype != newspm->fmttype ) ||
        ( spm->nnzexp  != newspm->nnzexp  ) )
    {
        return newspm;
    }
    else {
        spmExit( newspm );
        free(newspm);
        return (pastix_spm_t*)spm;
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Create a copy of the spm
 *
 * Duplicate the spm data structure given as parameter. All new arrays are
 * allocated and copied from the original matrix. Both matrices need to be
 * freed.
 *
 *******************************************************************************
 *
 * @param[in] spm
 *          The sparse matrix to copy.
 *
 *******************************************************************************
 *
 * @return
 *          The copy of the sparse matrix.
 *
 *******************************************************************************/
pastix_spm_t *
spmCopy( const pastix_spm_t *spm )
{
    pastix_spm_t *newspm = (pastix_spm_t*)malloc(sizeof(pastix_spm_t));
    pastix_int_t colsize, rowsize, valsize, dofsize;

    memcpy( newspm, spm, sizeof(pastix_spm_t));

    switch(spm->fmttype){
    case PastixCSC:
        colsize = spm->n + 1;
        rowsize = spm->nnz;
        valsize = spm->nnzexp;
        dofsize = spm->n + 1;
        break;
    case PastixCSR:
        colsize = spm->nnz;
        rowsize = spm->n + 1;
        valsize = spm->nnzexp;
        dofsize = spm->n + 1;
        break;
    case PastixIJV:
        colsize = spm->nnz;
        rowsize = spm->nnz;
        valsize = spm->nnzexp;
        dofsize = spm->n + 1;
    }

    if(spm->colptr != NULL) {
        newspm->colptr = (pastix_int_t*)malloc( colsize * sizeof(pastix_int_t));
        memcpy( newspm->colptr, spm->colptr, colsize * sizeof(pastix_int_t));
    }
    if(spm->rowptr != NULL) {
        newspm->rowptr = (pastix_int_t*)malloc(rowsize * sizeof(pastix_int_t));
        memcpy( newspm->rowptr, spm->rowptr, rowsize * sizeof(pastix_int_t));
    }
    if(spm->loc2glob != NULL) {
        newspm->loc2glob = (pastix_int_t*)malloc(dofsize * sizeof(pastix_int_t));
        memcpy( newspm->loc2glob, spm->loc2glob, dofsize * sizeof(pastix_int_t));
    }
    if(spm->dofs != NULL) {
        newspm->dofs = (pastix_int_t*)malloc(dofsize * sizeof(pastix_int_t));
        memcpy( newspm->dofs, spm->dofs, dofsize * sizeof(pastix_int_t));
    }
    if(spm->values != NULL) {
        valsize = valsize * pastix_size_of( spm->flttype );
        newspm->values = malloc(valsize);
        memcpy( newspm->values, spm->values, valsize);
    }
    return newspm;
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Expand a multi-dof spm matrix into an spm with constant dof to 1.
 *
 * Duplicate the spm data structure given as parameter. All new arrays are
 * allocated and copied from the original matrix. Both matrices need to be
 * freed.
 *
 *******************************************************************************
 *
 * @param[in] spm
 *          The sparse matrix to copy.
 *
 *******************************************************************************
 *
 * @return
 *          The copy of the sparse matrix.
 *
 *******************************************************************************/
pastix_spm_t *
spmExpand(const pastix_spm_t* spm)
{
    switch(spm->flttype)
    {
    case PastixPattern:
        return p_spmExpand(spm);
        break;
    case PastixFloat:
        return s_spmExpand(spm);
        break;
    case PastixComplex32:
        return c_spmExpand(spm);
        break;
    case PastixComplex64:
        return z_spmExpand(spm);
        break;
    case PastixDouble:
    default:
        return d_spmExpand(spm);
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Compute a matrix-vector product
 *
 * Compute the matrix vector product:
 *
 *    y = alpha * op(A) * x + beta * y, where op(A) is one of
 *
 *    op( A ) = A  or op( A ) = A' or op( A ) = conjg( A' )
 *
 *  alpha and beta are scalars, and x and y are vectors.
 *
 *******************************************************************************
 *
 * @param[in] trans
 *          Specifies whether the matrix spm is transposed, not transposed or
 *          conjugate transposed:
 *          = PastixNoTrans:   A is not transposed;
 *          = PastixTrans:     A is transposed;
 *          = PastixConjTrans: A is conjugate transposed.
 *
 * @param[in] alpha
 *          alpha specifies the scalar alpha.
 *
 * @param[in] spm
 *          The PastixGeneral spm.
 *
 * @param[in] x
 *          The vector x.
 *
 * @param[in] beta
 *          beta specifies the scalar beta.
 *
 * @param[in,out] y
 *          The vector y.
 *
 *******************************************************************************
 *
 * @return
 *      \retval PASTIX_SUCCESS if the y vector has been computed succesfully,
 *      \retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
int
spmMatVec(const pastix_trans_t trans,
          const void          *alpha,
          const pastix_spm_t  *spm,
          const void          *x,
          const void          *beta,
                void          *y )
{
    if ( spm->fmttype != PastixCSC ) {
        return PASTIX_ERR_BADPARAMETER;
    }

    switch (spm->flttype) {
    case PastixFloat:
        return s_spmCSCMatVec( trans, alpha, spm, x, beta, y );
    case PastixComplex32:
        return c_spmCSCMatVec( trans, alpha, spm, x, beta, y );
    case PastixComplex64:
        return z_spmCSCMatVec( trans, alpha, spm, x, beta, y );
    case PastixDouble:
    default:
        return d_spmCSCMatVec( trans, alpha, spm, x, beta, y );
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Generate right hand sides.
 *
 * Generate nrhs right hand side vectors associated to a given matrix to test a
 * problem with a solver. The vectors can be initialized randomly, or to get a
 * specific solution.
 *
 *******************************************************************************
 *
 * @param[in] type
 *          Defines how to compute the vector b.
 *          - PastixRhsOne:  b is computed such that x = 1 [ + I ]
 *          - PastixRhsI:    b is computed such that x = i [ + i * I ]
 *          - PastixRhsRndX: b is computed by matrix-vector product, such that
 *            is a random vector in the range [-0.5, 0.5]
 *          - PastixRhsRndB: b is computed randomly and x is not computed.
 *
 * @param[in] nrhs
 *          Defines the number of right hand side that must be generated.
 *
 * @param[in] spm
 *          The sparse matrix uses to generate the right hand side, and the
 *          solution of the full problem.
 *
 * @param[out] x
 *          On exit, if x != NULL, then the x vector(s) generated to compute b
 *          is returned. Must be of size at least ldx * spm->n.
 *
 * @param[in] ldx
 *          Defines the leading dimension of x when multiple right hand sides
 *          are available. ldx >= spm->n.
 *
 * @param[in,out] b
 *          b must be an allocated matrix of size at least ldb * nrhs.
 *          On exit, b is initialized as defined by the type parameter.
 *
 * @param[in] ldb
 *          Defines the leading dimension of b when multiple right hand sides
 *          are available. ldb >= spm->n.
 *
 *******************************************************************************
 *
 * @return
 *      \retval PASTIX_SUCCESS if the b vector has been computed succesfully,
 *      \retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
int
spmGenRHS( int type, int nrhs,
           const pastix_spm_t  *spm,
           void                *x, int ldx,
           void                *b, int ldb )
{
    static int (*ptrfunc[4])(int, int,
                             const pastix_spm_t *,
                             void *, int, void *, int) =
        {
            s_spmGenRHS, d_spmGenRHS, c_spmGenRHS, z_spmGenRHS
        };

    int id = spm->flttype - PastixFloat;
    if ( (id < 0) || (id > 3) ) {
        return PASTIX_ERR_BADPARAMETER;
    }
    else {
        return ptrfunc[id](type, nrhs, spm, x, ldx, b, ldb );
    }
}

/**
 *******************************************************************************
 *
 * @ingroup pastix_spm
 *
 * @brief Check backward and forward errors
 *
 * Check the backward error, and the forward error if x0 is provided.
 *
 *******************************************************************************
 *
 * @param[in] nrhs
 *          Defines the number of right hand side that must be generated.
 *
 * @param[in] spm
 *          The sparse matrix uses to generate the right hand side, and the
 *          solution of the full problem.
 *
 * @param[in,out] x0
 *          If x0 != NULL, the forward error is computed.
 *          On exit, x0 stores (x0-x)
 *
 * @param[in] ldx0
 *          Defines the leading dimension of x0 when multiple right hand sides
 *          are available. ldx0 >= spm->n.
 *
 * @param[in,out] b
 *          b is a matrix of size at least ldb * nrhs.
 *          On exit, b stores Ax-b.
 *
 * @param[in] ldb
 *          Defines the leading dimension of b when multiple right hand sides
 *          are available. ldb >= spm->n.
 *
 * @param[in] x
 *          Contains the solution computed by the solver.
 *
 * @param[in] ldx
 *          Defines the leading dimension of x when multiple right hand sides
 *          are available. ldx >= spm->n.
 *
 *******************************************************************************
 *