starpu_dpotrf.c

Go to the documentation of this file.

/**
 *
 * @file starpu_dpotrf.c
 *
 * PaStiX dpotrf StarPU wrapper.
 *
 * @copyright 2016-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 6.3.2
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @date 2023-07-21
 * @generated from /builds/solverstack/pastix/sopalin/starpu/starpu_zpotrf.c, normal z -> d, Wed Dec 13 12:09:29 2023
 *
 * @addtogroup starpu_potrf
 * @{
 *
 **/
#include "common.h"
#include "blend/solver.h"
#include "sopalin/sopalin_data.h"
#include "pastix_dcores.h"
#include "pastix_starpu.h"
#include "pastix_dstarpu.h"
 
/**
 *******************************************************************************
 *
 * @brief TODO
 *
 *******************************************************************************
 *
 * @param[in] sopalin_data
 *          TODO
 *
 * @param[in] cblk
 *          TODO
 *
 * @param[in] fcblk
 *          TODO
 *
 * @param[in] blokA
 *          TODO
 *
 * @param[in] prio
 *          TODO
 *
 *******************************************************************************/
void
starpu_task_potrf_dgemmsp( sopalin_data_t   *sopalin_data,
                           SolverCblk       *cblk,
                           const SolverBlok *blokB,
                           SolverCblk       *fcblk,
                           int               prio )
{
    const SolverBlok *blokA, *lblk;
 
    if ( cblk->cblktype & CBLK_TASKS_2D ) {
        lblk = cblk[1].fblokptr;
        for ( blokA = blokB; blokA < lblk; blokA++ ) {
 
            starpu_task_blok_dgemmsp( sopalin_data, PastixLCoef, PastixLCoef, PastixTrans,
                                      cblk, fcblk, blokA, blokB, prio );
 
            /* Skip B blocks facing the same cblk */
            while ( ( blokA < lblk ) &&
                    ( blokA[0].fcblknm == blokA[1].fcblknm ) &&
                    ( blokA[0].lcblknm == blokA[1].lcblknm ) )
            {
                blokA++;
            }
        }
    }
    else {
        starpu_task_cblk_dgemmsp( sopalin_data, PastixLCoef, PastixLCoef, PastixTrans,
                                  cblk, blokB, fcblk, prio );
    }
}
 
/**
 *******************************************************************************
 *
 * @brief TODO
 *
 *******************************************************************************
 *
 * @param[in] sopalin_data
 *          TODO
 *
 * @param[in] cblk
 *          TODO
 *
 * @param[in] prio
 *          TODO
 *
 *******************************************************************************/
void
starpu_task_dpotrfsp( sopalin_data_t *sopalin_data,
                      SolverCblk     *cblk,
                      int             prio )
{
    SolverBlok   *lblk, *blok;
    pastix_int_t  m;
 
    if ( cblk->cblktype & CBLK_TASKS_2D ) {
        starpu_task_blok_dpotrf( sopalin_data, cblk, prio );
 
        lblk = cblk[1].fblokptr;
        for ( blok = cblk->fblokptr + 1, m = 0; blok < lblk; blok++, m++ ) {
 
            starpu_task_blok_dtrsmsp( sopalin_data, PastixLCoef, PastixRight, PastixLower,
                                      PastixTrans, PastixNonUnit,
                                      cblk, blok, prio );
 
            /* Skip blocks facing the same cblk */
            while ( ( blok < lblk ) &&
                    ( blok[0].fcblknm == blok[1].fcblknm ) &&
                    ( blok[0].lcblknm == blok[1].lcblknm ) )
            {
                blok++;
            }
        }
    }
    else {
        starpu_task_cblk_dpotrfsp( sopalin_data, cblk, prio );
    }
}
 
/**
 *******************************************************************************
 *
 * @brief Perform a sparse Cholesky factorization with 1D kernels.
 *
 * The function performs the Cholesky factorization of a sparse symmetric
 * positive definite (or Symmetric positive definite in the real case) matrix
 * A.
 * The factorization has the form
 *
 *    \f[ A = L\times L^T \f]
 *
 * where L is a sparse lower triangular matrix.
 *
 *******************************************************************************
 *
 * @param[inout] sopalin_data
 *          Solver matrix information structure that will guide the algorithm.
 *
 * @param[inout] desc
 *          StarPU descriptor of the sparse matrix.
 *
 ******************************************************************************/
void
starpu_dpotrf_sp1dplus_rl( sopalin_data_t              *sopalin_data,
                           starpu_sparse_matrix_desc_t *desc )
{
    const SolverMatrix *solvmtx = sopalin_data->solvmtx;
    SolverCblk         *cblk, *fcblk;
    SolverBlok         *blok, *lblk;
    pastix_int_t        k, m, cblknbr, cblk_n;
 
    cblknbr = solvmtx->cblknbr;
    cblk = solvmtx->cblktab;
    for (k=0; k<solvmtx->cblknbr; k++, cblk++){
 
        if ( cblk->cblktype & CBLK_IN_SCHUR ) {
            break;
        }
 
        starpu_task_cblk_dpotrfsp( sopalin_data, cblk,
                                   cblknbr - k );
 
        blok = cblk->fblokptr + 1; /* this diagonal block     */
        lblk = cblk[1].fblokptr;   /* the next diagonal block */
 
        /* if there are off-diagonal supernodes in the column */
        for(m=0; blok < lblk; blok++, m++ )
        {
            fcblk = (solvmtx->cblktab + blok->fcblknm);
            cblk_n = fcblk - solvmtx->cblktab;
            starpu_task_cblk_dgemmsp( sopalin_data, PastixLCoef, PastixLCoef, PastixTrans,
                                      cblk, blok, fcblk,
                                      cblknbr - pastix_imin( k + m, cblk_n ) );
        }
        starpu_sparse_cblk_wont_use( PastixLCoef, cblk );
    }
    (void)desc;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform a sparse Cholesky factorization with 1D kernels.
 *
 * The function performs the Cholesky factorization of a sparse symmetric
 * positive definite (or Symmetric positive definite in the real case) matrix
 * A.
 * The factorization has the form
 *
 *    \f[ A = L\times L^T \f]
 *
 * where L is a sparse lower triangular matrix.
 *
 *******************************************************************************
 *
 * @param[inout] sopalin_data
 *          Solver matrix information structure that will guide the algorithm.
 *
 * @param[inout] desc
 *          StarPU descriptor of the sparse matrix.
 *
 ******************************************************************************/
void
starpu_dpotrf_sp1dplus_ll( sopalin_data_t              *sopalin_data,
                           starpu_sparse_matrix_desc_t *desc )
{
    const SolverMatrix *solvmtx = sopalin_data->solvmtx;
    SolverCblk         *cblk, *fcblk;
    SolverBlok         *blok;
    pastix_int_t        k, m, cblknbr;
 
    cblknbr = solvmtx->cblknbr;
    fcblk   = solvmtx->cblktab;
    for ( k = 0; k < solvmtx->cblknbr; k++, fcblk++ ) {
 
        for ( m = fcblk[0].brownum; m < fcblk[1].brownum; m++ ) {
            blok = solvmtx->bloktab + solvmtx->browtab[m];
            cblk = solvmtx->cblktab + blok->lcblknm;
 
            if ( cblk->cblktype & CBLK_IN_SCHUR ) {
                break;
            }
 
            starpu_task_cblk_dgemmsp( sopalin_data, PastixLCoef, PastixLCoef, PastixTrans,
                                      cblk, blok, fcblk, cblknbr - k );
        }
 
        if ( fcblk->cblktype & CBLK_IN_SCHUR ) {
            continue;
        }
 
        starpu_task_cblk_dpotrfsp( sopalin_data, fcblk,
                                   cblknbr - k );
    }
 
    cblk = solvmtx->cblktab;
    for ( k = 0; k < solvmtx->cblknbr; k++, cblk++ ) {
        starpu_sparse_cblk_wont_use( PastixLCoef, cblk );
    }
    (void)desc;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform a sparse Cholesky factorization with 1D and 2D kernels.
 *
 * The function performs the Cholesky factorization of a sparse symmetric
 * positive definite (or Symmetric positive definite in the real case) matrix
 * A.
 * The factorization has the form
 *
 *    \f[ A = L\times L^T \f]
 *
 * where L is a sparse lower triangular matrix.
 *
 *******************************************************************************
 *
 * @param[inout] sopalin_data
 *          Solver matrix information structure that will guide the algorithm.
 *
 * @param[inout] desc
 *          StarPU descriptor of the sparse matrix.
 *
 ******************************************************************************/
void
starpu_dpotrf_sp2d_rl( sopalin_data_t              *sopalin_data,
                       starpu_sparse_matrix_desc_t *desc )
{
    const SolverMatrix *solvmtx = sopalin_data->solvmtx;
    SolverCblk         *cblk, *fcblk;
    SolverBlok         *blok, *lblk;
    pastix_int_t        k, m, cblknbr, cblk_n;
 
    cblknbr = solvmtx->cblknbr;
 
    /* Let's submit all 1D tasks first */
    cblk = solvmtx->cblktab;
    for ( k = 0; k <= solvmtx->cblkmax1d; k++, cblk++ ) {
 
        if ( cblk->cblktype & CBLK_IN_SCHUR ) {
            break;
        }
 
        if ( cblk->cblktype & CBLK_TASKS_2D ) {
            continue;
        }
 
        starpu_task_dpotrfsp( sopalin_data, cblk, cblknbr - k );
 
        lblk = cblk[1].fblokptr; /* the next diagonal block */
 
        for ( blok = cblk->fblokptr + 1, m = 0; blok < lblk; blok++, m++ ) {
            fcblk  = ( solvmtx->cblktab + blok->fcblknm );
            cblk_n = ( cblk->cblktype & CBLK_TASKS_2D ) ? blok->fcblknm : fcblk - solvmtx->cblktab;
 
            starpu_task_potrf_dgemmsp( sopalin_data, cblk, blok, fcblk,
                                       cblknbr - pastix_imin( k + m, cblk_n ) );
        }
        starpu_sparse_cblk_wont_use( PastixLCoef, cblk );
    }
 
    /* Let's submit all 2D tasks */
    cblk = solvmtx->cblktab + solvmtx->cblkmin2d;
    for ( k = solvmtx->cblkmin2d; k < solvmtx->cblknbr; k++, cblk++ ) {
 
        if ( cblk->cblktype & CBLK_IN_SCHUR ) {
            continue;
        }
 
        if ( ! ( cblk->cblktype & CBLK_TASKS_2D ) ) {
            continue;
        }
 
        starpu_task_dpotrfsp( sopalin_data, cblk, cblknbr - k );
 
        lblk = cblk[1].fblokptr; /* the next diagonal block */
 
        for ( blok = cblk->fblokptr + 1, m = 0; blok < lblk; blok++, m++ ) {
            fcblk  = ( solvmtx->cblktab + blok->fcblknm );
            cblk_n = ( cblk->cblktype & CBLK_TASKS_2D ) ? blok->fcblknm : fcblk - solvmtx->cblktab;
 
            starpu_task_potrf_dgemmsp( sopalin_data, cblk, blok, fcblk,
                                       cblknbr - pastix_imin( k + m, cblk_n ) );
 
            /* Skip A blocks facing the same cblk */
            while( ( blok < lblk ) &&
                    ( blok[0].fcblknm == blok[1].fcblknm ) &&
                    ( blok[0].lcblknm == blok[1].lcblknm ) )
            {
                blok++;
            }
        }
        starpu_sparse_cblk_wont_use( PastixLCoef, cblk );
    }
 
    (void)desc;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform a sparse Cholesky factorization with 1D and 2D kernels.
 *
 * The function performs the Cholesky factorization of a sparse symmetric
 * positive definite (or Symmetric positive definite in the real case) matrix
 * A.
 * The factorization has the form
 *
 *    \f[ A = L\times L^T \f]
 *
 * where L is a sparse lower triangular matrix.
 *
 *******************************************************************************
 *
 * @param[inout] sopalin_data
 *          Solver matrix information structure that will guide the algorithm.
 *
 * @param[inout] desc
 *          StarPU descriptor of the sparse matrix.
 *
 ******************************************************************************/
void
starpu_dpotrf_sp2d_ll( sopalin_data_t              *sopalin_data,
                       starpu_sparse_matrix_desc_t *desc )
{
    const SolverMatrix *solvmtx = sopalin_data->solvmtx;
    SolverCblk         *cblk, *fcblk;
    SolverBlok         *blok = NULL;
    SolverBlok         *blok_prev;
    pastix_int_t        k, m, cblknbr;
 
    cblknbr = solvmtx->cblknbr;
    fcblk   = solvmtx->cblktab;
    for ( k = 0; k < cblknbr; k++, fcblk++ ) {
 
        for ( m = fcblk[0].brownum; m < fcblk[1].brownum; m++ ) {
            blok_prev = blok;
            blok = solvmtx->bloktab + solvmtx->browtab[m];
            cblk = solvmtx->cblktab + blok->lcblknm;
 
            if ( cblk->cblktype & CBLK_IN_SCHUR ) {
                continue;
            }
 
            if( ( cblk->cblktype & CBLK_TASKS_2D )      &&
                ( blok_prev          != NULL          ) &&
                ( blok_prev->fcblknm == blok->fcblknm ) &&
                ( blok_prev->lcblknm == blok->lcblknm ) )
            {
                continue;
            }
 
            starpu_task_potrf_dgemmsp( sopalin_data, cblk, blok, fcblk,
                                       cblknbr - k );
        }
 
        if ( fcblk->cblktype & CBLK_IN_SCHUR ) {
            continue;
        }
 
       starpu_task_dpotrfsp( sopalin_data, fcblk, cblknbr - k );
    }
 
    cblk = solvmtx->cblktab;
    for ( k = 0; k < solvmtx->cblknbr; k++, cblk++ ) {
        starpu_sparse_cblk_wont_use( PastixLCoef, cblk );
    }
 
    (void)desc;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform a sparse Cholesky factorization using StarPU runtime.
 *
 * The function performs the Cholesky factorization of a sparse symmetric
 * positive definite (or Symmetric positive definite in the real case) matrix
 * A.
 * The factorization has the form
 *
 *    \f[ A = L\times L^T \f]
 *
 * where L is a sparse lower triangular matrix.
 *
 * The algorithm is automatically chosen between the 1D and 2D version based on
 * the API parameter IPARM_TASKS2D_LEVEL. If IPARM_TASKS2D_LEVEL != 0
 * the 2D scheme is applied, the 1D otherwise.
 *
 *******************************************************************************
 *
 * @param[inout] pastix_data
 *          The pastix_data structure that describes the solver instance.
 *
 * @param[inout] sopalin_data
 *          Solver matrix information structure that will guide the algorithm.
 *
 ******************************************************************************/
void
starpu_dpotrf( pastix_data_t  *pastix_data,
               sopalin_data_t *sopalin_data )
{
    starpu_sparse_matrix_desc_t *sdesc = sopalin_data->solvmtx->starpu_desc;
    double sub = 0.;
    double com = 0.;
 
    /*
     * Start StarPU if not already started
     */
    if (pastix_data->starpu == NULL) {
        int argc = 0;
        pastix_starpu_init( pastix_data, &argc, NULL, NULL );
    }
 
    if ( sdesc == NULL ) {
        /* Create the matrix descriptor */
        starpu_sparse_matrix_init( sopalin_data->solvmtx,
                                   PastixSymmetric,
                                   pastix_data->inter_node_procnbr,
                                   pastix_data->inter_node_procnum,
                                   PastixDouble );
        sdesc = sopalin_data->solvmtx->starpu_desc;
    }
 
    starpu_profiling_status_set(STARPU_PROFILING_ENABLE);
#if defined(STARPU_USE_FXT)
    if (pastix_data->iparm[IPARM_TRACE] & PastixTraceNumfact) {
        starpu_fxt_start_profiling();
    }
#endif
#if defined(PASTIX_STARPU_STATS)
    clockStart( sub );
#else
    starpu_resume();
#endif
    /*
     * Select 1D or 2D algorithm based on 2d tasks level
     */
    if ( pastix_data->iparm[IPARM_TASKS2D_LEVEL] != 0 )
    {
        if ( pastix_data->iparm[IPARM_FACTO_LOOK_SIDE] == PastixFactLeftLooking ) {
            starpu_dpotrf_sp2d_ll( sopalin_data, sdesc );
        }
        else {
            starpu_dpotrf_sp2d_rl( sopalin_data, sdesc );
        }
    }
    else
    {
        if ( pastix_data->iparm[IPARM_FACTO_LOOK_SIDE] == PastixFactLeftLooking ) {
            starpu_dpotrf_sp1dplus_ll( sopalin_data, sdesc );
        }
        else {
            starpu_dpotrf_sp1dplus_rl( sopalin_data, sdesc );
        }
    }
 
    starpu_sparse_matrix_getoncpu( sdesc );
#if defined(PASTIX_STARPU_STATS)
    clockStop( sub );
    clockStart( com );
    starpu_resume();
#endif
    starpu_task_wait_for_all();
#if defined(PASTIX_WITH_MPI)
    starpu_mpi_wait_for_all( pastix_data->pastix_comm );
    starpu_mpi_barrier( pastix_data->inter_node_comm );
#endif
    starpu_pause();
#if defined(STARPU_USE_FXT)
    if (pastix_data->iparm[IPARM_TRACE] & PastixTraceNumfact) {
        starpu_fxt_stop_profiling();
    }
#endif
    starpu_profiling_status_set(STARPU_PROFILING_DISABLE);
#if defined(PASTIX_STARPU_STATS)
    clockStop( com );
    print_stats( sub, com, pastix_data->solvmatr );
#endif
 
    (void)com;
    (void)sub;
    return;
}
 
/**
 *@}
 */