z_spm.c

/**
 *
 * @file z_spm.c
 *
 * SParse Matrix package precision dependent routines.
 *
 * @copyright 2016-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 1.0.0
 * @author Mathieu Faverge
 * @date 2013-06-24
 *
 * @precisions normal z -> c d s p
 *
 **/
#include "common.h"
#include "z_spm.h"
#include <string.h>

/**
 *******************************************************************************
 *
 * @ingroup spm_dev_check
 *
 * @brief This routine sorts the subarray of edges of each vertex in a
 * centralized spm stored in CSC or CSR format.
 *
 * Nothing is performed if IJV  format is used.
 *
 * @warning This function should NOT be called if dof is greater than 1.
 *
 *******************************************************************************
 *
 * @param[inout] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the same sparse matrix with subarrays of edges sorted by
 *          ascending order.
 *
 *******************************************************************************/
void
z_spmSort( spmatrix_t *spm )
{
    spm_int_t       *colptr = spm->colptr;
    spm_int_t       *rowptr = spm->rowptr;
    spm_complex64_t *values = spm->values;
    void *sortptr[2];
    spm_int_t n = spm->n;
    spm_int_t i, size;
    (void)sortptr;

    if (spm->dof > 1){
        fprintf(stderr, "z_spmSort: Number of dof (%d) greater than one not supported\n", (int)spm->dof);
        exit(1);
    }

    /* Sort in place each subset */
    if ( spm->fmttype == SpmCSC ) {
        for (i=0; i<n; i++, colptr++)
        {
            size = colptr[1] - colptr[0];

#if defined(PRECISION_p)
            spmIntSort1Asc1( rowptr, size );
#else
            sortptr[0] = rowptr;
            sortptr[1] = values;
            z_spmIntFltSortAsc( sortptr, size );
#endif
            rowptr += size;
            values += size;
        }
    }
    else if ( spm->fmttype == SpmCSR ) {
        for (i=0; i<n; i++, rowptr++)
        {
            size = rowptr[1] - rowptr[0];

#if defined(PRECISION_p)
            spmIntSort1Asc1( rowptr, size );
#else
            sortptr[0] = colptr;
            sortptr[1] = values;
            z_spmIntFltSortAsc( sortptr, size );
#endif
            colptr += size;
            values += size;
        }
    }
}

/**
 *******************************************************************************
 *
 * @ingroup spm_dev_check
 *
 * @brief This routine merge the multiple entries in a sparse
 * matrix by suming their values together.
 *
 * The sparse matrix needs to be sorted  first (see z_spmSort()).
 *
 * @warning This function should NOT be called if dof is greater than 1.
 *
 *******************************************************************************
 *
 * @param[inout] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the reduced sparse matrix of multiple entries were present
 *          in it. The multiple values for a same vertex are sum up together.
 *
 ********************************************************************************
 *
 * @return The number of vertices that were merged
 *
 *******************************************************************************/
spm_int_t
z_spmMergeDuplicate( spmatrix_t *spm )
{
    spm_int_t       *colptr = spm->colptr;
    spm_int_t       *oldrow = spm->rowptr;
    spm_int_t       *newrow = spm->rowptr;
    spm_complex64_t *newval = spm->values;
    spm_complex64_t *oldval = spm->values;
    spm_int_t n       = spm->n;
    spm_int_t baseval = spm->colptr[0];
    spm_int_t dof2    = spm->dof * spm->dof;
    spm_int_t idx, i, j, d, size;
    spm_int_t merge = 0;

    if ( spm->fmttype == SpmCSC ) {
        idx = 0;
        for (i=0; i<n; i++, colptr++)
        {
            size = colptr[1] - colptr[0];
            for (j = 0; j < size;
                 j++, oldrow++, oldval+=dof2, newrow++, newval+=dof2, idx++)
            {
                if ( newrow != oldrow ) {
                    newrow[0] = oldrow[0];
#if !defined(PRECISION_p)
                    memcpy( newval, oldval, dof2 * sizeof(spm_complex64_t) );
#endif
                }

                while( ((j+1) < size) && (newrow[0] == oldrow[1]) ) {
                    j++; oldrow++; oldval+=dof2;
#if !defined(PRECISION_p)
                    /* Merge the two sets of values */
                    for (d = 0; d < dof2; d++) {
                        newval[d] += oldval[d];
                    }
#endif
                    merge++;
                }
            }
            assert( ((merge == 0) && ( colptr[1] == idx+baseval)) ||
                    ((merge != 0) && ( colptr[1] >  idx+baseval)) );

            colptr[1] = idx + baseval;
        }
        assert( ((merge == 0) && (spm->nnz         == idx)) ||
                ((merge != 0) && (spm->nnz - merge == idx)) );

        if (merge > 0) {
            spm->nnz = spm->nnz - merge;
            spm->gnnz = spm->nnz;

            newrow = malloc( spm->nnz * sizeof(spm_int_t) );
            memcpy( newrow, spm->rowptr, spm->nnz * sizeof(spm_int_t) );
            free(spm->rowptr);
            spm->rowptr = newrow;

#if !defined(PRECISION_p)
            newval = malloc( spm->nnz * dof2 * sizeof(spm_int_t) );
            memcpy( newval, spm->values, spm->nnz * dof2 * sizeof(spm_complex64_t) );
            free(spm->values);
            spm->values = newval;
#endif
        }
    }

    (void)d;
    return merge;
}

/**
 *******************************************************************************
 *
 * @ingroup spm_dev_check
 *
 * @brief This routine corrects the sparse matrix structure if it's
 * pattern is not symmetric.
 *
 * It returns the new symmetric pattern with zeroes on
 * the new entries.
 *
 *******************************************************************************
 *
 * @param[inout] spm
 *          On entry, the pointer to the sparse matrix structure.
 *          On exit, the same sparse matrix with extra entries that makes it
 *          pattern symmetric.
 *
 *******************************************************************************
 *
 * @retval Return the number of entries added to the matrix.
 *
 *******************************************************************************/
spm_int_t
z_spmSymmetrize( spmatrix_t *spm )
{
    spm_complex64_t *oldval, *valtmp, *newval = NULL;
    spm_int_t *oldcol, *coltmp, *newcol = NULL;
    spm_int_t *oldrow, *rowtmp, *newrow = NULL;
    spm_int_t *toaddtab = NULL;
    spm_int_t *toaddtmp, toaddcnt, toaddsze;
    spm_int_t  n       = spm->n;
    spm_int_t  dof2    = spm->dof * spm->dof;
    spm_int_t  i, j, k, r, size;
    spm_int_t  baseval;

    toaddcnt = 0;
    toaddsze = 0;

    if ( (spm->fmttype == SpmCSC) || (spm->fmttype == SpmCSR) ) {
        if (spm->fmttype == SpmCSC) {
            oldcol = spm->colptr;
            coltmp = spm->colptr;
            oldrow = spm->rowptr;
            rowtmp = spm->rowptr;
        }
        else {
            oldcol = spm->rowptr;
            coltmp = spm->rowptr;
            oldrow = spm->colptr;
            rowtmp = spm->colptr;
        }

        baseval  = oldcol[0];
        for (i=0; i<n; i++, coltmp++)
        {
            size = coltmp[1] - coltmp[0];
            for (r=0; r<size; r++, rowtmp++ )
            {
                j = rowtmp[0]-baseval;
                if ( i != j ) {
                    /* Look for the element (j, i) */
                    spm_int_t frow = oldcol[ j ] - baseval;
                    spm_int_t lrow = oldcol[ j+1 ] - baseval;
                    int found = 0;

                    for (k = frow; (k < lrow); k++)
                    {
                        if (i == (oldrow[k]-baseval))
                        {
                            found = 1;
                            break;
                        }
                        else if ( i < (oldrow[k]-baseval))
                        {
                            /* The spm is sorted so we will not find it later */
                            break;
                        }
                    }

                    if ( !found ) {
                        if ( newcol == NULL ) {
                            newcol = malloc( (spm->n+1) * sizeof(spm_int_t) );
                            for (k=0; k<spm->n; k++) {
                                newcol[k] = oldcol[k+1] - oldcol[k];
                            }

                            /* Let's start with a buffer that will contain 5% of extra elements */
                            toaddsze = spm_imax( 0.05 * (double)spm->nnz, 1 );
                            toaddtab = malloc( 2*toaddsze * sizeof(spm_int_t) );
                        }

                        if (toaddcnt >= toaddsze) {
                            toaddsze *= 2;
                            toaddtab = (spm_int_t*)realloc( toaddtab, 2*toaddsze*sizeof(spm_int_t) );
                        }

                        /* Newcol stores the number of elements per column */
                        newcol[ j ]++;
                        /* toaddtab stores the couples (j, i) to be added in the final sparse matrix */
                        toaddtab[ toaddcnt * 2     ] = j;
                        toaddtab[ toaddcnt * 2 + 1 ] = i;
                        toaddcnt++;
                    }
                }
            }
        }

        if (toaddcnt > 0) {
            spm_int_t newsize, oldsize;

            /* Sort the array per column */
            spmIntSort2Asc1(toaddtab, toaddcnt);

            spm->nnz = spm->nnz + toaddcnt;
            spm->gnnz = spm->nnz;

            newrow = malloc( spm->nnz        * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
            newval = malloc( spm->nnz * dof2 * sizeof(spm_complex64_t) );
#endif
            coltmp = newcol;
            rowtmp = newrow;
            valtmp = newval;
            oldval = spm->values;
            toaddtmp = toaddtab;

            newsize = coltmp[0];
            coltmp[0] = baseval;
            for (i=0; i<n; i++, coltmp++, oldcol++)
            {
                /* Copy the existing elements */
                oldsize = oldcol[1] - oldcol[0];
                memcpy( rowtmp, oldrow, oldsize * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
                memcpy( valtmp, oldval, oldsize * dof2 * sizeof(spm_complex64_t) );
#endif

                oldrow += oldsize;
                rowtmp += oldsize;
                oldval += oldsize * dof2;
                valtmp += oldsize * dof2;

                /* Some elements have been added */
                assert( newsize >= oldsize );
                if ( newsize > oldsize ) {
                    int added = 0;
                    int tobeadded = newsize - oldsize;

                    /* At least one element is equal to i */
                    assert( toaddtmp[0] == i );

                    /* Let's add the new vertices */
                    while( (added < tobeadded) && (toaddtmp[0] == i) ) {
                        rowtmp[0] = toaddtmp[1] + baseval;
                        rowtmp++; toaddtmp+=2; added++;
                    }
                    assert( added == tobeadded );

#if !defined(PRECISION_p)
                    /* Add 0.s for the associated values */
                    memset( valtmp, 0, added * dof2 * sizeof(spm_complex64_t) );
                    valtmp += added * dof2;
#endif
                }

                /* Use oldsize as temporary variable to update the new colptr */
                oldsize = newsize;
                newsize = coltmp[1];
                coltmp[1] = coltmp[0] + oldsize;
            }

            assert( coltmp[0]-baseval == spm->nnz );
            free( toaddtab );
            free( spm->colptr );
            free( spm->rowptr );
            free( spm->values );
            if (spm->fmttype == SpmCSC) {
                spm->colptr = newcol;
                spm->rowptr = newrow;
            }
            else {
                spm->colptr = newrow;
                spm->rowptr = newcol;
            }
            spm->values = newval;
        }
    }

    return toaddcnt;
}