spm.c

           const pastix_spm_t  *spm,
           const void          *x,
           const void          *beta,
                 void          *y )
{
    pastix_spm_t *espm = (pastix_spm_t*)spm;
    int rc = PASTIX_SUCCESS;

    if ( spm->fmttype != PastixCSC ) {
        return PASTIX_ERR_BADPARAMETER;
    }

    if ( spm->dof != 1 ) {
        espm = spmExpand( spm );
    }
    switch (spm->flttype) {
    case PastixFloat:
        rc = s_spmCSCMatVec( trans, alpha, espm, x, beta, y );
        break;
    case PastixComplex32:
        rc = c_spmCSCMatVec( trans, alpha, espm, x, beta, y );
        break;
    case PastixComplex64:
        rc = z_spmCSCMatVec( trans, alpha, espm, x, beta, y );
        break;
    case PastixDouble:
    default:
        rc = d_spmCSCMatVec( trans, alpha, espm, x, beta, y );
    }

    if ( spm != espm ) {
        spmExit( espm );
        free(espm);
    }
    return rc;
}

/**
 *******************************************************************************
 *
 * @brief Compute a matrix-matrix product.
 *
 *    y = alpha * op(A) * B + beta * C
 *
 * 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:
 *          - PastixTrans
 *          - PastixNoTrans
 *          - PastixConjTrans
 *
 * @param[in] n
 *          The number of columns of the matrices B and C.
 *
 * @param[in] alpha
 *          alpha specifies the scalar alpha.
 *
 * @param[in] A
 *          The square sparse matrix A
 *
 * @param[in] B
 *          The matrix B of size ldb-by-n
 *
 * @param[in] ldb
 *          The leading dimension of the matrix B. ldb >= A->n
 *
 * @param[in] beta
 *          beta specifies the scalar beta.
 *
 * @param[inout] C
 *          The matrix C of size ldc-by-n
 *
 * @param[in] ldc
 *          The leading dimension of the matrix C. ldc >= A->n
 *
 *******************************************************************************
 *
 * @retval PASTIX_SUCCESS if the y vector has been computed successfully,
 * @retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
int
spmMatMat(       pastix_trans_t trans,
                 pastix_int_t   n,
           const void          *alpha,
           const pastix_spm_t  *A,
           const void          *B,
                 pastix_int_t   ldb,
           const void          *beta,
                 void          *C,
                 pastix_int_t   ldc )
{
    pastix_spm_t *espm = (pastix_spm_t*)A;
    int rc = PASTIX_SUCCESS;

    if ( A->fmttype != PastixCSC ) {
        return PASTIX_ERR_BADPARAMETER;
    }

    if ( A->dof != 1 ) {
        espm = spmExpand( A );
    }
    switch (A->flttype) {
    case PastixFloat:
        rc = s_spmCSCMatMat( trans, n, alpha, espm, B, ldb, beta, C, ldc );
        break;
    case PastixComplex32:
        rc = c_spmCSCMatMat( trans, n, alpha, espm, B, ldb, beta, C, ldc );
        break;
    case PastixComplex64:
        rc = z_spmCSCMatMat( trans, n, alpha, espm, B, ldb, beta, C, ldc );
        break;
    case PastixDouble:
    default:
        rc = d_spmCSCMatMat( trans, n, alpha, espm, B, ldb, beta, C, ldc );
        break;
    }

    if ( A != espm ) {
        spmExit( espm );
        free(espm);
    }
    return rc;
}

/**
 *******************************************************************************
 *
 * @brief Generate right hand side vectors associated to a given matrix.
 *
 * 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[inout] 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.
 *
 *******************************************************************************
 *
 * @retval PASTIX_SUCCESS if the b vector has been computed successfully,
 * @retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
int
spmGenRHS( pastix_rhstype_t type, pastix_int_t nrhs,
           const pastix_spm_t  *spm,
           void                *x, pastix_int_t ldx,
           void                *b, pastix_int_t ldb )
{
    static int (*ptrfunc[4])(pastix_rhstype_t, 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 );
    }
}

/**
 *******************************************************************************
 *
 * @brief 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 used to generate the right hand side, and the
 *          solution of the full problem.
 *
 * @param[inout] 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[inout] 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.
 *
 *******************************************************************************
 *
 * @retval PASTIX_SUCCESS if the tests are succesfull
 * @retval PASTIX_ERR_BADPARAMETER if the input matrix is incorrect
 * @retval 1, if one of the test failed
 *
 *******************************************************************************/
int
spmCheckAxb( pastix_int_t nrhs,
             const pastix_spm_t  *spm,
                   void *x0, pastix_int_t ldx0,
                   void *b,  pastix_int_t ldb,
             const void *x,  pastix_int_t 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 );
    }
}

/**
 *******************************************************************************
 *
 * @brief Scale the spm.
 *
 * A = alpha * A
 *
 *******************************************************************************
 *
 * @param[in] alpha
 *           The scaling parameter.
 *
 * @param[inout] spm
 *          The sparse matrix to scal.
 *
 *******************************************************************************/
void
spmScalMatrix(double alpha, pastix_spm_t* spm)
{
    switch(spm->flttype)
    {
    case PastixPattern:
        break;
    case PastixFloat:
        s_spmScal((float)alpha, spm);
        break;
    case PastixComplex32:
        c_spmScal((float)alpha, spm);
        break;
    case PastixComplex64:
        z_spmScal(alpha, spm);
        break;
    case PastixDouble:
    default:
        d_spmScal(alpha, spm);
    }
}

/**
 *******************************************************************************
 *
 * @brief Scale a vector according to the spm type.
 *
 * x = alpha * x
 *
 *******************************************************************************
 *
 * @param[in] flt
 *          Datatype of the elements in the vector that must be:
 *          @arg PastixFloat
 *          @arg PastixDouble
 *          @arg PastixComplex32
 *          @arg PastixComplex64
 *
 * @param[in] n
 *          Number of elements in the input vectors
 *
 * @param[in] alpha
 *           The scaling parameter.
 *
 * @param[inout] x
 *          The vector to scal of size ( 1 + (n-1) * abs(incx) ), and of type
 *          defined by flt.
 *
 * @param[in] incx
 *          Storage spacing between elements of x.
 *
 *******************************************************************************/
void
spmScalVector( pastix_coeftype_t flt,
               double            alpha,
               pastix_int_t      n,
               void             *x,
               pastix_int_t      incx )
{
    switch( flt )
    {
    case PastixPattern:
        break;
    case PastixFloat:
        cblas_sscal( n, (float)alpha, x, incx );
        break;
    case PastixComplex32:
        cblas_csscal( n, (float)alpha, x, incx );
        break;
    case PastixComplex64:
        cblas_zdscal( n, alpha, x, incx );
        break;
    case PastixDouble:
    default:
        cblas_dscal( n, alpha, x, incx );
    }
}

/**
 * @}
 */