spm.c

 *
 * @param[in] f
 *          File to print the spm matrix.
 *
 *******************************************************************************/
void
spmPrint( const pastix_spm_t* spm, FILE *stream )
{
    if (stream == NULL) {
        stream = stdout;
    }

    switch(spm->flttype)
    {
    case PastixPattern:
        //return p_f, spmPrint(f, spm);
        break;
    case PastixFloat:
        s_spmPrint(stream, spm);
        break;
    case PastixComplex32:
        c_spmPrint(stream, spm);
        break;
    case PastixComplex64:
        z_spmPrint(stream, spm);
        break;
    case PastixDouble:
    default:
        d_spmPrint(stream, spm);
    }
}

/**
 *******************************************************************************
 *
 * @brief Expand a multi-dof spm matrix into an spm with constant dof set 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);
    }
}

/**
 *******************************************************************************
 *
 * @brief Compute a 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:
 *          - PastixTrans
 *          - PastixNoTrans
 *          - PastixConjTrans
 *
 * @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[inout] y
 *          The vector y.
 *
 *******************************************************************************
 *
 * @retval PASTIX_SUCCESS if the y vector has been computed successfully,
 * @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 )
{
    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 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 b vector has been computed successfully,
 * @retval PASTIX_ERR_BADPARAMETER otherwise.
 *
 *******************************************************************************/
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(const pastix_complex64_t alpha, pastix_spm_t* spm)
{
    switch(spm->flttype)
    {
    case PastixPattern:
        break;
    case PastixFloat:
        s_spmScal(alpha, spm);
        break;
    case PastixComplex32:
        c_spmScal(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] alpha
 *           The scaling parameter.
 *
 * @param[in] spm
 *          The spm structure to know the type of the vector.
 *
 * @param[inout] x
 *          The vector to scal.
 *
 *******************************************************************************/
void
spmScalVector(const double alpha, pastix_spm_t* spm, void *x)
{
    switch(spm->flttype)
    {
    case PastixPattern:
        break;
    case PastixFloat:
        cblas_sscal(spm->n, alpha, x, 1);
        break;
    case PastixComplex32:
        cblas_csscal(spm->n, alpha, x, 1);
        break;
    case PastixComplex64:
        cblas_zdscal(spm->n, alpha, x, 1);
        break;
    case PastixDouble:
    default:
        cblas_dscal(spm->n, alpha, x, 1);
    }
}

/**
 * @}
 */