d_spm_get_diagonal_zeros.c

Go to the documentation of this file.

/**
 *
 * @file d_spm_get_diagonal_zeros.c
 *
 * SParse Matrix package precision dependent routines to detect zero diagonal
 * entries.
 *
 * @copyright 2016-2026 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 1.3.0
 * @author Fatima-Zahra Salhi
 * @date 2026-07-09
 *
 * @generated from /builds/fpruvost/spm/src/z_spm_get_diagonal_zeros.c, normal z -> d, Fri Oct 23 13:04:59 2026
 *
 * @ingroup spm_dev_zeros
 * @{
 *
 **/
#include "common.h"
 
/**
 *******************************************************************************
 *
 * @ingroup spm_dev_zeros
 *
 * @brief Scan a CSC or CSR matrix for zero or missing diagonal entries.
 *
 *******************************************************************************
 *
 * @param[in] spm
 *          The sparse matrix, in SpmCSC or SpmCSR format.
 *
 * @param[inout] zeros_n
 *          On entry, the capacity of zeros_list. On exit, the number of zero
 *          or missing diagonal entries.
 *
 * @param[inout] zeros_list
 *          On entry, a pre-allocated array of size *zeros_n. On exit, the
 *          first *zeros_n entries contain the indices of the zero or missing
 *          diagonal entries. The array may be reallocated.
 *
 *******************************************************************************
 *
 * @retval SPM_SUCCESS
 *          On success.
 *
 * @retval SPM_ERR_OUTOFMEMORY
 *          If zeros_list could not be enlarged.
 *
 *******************************************************************************/
 
static inline int
d_spm_get_diagonal_zeros_csx( const spmatrix_t *spm,
                              spm_int_t        *zeros_n,
                              spm_int_t       **zeros_list )
{
    const double *values = spm->values;
    const spm_int_t       *colptr;
    const spm_int_t       *rowptr;
    spm_int_t              baseval = spmFindBase( spm );
    spm_int_t              size    = *zeros_n;
    spm_int_t              i, j;
 
    assert( (spm->fmttype == SpmCSC) ||
            (spm->fmttype == SpmCSR) );
    if ( spm->fmttype == SpmCSC ) {
        colptr = spm->colptr;
        rowptr = spm->rowptr;
    }
    else {
        colptr = spm->rowptr;
        rowptr = spm->colptr;
    }
 
    *zeros_n = 0;
    for ( i=0; i<spm->n; i++ )
    {
        int has_nonzero_diag = 0;
 
        spm_int_t begin = colptr[i] - baseval;
        spm_int_t end   = colptr[i + 1] - baseval;
        spm_int_t diag  = i + baseval;
 
        for ( j=begin; j<end; j++ )
        {
            if ( rowptr[j] == diag ) {
                has_nonzero_diag = (values[j] != (double)0.0);
                break;
            }
        }
 
        /*
         * If the diagonal entry is missing, has_nonzero_diag remains false.
         */
        if ( !has_nonzero_diag ) {
            if ( *zeros_n == size ) {
                spm_int_t  newsize = spm_imin( 2 * size, spm->n );
                spm_int_t *newlist;
 
                newlist = realloc( *zeros_list,
                                   newsize * sizeof( spm_int_t ) );
                if ( newlist == NULL ) {
                    return SPM_ERR_OUTOFMEMORY;
                }
                *zeros_list = newlist;
                size        = newsize;
            }
            (*zeros_list)[(*zeros_n)++] = diag;
        }
    }
 
    return SPM_SUCCESS;
}
 
/**
 *******************************************************************************
 *
 * @ingroup spm_dev_zeros
 *
 * @brief Detect unknowns with a zero or missing diagonal entry
 *
 * @param[in] spm
 *          The sparse matrix, in SpmCSC or SpmCSR format.
 *
 * @param[out] zeros_n
 *          On exit, the number of unknowns whose diagonal entry is zero
 *          or absent.
 *
 * @param[out] zeros_list
 *          On exit, pointer to an array containing the baseval-indexed
 *          indices of the affected unknowns. Set to NULL when zeros_n is 0
 *          or when the function fails. The caller is responsible for freeing
 *          this array.
 *
 *******************************************************************************
 *
 * @retval SPM_SUCCESS
 *          On success.
 *
 * @retval SPM_ERR_BADPARAMETER
 *          If spm->fmttype is neither SpmCSC nor SpmCSR.
 *
 * @retval SPM_ERR_OUTOFMEMORY
 *          If the allocation of zeros_list failed.
 *
 *******************************************************************************/
 
int
d_spmGetDiagonalZeros( const spmatrix_t *spm,
                       spm_int_t        *zeros_n,
                       spm_int_t       **zeros_list )
{
    spm_int_t *newlist;
    int        rc;
 
    *zeros_n    = 0;
    *zeros_list = NULL;
 
    if ( (spm->fmttype != SpmCSC) &&
         (spm->fmttype != SpmCSR) )
    {
        return SPM_ERR_BADPARAMETER;
    }
 
    if ( spm->n == 0 ) {
        return SPM_SUCCESS;
    }
 
    *zeros_n    = spm_imax( 1, spm->n / 20 );
    *zeros_list = malloc( (*zeros_n) * sizeof( spm_int_t ) );
    if ( *zeros_list == NULL ) {
        *zeros_n = 0;
        return SPM_ERR_OUTOFMEMORY;
    }
 
    switch (spm->fmttype) {
    case SpmCSC:
    case SpmCSR:
        rc = d_spm_get_diagonal_zeros_csx( spm, zeros_n, zeros_list );
        break;
    case SpmIJV:
        fprintf( stderr, "spmGetDiagonalZeros: Unsupported for IJV matrices\n");
        rc = SPM_ERR_NOTIMPLEMENTED;
        break;
    }
 
    if ( rc != SPM_SUCCESS ) {
        free( *zeros_list );
        *zeros_list = NULL;
        *zeros_n    = 0;
        return rc;
    }
 
    /* Resize zeros_list to the exact number of entries */
    if ( *zeros_n == 0 ) {
        free( *zeros_list );
        *zeros_list = NULL;
    }
    else {
        newlist = realloc( *zeros_list,
                           *zeros_n * sizeof( spm_int_t ) );
        if ( newlist == NULL ) {
            free( *zeros_list );
            *zeros_list = NULL;
            *zeros_n    = 0;
            return SPM_ERR_OUTOFMEMORY;
        }
        *zeros_list = newlist;
    }
 
    return SPM_SUCCESS;
}
 
/**
 * @}
 */