core_slr2xx.c

Go to the documentation of this file.

/**
 *
 * @file core_slr2xx.c
 *
 * PaStiX low-rank kernel routines that perform the addition of AB into C.
 *
 * @copyright 2016-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 6.4.0
 * @author Mathieu Faverge
 * @author Gregoire Pichon
 * @author Pierre Ramet
 * @date 2024-07-05
 * @generated from /builds/solverstack/pastix/kernels/core_zlr2xx.c, normal z -> s, Thu Aug 29 14:20:19 2024
 *
 **/
#include "common.h"
#include <cblas.h>
#include "kernels_trace.h"
#include "blend/solver.h"
#include "pastix_scores.h"
#include "pastix_slrcores.h"
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
static float sone  =  1.0;
static float szero =  0.0;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
 
/**
 *******************************************************************************
 *
 * @brief Perform the addition of the low-rank matrix AB and the full-rank
 * matrix C.
 *
 *******************************************************************************
 *
 * @param[inout] params
 *          The LRMM structure that stores all the parameters used in the LRMM
 *          functions family.
 *          On exit, the C matrix is udpated with the addition of AB.
 *          @sa core_slrmm_t
 *
 * @param[in] AB
 *          The low-rank structure of the AB matrix to apply to C.
 *
 * @param[in] transV
 *          Specify if AB->v is stored normally or transposed.
 *          - If PastixNoTrans, AB->v is stored normally for low-rank format.
 *          - If PastixTrans, AB->v is stored transposed.
 *          - If PastixTrans, AB->v is stored transposed, and () must be
 *          applied to the matrix.
 *
 *******************************************************************************
 *
 * @return The number of flops required to perform the operation.
 *
 *******************************************************************************/
static inline pastix_fixdbl_t
core_slr2fr( core_slrmm_t           *params,
             const pastix_lrblock_t *AB,
             pastix_trans_t          transV )
{
    PASTE_CORE_SLRMM_PARAMS( params );
    pastix_int_t ldabu = M;
    pastix_int_t ldabv = (transV == PastixNoTrans) ? AB->rkmax : N;
    pastix_fixdbl_t flops = 0.;
    float *Cfr = C->u;
    Cfr += Cm * offy + offx;
 
    assert( C->rk == -1 );
 
    /* TODO: find a suitable name to trace this kind of kernel. */
    if ( AB->rk == -1 ) {
        flops = 2 * M * N;
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_updateCfr );
        core_sgeadd( PastixNoTrans, M, N,
                     alpha, AB->u, ldabu,
                     beta,  Cfr,   Cm );
        kernel_trace_stop_lvl2( flops );
    }
    else {
        flops = FLOPS_SGEMM( M, N, AB->rk );
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_updateCfr );
        cblas_sgemm( CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)transV,
                     M, N, AB->rk,
                     (alpha), AB->u, ldabu,
                                         AB->v, ldabv,
                     (beta),  Cfr,   Cm );
        kernel_trace_stop_lvl2( flops );
    }
 
    PASTE_CORE_SLRMM_VOID;
    return flops;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform the addition of the low-rank matrix AB and the low-rank
 * matrix C.
 *
 *******************************************************************************
 *
 * @param[inout] params
 *          The LRMM structure that stores all the parameters used in the LRMM
 *          functions family.
 *          On exit, the C matrix is udpated with the addition of AB.
 *          @sa core_slrmm_t
 *
 * @param[in] AB
 *          The low-rank structure of the AB matrix to apply to C.
 *
 * @param[in] transV
 *          Specify if AB->v is stored normally or transposed.
 *          - If PastixNoTrans, AB->v is stored normally for low-rank format.
 *          - If PastixTrans, AB->v is stored transposed.
 *          - If PastixTrans, AB->v is stored transposed, and () must be
 *          applied to the matrix.
 *
 *******************************************************************************
 *
 * @return The number of flops required to perform the operation.
 *
 *******************************************************************************/
static inline pastix_fixdbl_t
core_slr2lr( core_slrmm_t           *params,
             const pastix_lrblock_t *AB,
             pastix_trans_t          transV )
{
    PASTE_CORE_SLRMM_PARAMS( params );
    pastix_int_t rklimit = core_get_rklimit( Cm, Cn );
    pastix_int_t rAB = ( AB->rk == -1 ) ? pastix_imin( M, N ) : AB->rk;
    pastix_int_t ldabu = M;
    pastix_int_t ldabv = (transV == PastixNoTrans) ? AB->rkmax : N;
    pastix_fixdbl_t total_flops = 0.;
    pastix_fixdbl_t flops = 0.;
 
    assert( (C->rk >= 0) && (C->rk <= C->rkmax) );
 
    /*
     * The rank is too big, we need to uncompress/compress C
     */
    if ( (C->rk + rAB) > rklimit )
    {
        float *Cfr, *Coff;
        int allocated = 0;
        if ( (Cfr = core_slrmm_getws( params, Cm * Cn )) == NULL ) {
            Cfr = malloc( Cm * Cn * sizeof(float) );
            allocated = 1;
        }
        Coff = Cfr + Cm * offy + offx;
 
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_uncompress );
        cblas_sgemm( CblasColMajor, CblasNoTrans, CblasNoTrans,
                     Cm, Cn, C->rk,
                     (sone),  C->u, Cm,
                                         C->v, C->rkmax,
                     (szero), Cfr,  Cm );
        flops = FLOPS_SGEMM( Cm, Cn, C->rk );
 
        /* Add A*B */
        if ( AB->rk == -1 ) {
            core_sgeadd( PastixNoTrans, M, N,
                         alpha, AB->u, M,
                         beta,  Coff,  Cm );
            flops += (2. * M * N);
        }
        else {
            cblas_sgemm( CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)transV,
                         M, N, AB->rk,
                         (alpha), AB->u, ldabu,
                                             AB->v, ldabv,
                         (beta),  Coff,  Cm );
            flops += FLOPS_SGEMM( M, N, AB->rk );
        }
        kernel_trace_stop_lvl2( flops );
        total_flops += flops;
 
        /* Try to recompress */
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_recompress );
        core_slrfree(C); // TODO: Can we give it directly to ge2lr as this
        flops = lowrank->core_ge2lr( lowrank->use_reltol, lowrank->tolerance, -1, Cm, Cn, Cfr, Cm, C );
        kernel_trace_stop_lvl2_rank( flops, C->rk );
        total_flops += flops;
 
        if (allocated) {
            free(Cfr);
        }
    }
    /*
     * The rank is not too large, we perform a low-rank update
     */
    else {
        total_flops += lowrank->core_rradd( lowrank, transV, &alpha,
                                            M,  N,  AB,
                                            Cm, Cn, C,
                                            offx, offy );
    }
 
    PASTE_CORE_SLRMM_VOID;
    return total_flops;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform the addition of the low-rank matrix AB into the null matrix C.
 *
 *******************************************************************************
 *
 * @param[inout] params
 *          The LRMM structure that stores all the parameters used in the LRMM
 *          functions family.
 *          On exit, the C matrix contains the product AB aligned with its own
 *          dimensions.
 *          @sa core_slrmm_t
 *
 * @param[in] AB
 *          The low-rank structure of the AB matrix to apply to C.
 *
 * @param[in] transV
 *          Specify if AB->v is stored normally or transposed.
 *          - If PastixNoTrans, AB->v is stored normally for low-rank format.
 *          - If PastixTrans, AB->v is stored transposed.
 *          - If PastixTrans, AB->v is stored transposed, and () must be
 *          applied to the matrix.
 *
 * @param[in] infomask
 *          Mask of informations returned by the core_sxx2lr() functions.
 *          If CORE_LRMM_ORTHOU is set, then AB.u is orthogonal, otherwise an
 *          orthogonalization step is added before adding it to C.
 *
 *******************************************************************************
 *
 * @return The number of flops required to perform the operation.
 *
 *******************************************************************************/
static inline pastix_fixdbl_t
core_slr2null(core_slrmm_t           *params,
              const pastix_lrblock_t *AB,
              pastix_trans_t          transV,
              int                     infomask )
{
    PASTE_CORE_SLRMM_PARAMS( params );
    pastix_int_t rklimit = core_get_rklimit( Cm, Cn );
    pastix_int_t ldabu = M;
    pastix_int_t ldabv = (transV == PastixNoTrans) ? AB->rkmax : N;
    pastix_fixdbl_t total_flops = 0.;
    pastix_fixdbl_t flops;
    int allocated = 0;
    PASTE_CORE_SLRMM_VOID;
 
    assert( C->rk == 0 );
 
    if ( AB->rk > rklimit ) {
        float *Cfr, *Coff;
        if ( (Cfr = core_slrmm_getws( params, Cm * Cn )) == NULL ) {
            Cfr = malloc( Cm * Cn * sizeof(float) );
            allocated = 1;
        }
        Coff = Cfr + Cm * offy + offx;
 
        /* Set to 0 if contribution smaller than C */
        if ( (M != Cm) || (N != Cn) ) {
            memset( Cfr, 0, Cm * Cn * sizeof(float) );
        }
 
        /* Uncompress the AB product into C */
        flops = FLOPS_SGEMM( M, N, AB->rk );
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_uncompress );
        cblas_sgemm( CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)transV,
                     M, N, AB->rk,
                     (alpha), AB->u, ldabu,
                                         AB->v, ldabv,
                     (beta),  Coff,  Cm );
        kernel_trace_stop_lvl2( flops );
        total_flops += flops;
 
        /* Try to recompress C */
        kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_recompress );
        flops = lowrank->core_ge2lr( lowrank->use_reltol, lowrank->tolerance, -1, Cm, Cn, Cfr, Cm, C );
        kernel_trace_stop_lvl2_rank( flops, C->rk );
        total_flops += flops;
 
        if ( allocated ) {
            free( work );
        }
    }
    else {
        /*
         * Let's chech that AB->u is orthogonal before copying it to C.u
         */
        int orthou = infomask & PASTIX_LRM3_ORTHOU;
        if ( !orthou ) {
            float *ABfr;
            pastix_lrblock_t    backup;
            int allocated = 0;
 
            kernel_trace_start_lvl2( PastixKernelLvl2_LR_add2C_orthou );
 
            if ( AB->rk > 0 ) {
                if ( (ABfr = core_slrmm_getws( params, M * N )) == NULL ) {
                    ABfr = malloc( M * N * sizeof(float) );
                    allocated = 1;
                }
 
                cblas_sgemm( CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)transV,
                             M, N, AB->rk,
                             (sone),  AB->u, ldabu,
                                                 AB->v, ldabv,
                             (szero), ABfr,  M );
                flops = FLOPS_SGEMM( M, N, AB->rk );
            }
            else {
                ABfr = AB->u;
                flops = 0.0;
            }
 
            flops += lowrank->core_ge2lr( lowrank->use_reltol, lowrank->tolerance, rklimit,
                                          M, N, ABfr, M, &backup );
 
            core_slrcpy( lowrank, PastixNoTrans, alpha,
                         M, N, &backup, Cm, Cn, C,
                         offx, offy );
 
            kernel_trace_stop_lvl2( flops );
            core_slrfree( &backup );
            total_flops += flops;
 
            if ( allocated ) {
                free( ABfr );
            }
        }
        /*
         * AB->u is orthogonal, we directly copy AB->u into C
         */
        else {
            core_slrcpy( lowrank, transV, alpha,
                         M, N, AB, Cm, Cn, C,
                         offx, offy );
        }
    }
 
    return total_flops;
}
 
/**
 *******************************************************************************
 *
 * @brief Perform the addition of two low-rank matrices.
 *
 *******************************************************************************
 *
 * @param[inout] params
 *          The LRMM structure that stores all the parameters used in the LRMM
 *          functions family.
 *          On exit, the C matrix contains the addition of C and A aligned with its own
 *          dimensions.
 *          @sa core_slrmm_t
 *
 * @param[in] A
 *          The low-rank structure of the A matrix to add to C.
 *
 * @param[in] transV
 *          Specify if A->v is stored normally or transposed.
 *          - If PastixNoTrans, AB->v is stored normally for low-rank format.
 *          - If PastixTrans, AB->v is stored transposed.
 *          - If PastixTrans, AB->v is stored transposed, and () must be
 *          applied to the matrix.
 *
 * @param[in] infomask
 *          Mask of informations returned by the core_sxx2lr() functions.
 *          If CORE_LRMM_ORTHOU is set, then A.u is orthogonal, otherwise an
 *          orthogonalization step is added before adding it to C.
 *
 *******************************************************************************
 *
 * @return The number of flops required to perform the operation.
 *
 *******************************************************************************/
pastix_fixdbl_t
core_slradd( core_slrmm_t           *params,
             const pastix_lrblock_t *A,
             pastix_trans_t          transV,
             int                     infomask )
{
    pastix_lrblock_t *C = params->C;
    pastix_fixdbl_t   flops = 0.;
 
    if ( A->rk != 0 ) {
        pastix_atomic_lock( params->lock );
        switch ( C->rk ) {
        case -1:
            /*
             * C became full rank
             */
            flops = core_slr2fr( params, A, transV );
            break;
 
        case 0:
            /*
             * C is still null
             */
            flops = core_slr2null( params, A, transV, infomask );
            break;
 
        default:
            /*
             * C is low-rank of rank k
             */
            flops = core_slr2lr( params, A, transV );
        }
        assert( C->rk <= C->rkmax);
        pastix_atomic_unlock( params->lock );
    }
 
    return flops;
}