core_dpotrfsp.c

Go to the documentation of this file.

/**
 *
 * @file core_dpotrfsp.c
 *
 * PaStiX kernel routines for Cholesky factorization.
 *
 * @copyright 2011-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 6.3.0
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @author Xavier Lacoste
 * @author Gregoire Pichon
 * @author Nolan Bredel
 * @date 2023-01-13
 * @generated from /builds/solverstack/pastix/kernels/core_zpotrfsp.c, normal z -> d, Mon Aug 28 13:40:33 2023
 *
 **/
#include "common.h"
#include "cblas.h"
#include "blend/solver.h"
#include "pastix_dcores.h"
#include "kernels_trace.h"
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
#define MAXSIZEOFBLOCKS 64
static double done  =  1.0;
static double mdone = -1.0;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
 
/**
 *******************************************************************************
 *
 * @ingroup kernel_blas_lapack_null
 *
 * @brief Compute the sequential static pivoting Cholesky factorization of the
 * matrix n-by-n A = L * L^t .
 *
 *******************************************************************************
 *
 * @param[in] n
 *          The number of rows and columns of the matrix A.
 *
 * @param[inout] A
 *          The matrix A to factorize with Cholesky factorization. The matrix
 *          is of size lda -by- n.
 *
 * @param[in] lda
 *          The leading dimension of the matrix A.
 *
 * @param[inout] nbpivots
 *          Pointer to the number of piovting operations made during
 *          factorization. It is updated during this call
 *
 * @param[in] criterion
 *          Threshold use for static pivoting. If diagonal value is under this
 *          threshold, its value is replaced by the threshold and the number of
 *          pivots is incremented.
 *
 *******************************************************************************
 *
 * @warning This routine will fail if it discovers a null or negative value on
 *          the diagonal during factorization.
 *
 *******************************************************************************/
static inline void
core_dpotf2sp( pastix_int_t        n,
               double *A,
               pastix_int_t        lda,
               pastix_int_t       *nbpivots,
               double              criterion )
{
    pastix_int_t k;
    double *Akk = A;   /* A [k  ][k] */
    double *Amk = A+1; /* A [k+1][k] */
    double  alpha;
 
    for (k=0; k<n; k++){
        if ( fabs(*Akk) < criterion ) {
            (*Akk) = (double)criterion;
            (*nbpivots)++;
        }
 
        /* Symmetric matrices, so imaginary part should be 0 */
        if ( (*Akk) < 0.0 )
        {
            pastix_print_error( "Negative diagonal term\n" );
        }
 
        *Akk = sqrt(*Akk);
        alpha = 1.0 / (*Akk);
 
        /* Scale the diagonal to compute L((k+1):n,k) */
        cblas_dscal(n-k-1, ( alpha ), Amk, 1 );
 
        /* Move to next Akk */
        Akk += (lda+1);
 
        cblas_dsyr(CblasColMajor, CblasLower,
                   n-k-1, -1.0,
                   Amk, 1,
                   Akk, lda);
 
        /* Move to next Amk */
        Amk = Akk+1;
    }
}
 
/**
 *******************************************************************************
 *
 * @brief Compute the block static pivoting Cholesky factorization of the matrix
 * n-by-n A = L * L^t .
 *
 *******************************************************************************
 *
 * @param[in] n
 *          The number of rows and columns of the matrix A.
 *
 * @param[inout] A
 *          The matrix A to factorize with Cholesky factorization. The matrix
 *          is of size lda -by- n.
 *
 * @param[in] lda
 *          The leading dimension of the matrix A.
 *
 * @param[inout] nbpivots
 *          Pointer to the number of piovting operations made during
 *          factorization. It is updated during this call
 *
 * @param[in] criterion
 *          Threshold use for static pivoting. If diagonal value is under this
 *          threshold, its value is replaced by the threshold and the nu,ber of
 *          pivots is incremented.
 *
 *******************************************************************************
 *
 * @warning This routine will fail if it discovers a null or negative value on
 *          the diagonal during factorization.
 *
 *******************************************************************************/
void
core_dpotrfsp( pastix_int_t        n,
               double *A,
               pastix_int_t        lda,
               pastix_int_t       *nbpivots,
               double              criterion )
{
    pastix_int_t k, blocknbr, blocksize, matrixsize;
    double *tmp,*tmp1,*tmp2;
 
    /* diagonal supernode is divided into MAXSIZEOFBLOCK-by-MAXSIZEOFBLOCKS blocks */
    blocknbr = pastix_iceil( n, MAXSIZEOFBLOCKS );
 
    for (k=0; k<blocknbr; k++) {
 
        blocksize = pastix_imin(MAXSIZEOFBLOCKS, n-k*MAXSIZEOFBLOCKS);
        tmp  = A+(k*MAXSIZEOFBLOCKS)*(lda+1);      /* Lk,k     */
 
        /* Factorize the diagonal block Akk*/
        core_dpotf2sp(blocksize, tmp, lda, nbpivots, criterion);
 
        if ((k*MAXSIZEOFBLOCKS+blocksize) < n) {
 
            tmp1 = tmp  + blocksize;       /* Lk+1,k   */
            tmp2 = tmp1 + blocksize * lda; /* Lk+1,k+1 */
 
            matrixsize = n-(k*MAXSIZEOFBLOCKS+blocksize);
 
            /* Compute the column L(k+1:n,k) = (L(k,k)D(k,k))^{-1}A(k+1:n,k)    */
            /* 1) Compute A(k+1:n,k) = A(k+1:n,k)L(k,k)^{-T} = D(k,k)L(k+1:n,k) */
                        /* input: L(k,k) in tmp, A(k+1:n,k) in tmp1   */
                        /* output: A(k+1:n,k) in tmp1                 */
            cblas_dtrsm(CblasColMajor,
                        CblasRight, CblasLower,
                        CblasTrans, CblasNonUnit,
                        matrixsize, blocksize,
                        (done), tmp,  lda,
                                           tmp1, lda);
 
            /* Update Ak+1k+1 = Ak+1k+1 - Lk+1k * Lk+1kT */
            cblas_dsyrk(CblasColMajor, CblasLower, CblasNoTrans,
                        matrixsize, blocksize,
                        (double)mdone, tmp1, lda,
                        (double)done,  tmp2, lda);
        }
    }
}
 
/**
 *******************************************************************************
 *
 * @brief Compute the Cholesky factorization of the diagonal block in a panel.
 *
 * @warning This routine will fail if it discovers a null or negative value on
 *          the diagonal during factorization.
 *
 *******************************************************************************
 *
 * @param[in] solvmtx
 *          Solver Matrix structure of the problem
 *
 * @param[in] cblk
 *          Pointer to the structure representing the panel to factorize in the
 *          cblktab array.  Next column blok must be accessible through cblk[1].
 *
 * @param[inout] dataL
 *          The pointer to the correct representation of the lower part of the data.
 *          - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
 *          - pastix_lr_block if the block is compressed.
 *
 *******************************************************************************
 *
 * @return The number of static pivoting performed during the diagonal block
 *         factorization.
 *
 *******************************************************************************/
int
cpucblk_dpotrfsp1d_potrf( SolverMatrix *solvmtx,
                          SolverCblk   *cblk,
                          void         *dataL )
{
    pastix_int_t  ncols, stride;
    pastix_int_t  nbpivots = 0;
    pastix_fixdbl_t time, flops;
    double *L;
    pastix_lrblock_t *lrL;
    double criterion = solvmtx->diagthreshold;
 
    time = kernel_trace_start( PastixKernelPOTRF );
 
    ncols  = cblk->lcolnum - cblk->fcolnum + 1;
    stride = (cblk->cblktype & CBLK_LAYOUT_2D) ? ncols : cblk->stride;
 
    /* check if diagonal column block */
    assert( cblk->fcolnum == cblk->fblokptr->frownum );
    assert( cblk->lcolnum == cblk->fblokptr->lrownum );
 
    if ( cblk->cblktype & CBLK_COMPRESSED ) {
        /* dataL is a LRblock */
        lrL = (pastix_lrblock_t *)dataL;
        L   = lrL->u;
        stride = ncols;
 
        assert( lrL->rk == -1 );
        assert( stride == lrL->rkmax );
    } else {
        L = (double *)dataL;
    }
 
    /* Factorize diagonal block */
    flops = FLOPS_DPOTRF( ncols );
    kernel_trace_start_lvl2( PastixKernelLvl2POTRF );
    core_dpotrfsp(ncols, L, stride, &nbpivots, criterion );
    kernel_trace_stop_lvl2( flops );
 
    kernel_trace_stop( cblk->fblokptr->inlast, PastixKernelPOTRF, ncols, 0, 0, flops, time );
 
    if ( nbpivots ) {
        pastix_atomic_add_32b( &(solvmtx->nbpivots), nbpivots );
    }
    return nbpivots;
}
 
/**
 *******************************************************************************
 *
 * @brief Compute the Cholesky factorization of one panel.
 *
 *******************************************************************************
 *
 * @param[in] solvmtx
 *          Solver Matrix structure of the problem
 *
 * @param[in] cblk
 *          Pointer to the structure representing the panel to factorize in the
 *          cblktab array.  Next column blok must be accessible through cblk[1].
 *
 * @param[inout] L
 *          The pointer to the correct representation of lower part of the data.
 *          - coeftab if the block is in full rank. Must be of size cblk.stride -by- cblk.width.
 *          - pastix_lr_block if the block is compressed.
 *
 *******************************************************************************
 *
 * @return The number of static pivoting during factorization of the diagonal
 * block.
 *
 *******************************************************************************/
int
cpucblk_dpotrfsp1d_panel( SolverMatrix *solvmtx,
                          SolverCblk   *cblk,
                          void         *L )
{
    pastix_int_t nbpivots;
    nbpivots = cpucblk_dpotrfsp1d_potrf( solvmtx, cblk, L );
 
    cpucblk_dtrsmsp( PastixRight, PastixLower,
                     PastixTrans, PastixNonUnit,
                     cblk, L, L, &(solvmtx->lowrank) );
    return nbpivots;
}
 
 
/**
 *******************************************************************************
 *
 * @brief Perform the Cholesky factorization of a given panel and apply all its
 * updates.
 *
 *******************************************************************************
 *
 * @param[in] solvmtx
 *          Solver Matrix structure of the problem
 *
 * @param[in] cblk
 *          Pointer to the structure representing the panel to factorize in the
 *          cblktab array.  Next column blok must be accessible through cblk[1].
 *
 * @param[in] work
 *          Temporary memory buffer.
 *
 * @param[in] lwork
 *          Temporary workspace dimension.
 *
 *******************************************************************************
 *
 * @return The number of static pivoting during factorization of the diagonal
 * block.
 *
 *******************************************************************************/
int
cpucblk_dpotrfsp1d( SolverMatrix       *solvmtx,
                    SolverCblk         *cblk,
                    double *work,
                    pastix_int_t        lwork )
{
    void        *L = cblk_getdataL( cblk );
    SolverCblk  *fcblk;
    SolverBlok  *blok, *lblk;
    pastix_int_t nbpivots;
 
    nbpivots = cpucblk_dpotrfsp1d_panel( solvmtx, cblk, L );
 
    blok = cblk->fblokptr + 1; /* First off-diagonal block */
    lblk = cblk[1].fblokptr;   /* Next diagonal block      */
 
    /* If there are off-diagonal blocks, perform the updates */
    for( ; blok < lblk; blok++ )
    {
        fcblk = solvmtx->cblktab + blok->fcblknm;
 
        if ( fcblk->cblktype & CBLK_FANIN ) {
            cpucblk_dalloc( PastixLCoef, fcblk );
        }
 
        cpucblk_dgemmsp( PastixLCoef, PastixTrans,
                         cblk, blok, fcblk,
                         L, L, cblk_getdataL( fcblk ),
                         work, lwork, &(solvmtx->lowrank) );
 
        cpucblk_drelease_deps( PastixLCoef, solvmtx, cblk, fcblk );
   }
 
    return nbpivots;
}