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 = PastixComplex64;
    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->colmajor  = 1;

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

/**
 *******************************************************************************
 *
 * @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->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;
    }

    if (spm->loc2glob != NULL) {
        for (i = 0; i < n; i++) {
            spm->loc2glob[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] ) {
        return conversionTable[spm->fmttype][ofmttype][spm->flttype]( spm );
    }
    else {
        return PASTIX_SUCCESS;
    }
}

/**
 *******************************************************************************
 *
 * @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 )
{
    double tmp;

    switch (spm->flttype) {
    case PastixFloat:
        tmp = (double)s_spmNorm( ntype, spm );
        return tmp;

    case PastixDouble:
        return d_spmNorm( ntype, spm );

    case PastixComplex32:
        tmp = (double)c_spmNorm( ntype, spm );
        return tmp;

    case PastixComplex64:
        return z_spmNorm( ntype, spm );

    default:
        return -1.;
    }
}

/**
 *******************************************************************************
 *
 * @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 )
{
    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 )
{
    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 )
{
    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 );

    /* 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->nnz     != newspm->nnz     ) )
    {
        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));

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

    if(spm->colptr != NULL) {
        newspm->colptr = (pastix_int_t*)malloc((spm->n+1) * sizeof(pastix_int_t));
        memcpy( newspm->colptr, spm->colptr, (spm->n+1) * sizeof(pastix_int_t));
    }
    if(spm->rowptr != NULL) {
        newspm->rowptr = (pastix_int_t*)malloc(spm->nnz * sizeof(pastix_int_t));
        memcpy( newspm->rowptr, spm->rowptr, spm->nnz * sizeof(pastix_int_t));
    }
    if(spm->loc2glob != NULL) {
        newspm->loc2glob = (pastix_int_t*)malloc(spm->n * sizeof(pastix_int_t));
        memcpy( newspm->loc2glob, spm->loc2glob, spm->n * sizeof(pastix_int_t));
    }
    if(spm->values != NULL) {
        size_t valsize = spm->nnz * pastix_size_of( spm->flttype );
        newspm->values = malloc(valsize);
        memcpy( newspm->values, spm->values, valsize);
    }
    return newspm;
}

/**
 *******************************************************************************
 *
 * @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.
 *
 *******************************************************************************/
/**
 * TODO: Maybe we should move down the cast of the parameters to the lowest
 * functions, and simplify this one to have identical calls to all subfunction
 */
int
spmMatVec(      int           trans,
          const void         *alpha,
          const pastix_spm_t *spm,
          const void         *x,
          const void         *beta,
                void         *y )
{
    switch (spm->mtxtype) {
    case PastixHermitian:
        switch (spm->flttype) {
        case PastixFloat:
            return s_spmSyCSCv( *((const float*)alpha), spm, (const float*)x, *((const float*)beta), (float*)y );
        case PastixComplex32:
            return c_spmHeCSCv( *((const pastix_complex32_t*)alpha), spm, (const pastix_complex32_t*)x, *((const pastix_complex32_t*)beta), (pastix_complex32_t*)y );
        case PastixComplex64:
            return z_spmHeCSCv( *((const pastix_complex64_t*)alpha), spm, (const pastix_complex64_t*)x, *((const pastix_complex64_t*)beta), (pastix_complex64_t*)y );
        case PastixDouble:
        default:
            return d_spmSyCSCv( *((const double*)alpha), spm, (const double*)x, *((const double*)beta), (double*)y );
        }
    case PastixSymmetric:
        switch (spm->flttype) {
        case PastixFloat:
            return s_spmSyCSCv( *((const float*)alpha), spm, (const float*)x, *((const float*)beta), (float*)y );
        case PastixComplex32:
            return c_spmSyCSCv( *((const pastix_complex32_t*)alpha), spm, (const pastix_complex32_t*)x, *((const pastix_complex32_t*)beta), (pastix_complex32_t*)y );
        case PastixComplex64:
            return z_spmSyCSCv( *((const pastix_complex64_t*)alpha), spm, (const pastix_complex64_t*)x, *((const pastix_complex64_t*)beta), (pastix_complex64_t*)y );
        case PastixDouble:
        default:
            return d_spmSyCSCv( *((const double*)alpha), spm, (const double*)x, *((const double*)beta), (double*)y );
        }
    case PastixGeneral:
    default:
        switch (spm->flttype) {
        case PastixFloat:
            return s_spmGeCSCv( trans, *((const float*)alpha), spm, (const float*)x, *((const float*)beta), (float*)y );
        case PastixComplex32:
            return c_spmGeCSCv( trans, *((const pastix_complex32_t*)alpha), spm, (const pastix_complex32_t*)x, *((const pastix_complex32_t*)beta), (pastix_complex32_t*)y );
        case PastixComplex64:
            return z_spmGeCSCv( trans, *((const pastix_complex64_t*)alpha), spm, (const pastix_complex64_t*)x, *((const pastix_complex64_t*)beta), (pastix_complex64_t*)y );
        case PastixDouble:
        default:
            return d_spmGeCSCv( trans, *((const double*)alpha), spm, (const double*)x, *((const double*)beta), (double*)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.
 *
 *******************************************************************************
 *
 * @return
 *      \retval PASTIX_SUCCESS if the b vector has been computed succesfully,
 *      \retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
int
spmCheckAxb( int nrhs,
             const pastix_spm_t  *spm,
                   void *x0, int ldx0,
                   void *b,  int ldb,
             const void *x,  int ldx )
{
    static int (*ptrfunc[4])(int, const pastix_spm_t *,
                             void *, int, void *, int, const void *, int) =
        {
            s_spmCheckAxb, d_spmCheckAxb, c_spmCheckAxb, z_spmCheckAxb
        };

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