laplacian.c

Go to the documentation of this file.

/**
 *
 * @file laplacian.c
 *
 * @copyright 2011-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 1.2.4
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @author Tony Delarue
 * @author Gregoire Pichon
 * @date 2024-05-29
 *
 **/
#include "common.h"
#include "spm_drivers.h"
#include "drivers/laplacian.h"
 
/**
 *
 * @ingroup spm_dev_driver
 *
 * @brief Print the usage information to generate correct Laplacian matrices.
 */
static inline void
laplacian_usage(void)
{
    fprintf(stderr,
            "Usage: genLaplacian( \"[<type>:]<dim1>[:<dim2>[:<dim3>[:<alpha>[:<beta>[:<dof>]]]]]\" )\n"
            "   Generate a Laplacian matrix M, of the form alpha * D - beta * A,\n"
            "   where D is the degree matrix, and A the adjacency matrix.\n"
            "   <type> p = pattern only\n"
            "          s = real simple\n"
            "          d = real double [default]\n"
            "          c = complex simple\n"
            "          z = complex double\n"
            "   <dim1> size of the first dimension of the laplacian\n"
            "   <dim2> size of the second dimension of the laplacian\n"
            "   <dim3> size of the third dimension of the laplacian\n"
            "   <dof>  size of the dof parameter to generate multi-dof matrices laplacian\n"
            "   Example:\n"
            "     genLaplacian( \"z:10:20\" )        generates a 2D complex double laplacian matrix of size 200.\n"
            "     genLaplacian( \"10:1:10:2.:0.5\" ) generates a 2D real double laplacian matrix of size 100 where M = 2. * D - 0.5 * A.\n"
            "     genLaplacian( \"s:10\" )           generates a 1D real single laplacian matrix of size 10.\n"
            );
}
 
/**
 *******************************************************************************
 *
 * @ingroup spm
 *
 * @brief Parse information given through the filename string to configure the
 * laplacian matrix to generate.
 *
 * The laplacian will be of size dim1 * dim2 * dim3, and will be equal to
 *     \f$ M = \alpha * D - \beta * A \f$
 *
 * where D is the degree matrix, and A the adjacency matrix.
 *
 *******************************************************************************
 *
 * @param[in] filename
 *          Configuration string of the Laplacian. See laplacian_usage() for
 *          more information.
 *
 * @param[out] flttype
 *          The floating type of the elements in the matrix.
 *
 * @param[out] dim1
 *          The first dimension of the laplacian
 *
 * @param[out] dim2
 *          The second dimension of the laplacian
 *
 * @param[out] dim3
 *          The third dimension of the laplacian
 *
 * @param[out] alpha
 *          The alpha coefficient for the degree matrix
 *
 * @param[out] beta
 *          The beta coefficient for the adjacency matrix
 *
 * @param[out] dof
 *          The degree of freedom of each unknown
 *
 *******************************************************************************
 *
 * @retval SPM_SUCCESS if the matrix has been generated successfully
 * @retval SPM_ERR_BADPARAMETER if the configuration string is incorrect
 *
 *******************************************************************************/
int
spmParseLaplacianInfo( const char     *filename,
                       spm_coeftype_t *flttype,
                       spm_int_t      *dim1,
                       spm_int_t      *dim2,
                       spm_int_t      *dim3,
                       double         *alpha,
                       double         *beta,
                       spm_int_t      *dof )
{
    double val1, val2;
    long tmp1, tmp2, tmp3, tmp4;
 
    *alpha = 1.;
    *beta  = 1.;
    *dof   = 1;
 
    /* Look for the datatype */
    {
        size_t nfilename = strlen( filename );
        char  flt, *tmpf;
        tmpf = calloc( spm_imax( nfilename, 256 ), sizeof(char) );
 
        if ( sscanf( filename, "%c:%255s", &flt, tmpf ) == 2 ) {
            filename += 2;
            switch( flt ){
            case 'Z':
            case 'z':
                *flttype = SpmComplex64;
                break;
 
            case 'C':
            case 'c':
                *flttype = SpmComplex32;
                break;
 
            case 'D':
            case 'd':
                *flttype = SpmDouble;
                break;
 
            case 'S':
            case 's':
                *flttype = SpmFloat;
                break;
 
            case 'P':
            case 'p':
                *flttype = SpmPattern;
                break;
 
            case '1':
            case '2':
            case '3':
            case '4':
            case '5':
            case '6':
            case '7':
            case '8':
            case '9':
                *flttype = SpmDouble;
                /*
                 * The first dimension is only one character long so we come
                 * back to the beginning of the string
                 */
                filename -= 2;
                break;
 
            default:
                laplacian_usage();
                free(tmpf);
                return SPM_ERR_BADPARAMETER;
            }
        }
        else {
            *flttype = SpmDouble;
        }
 
        free(tmpf);
    }
 
    /* Scan the dimensions */
    *dim1 = *dim2 = *dim3 = 1;
 
    if ( sscanf( filename, "%ld:%ld:%ld:%lf:%lf:%ld",
                 &tmp1, &tmp2, &tmp3, &val1, &val2, &tmp4 ) == 6 ) {
        *dim1 = (spm_int_t)tmp1;
        *dim2 = (spm_int_t)tmp2;
        *dim3 = (spm_int_t)tmp3;
        *alpha = val1;
        *beta  = val2;
        *dof   = (spm_int_t)tmp4;
    }
    else if ( sscanf( filename, "%ld:%ld:%ld:%lf:%lf", &tmp1, &tmp2, &tmp3, &val1, &val2 ) == 5 ) {
        *dim1 = (spm_int_t)tmp1;
        *dim2 = (spm_int_t)tmp2;
        *dim3 = (spm_int_t)tmp3;
        *alpha = val1;
        *beta  = val2;
    }
    else if ( sscanf( filename, "%ld:%ld:%ld:%lf", &tmp1, &tmp2, &tmp3, &val1 ) == 4 ) {
        *dim1 = (spm_int_t)tmp1;
        *dim2 = (spm_int_t)tmp2;
        *dim3 = (spm_int_t)tmp3;
        *alpha = val1;
    }
    else if ( sscanf( filename, "%ld:%ld:%ld", &tmp1, &tmp2, &tmp3 ) == 3 ) {
        *dim1 = (spm_int_t)tmp1;
        *dim2 = (spm_int_t)tmp2;
        *dim3 = (spm_int_t)tmp3;
    }
    else if ( sscanf( filename, "%ld:%ld", &tmp1, &tmp2 ) == 2 ) {
        *dim1 = (spm_int_t)tmp1;
        *dim2 = (spm_int_t)tmp2;
    }
    else if ( sscanf( filename, "%ld", &tmp1 ) == 1 ) {
        *dim1 = (spm_int_t)tmp1;
    }
    else {
        laplacian_usage();
        return SPM_ERR_BADPARAMETER;
    }
 
    assert( *dof != 0 );
 
    /* One of the dimension was set to 0 */
    if ( (*dim1 == 0) || (*dim2 == 0) || (*dim3 == 0) ) {
        laplacian_usage();
        return SPM_ERR_BADPARAMETER;
    }
 
    return SPM_SUCCESS;
}
 
/**
 * @ingroup spm_dev_driver
 *
 * @brief Pointers to the 7 points Laplacian generator functions
 */
static void (*laplacian_7points[6])(spmatrix_t *, spm_int_t, spm_int_t, spm_int_t, spm_fixdbl_t, spm_fixdbl_t) =
{
    p_spmLaplacian_7points,
    NULL,
    s_spmLaplacian_7points,
    d_spmLaplacian_7points,
    c_spmLaplacian_7points,
    z_spmLaplacian_7points
};
 
/**
 * @ingroup spm_dev_driver
 *
 * @brief Pointers to the 27 points Laplacian generator functions
 */
static void (*laplacian_27points[6])(spmatrix_t *, spm_int_t, spm_int_t, spm_int_t, spm_fixdbl_t, spm_fixdbl_t) =
{
    p_spmLaplacian_27points,
    NULL,
    s_spmLaplacian_27points,
    d_spmLaplacian_27points,
    c_spmLaplacian_27points,
    z_spmLaplacian_27points
};
 
/**
 *******************************************************************************
 *
 * @ingroup spm_dev_driver
 *
 * @brief Generate a Laplacian of size spm->n
 *
 *******************************************************************************
 *
 * @param[in] filename
 *          Configuration string of the Laplacian.
 *          [type:]dim1[:dim2[:dim3]]
 *             type p = pattern only\n"
 *                  s = real simple\n"
 *                  d = real double [default]\n"
 *                  c = complex simple\n"
 *                  z = complex double\n"
 *             dim1 size of the first dimension of the 1D|2D|3D laplacian\n"
 *             dim2 size of the second dimension of the 2D|3D laplacian\n"
 *             dim3 size of the third dimension of the 3D laplacian\n"
 *
 * @param[inout] spm
 *          At start, an allocated spm structure.
 *          At exit, contains a laplacian matrix in the spm format.
 *
 *******************************************************************************
 *
 * @retval SPM_SUCCESS if the matrix has been generated successfully
 * @retval SPM_ERR_BADPARAMETER if the configuration string is incorrect
 *
 *******************************************************************************/
int
genLaplacian( const char *filename,
              spmatrix_t *spm )
{
    spm_coeftype_t flttype;
    spm_int_t dim1, dim2, dim3, dof;
    double alpha = 1.;
    double beta = 1.;
    int rc;
 
    rc = spmParseLaplacianInfo( filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta, &dof );
    if (rc != SPM_SUCCESS)
        return rc;
 
    /* Build global SPM */
    spm->flttype = flttype;
    spm->gN      = dim1 * dim2 * dim3;
    spm->dof     = 1;
    laplacian_7points[spm->flttype]( spm, dim1, dim2, dim3, alpha, beta );
 
    /* Update the final values */
    spmUpdateComputedFields( spm );
 
    /* Extend the spm values if necessary */
    if ( dof != 1 ) {
        spmatrix_t spmtmp;
        int        rc;
        if ( dof < 1 ) {
            rc = spmDofExtend( spm, 1, -dof, &spmtmp );
        }
        else {
            rc = spmDofExtend( spm, 0, dof, &spmtmp );
        }
        if ( rc != SPM_SUCCESS ) {
            fprintf( stderr, "Issue while extending the matrix to multi-dof\n" );
            assert(0);
            return rc;
        }
        spmExit( spm );
        memcpy( spm, &spmtmp, sizeof(spmatrix_t) );
    }
 
    return SPM_SUCCESS;
}
 
/**
 *******************************************************************************
 *
 * @ingroup spm_dev_driver
 *
 * @brief Generate a extended Laplacian of size spm->n
 *
 *******************************************************************************
 *
 * @param[in] filename
 *          Configuration string of the Laplacian.
 *          [type:]dim1[:dim2[:dim3]]
 *             type p = pattern only
 *                    s = real simple
 *                    d = real double [default]
 *                    c = complex simple
 *                    z = complex double
 *             dim1 size of the first dimension of the 1D|2D|3D laplacian
 *             dim2 size of the second dimension of the 2D|3D laplacian
 *             dim3 size of the third dimension of the 3D laplacian
 *
 * @param[inout] spm
 *          At start, an allocated spm structure.
 *          At exit, contains a laplacian matrix in the spm format.
 *
 *******************************************************************************
 *
 * @retval SPM_SUCCESS if the matrix has been generated successfully
 * @retval SPM_ERR_BADPARAMETER if the configuration string is incorrect
 *
 *******************************************************************************/
int
genExtendedLaplacian( const char *filename,
                      spmatrix_t *spm )
{
    spm_coeftype_t flttype;
    spm_int_t dim1, dim2, dim3, dof;
    double alpha = 1.;
    double beta = 1.;
    int rc;
 
    rc = spmParseLaplacianInfo(filename, &flttype, &dim1, &dim2, &dim3, &alpha, &beta, &dof );
    if (rc != SPM_SUCCESS)
        return rc;
 
    /* Build global SPM */
    spm->flttype = flttype;
    spm->gN      = dim1 * dim2 * dim3;
    spm->dof     = 1;
    laplacian_27points[spm->flttype]( spm, dim1, dim2, dim3, alpha, beta );
 
    /* Update the info with dof if necessary */
    spmUpdateComputedFields( spm );
 
    /* Extend the spm values if necessary */
    if ( dof != 1 ) {
        spmatrix_t spmtmp;
        int        rc;
        if ( dof < 1 ) {
            rc = spmDofExtend( spm, 1, -dof, &spmtmp );
        }
        else {
            rc = spmDofExtend( spm, 0, dof, &spmtmp );
        }
        if ( rc != SPM_SUCCESS ) {
            fprintf( stderr, "Issue while extending the matrix to multi-dof\n" );
            assert(0);
            return SPM_ERR_UNKNOWN;
        }
        spmExit( spm );
        memcpy( spm, &spmtmp, sizeof(spmatrix_t) );
    }
 
    return SPM_SUCCESS;
}